Julien Derr
email: julien.derr@univ-paris-diderot.fr

Andrew D Rutenberg
email: andrew.rutenberg@dal.ca

Monodisperse domains by proteolytic control of the coarsening instability

Keywords: numbers: 87.16.A-, 87.16.dt, 64.75.Jk

The coarsening instability typically disrupts steady-state cluster-size distributions. We show that degradation coupled to the cluster size, such as arising from biological proteolysis, leads to a novel fixed-point cluster size. Stochastic evaporative and condensative fluxes determine the width of the fixed-point size distribution. At the fixed-point, we show how the peak size and width depend on number, interactions, and proteolytic rate. This proteolytic size-control mechanism is consistent with the phenomenology of pseudo-pilus length control in the general secretion pathway of bacteria.

I. INTRODUCTION

Living cells control the sizes of subcellular structures. Mechanisms of size control include molecular rulers for bacterial injectisome length [1], measuring cups for the length of the bacterial flagellar hook [2], counting for telomere ends [3], and equilibrium energetics for actin bundle radius [4]. Size control is also evident in the length of eukaryotic flagella [5] or the size of lipid rafts [6].

Size control is challenging in bacteria because of the strong stochastic effects expected in such small cells. It is especially interesting how bacteria control the size of extracellular macromolecular assemblies, such as bacterial secretion systems and pili. In this paper, we investigate a novel length-control mechanism that may apply to the pseudopilus (here "ppilus") of the type-II secretion system (T2SS) of Gram-negative bacteria [7]. In the general secretory pathway, proteins are first secreted across the bacterial inner membrane by the Sec or Tat systems then across the outer membrane by the T2SS. The T2SS pushes folded proteins across the periplasm, and out a secreton in the outer membrane, using an assembling and disassembling ppilus that is thought to function as a piston or plunger [7]. The ppilus assembles from the energized inner membrane and spans the periplasmic space. The primary pilin subunit that assembles into the pseudopilus of the T2SS is variously called PulG in Klebsiella oxytoca, XcpT in Pseudomonas aeruginosa, or more generally GspG (here "G") and is homologous to the PilA pilin of the type-IV pilus used in twitching motility in, e.g., P. aeruginosa or Myxococcus xanthus [8].

For a functioning T2SS the ppilus length should span the periplasm, which is approximately 21 nm across [9] -or 85 G monomers in the ppilus structure [10]. The ppilus is not normally seen outside the cell [7], hence the "pseudo" prefix, indicating an effective length-control mechanism. Overexpression of G leads to visible extracellular ppili [11][12][13] which rules out the (fixed-size) molecular ruler or measuring cup mechanisms of size regulation.

Underexpression of a minor pilin (XcpX [11], PulK [13], or GspK -here "K") leads to long extracellular ppili. So, it has been proposed that simple stoichiometric control applies to T2SS ppilus length control [14]. However, stoichiometric control faces the inherent challenges of precisely controlled protein expression [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF]. Furthermore, while a single "stoichiometric" ppilus with a fixed size pool of G monomers can have a narrow length distribution peaked at the pool size [START_REF]With a pool of Ng monomers and a ratio t of termination rate to monomer addition rate, the length distribution of a single ppilus is p1(l) = tΓ(Ng + 1)Γ(Ng + tl)/[Γ(Ng + t + 1)Γ(Ngl + 1)]. For small t this is peaked at Ng[END_REF], this is not true of multiple ppili sharing a common pool of G monomers. With multiple ppili stoichiometric length distributions are exponential, as seen in stochastic simulations with more than one ppili (see Appendix A). Since five to ten ppili are present on each individual bacteria [START_REF] Buddelmeijer | [END_REF], stoichiometric length control of each individual ppilus requires an additional mechanism to partition G proteins equally between the ppili.

II. MODEL

The observations of G-G interactions [18] and of Gclusters in individual bacteria [18,[START_REF] Scott | Proc. Natl. Acad. Sci[END_REF] are consistent with clusters of G that could be associated with each ppilus, as illustrated by the dashed regions in Fig. 1. The size of each cluster, i.e. the number of G monomers, would then determine the maximal length of the associated ppilus, converting a ppilus length-control problem into a G-cluster size-control problem. Nevertheless, thermally driven evaporation, condensation, and diffusion (Fig. 1 (e), (c), and (D) respectively) will destabilize spontaneous partitioning of G among many clusters. The subsequent coarsening of the size-distribution, treated by Lifshitz, Slyozov and Wagner (LSW) [START_REF] Bray | [END_REF], would lead to a (-)

(+) (D) (E) (C) OM IM
FIG. 1: (Color online) Cartoon of the bacterial periplasmic space showing the inner (IM) and outer (OM) membrane, with pilin "G" proteins (red circles). Two pseudo-pili are shown, extending from their inner membrane base towards their outer membrane secretons (green squares). The IM clusters associated with the ppili are indicated by the dashed circles. The number of G proteins in the ppilus together with its associated cluster is unchanged by ppilus assembly or disassembly (not shown). Numbers change by the physical processes of condensation (c), evaporation (e), proteolysis or recycling (-), and insertion (+). Diffusion (D) of monomers couples the clusters.

single large cluster (condensed phase) of G in equilibrium with small clusters of G (vapour) in the bacterial inner membrane. The net growth of large diffusively-coupled 2d clusters can be expressed in terms of their radius R by dR/dt = A (1/R c -1/R) /R, where A is a constant. Cluster growth is the net result of a condensation flux (due to the supersaturation of G monomers in the inner membrane) and an evaporative flux (due to the curved boundary of the cluster within the membrane). R c is the critical cluster size above which clusters grow and below which they shrink, and it typically grows with time due to the decreasing supersaturation associated with increasing average cluster radius via the Gibbs-Thompson effect. Growth of the average cluster size is also associated with a decreasing number of clusters, and with a broad distribution of cluster sizes [START_REF] Bray | [END_REF]. This is qualitatively unchanged for collections of smaller or irregularlyshaped clusters.

In a biological system we should add both protein synthesis and proteolysis. Protein synthesis will simply contribute to the supersaturation, but proteolysis will add a new degradative term to the LSW dynamics proportional to the number of monomers in the cluster. Bacterial proteolytic mechanisms include cytoplasmic ATP-dependent proteosomes [21], ubiquitin-like targeting systems [22], and periplasmic ATP-independent proteases [23]. Nondegradative recycling of components away from the membrane has an equivalent effect [24][25][26]. In terms of the number of monomers in a large 2d cluster of size N we will then have

dN dt = S - E √ N -αN, (1)
where S corresponds to condensation due to supersaturation, E corresponds to evaporation, and α is the proteolytic rate. Eqn. 1 is an approximate mean-field equation, without stochastic effects or higher-order curvature corrections expected for smaller clusters. However, it illustrates (see inset of Fig. 2) how intermediate values of the proteolytic term should generically stabilize coarsening to a steady-state cluster size (red dot in inset), while for large enough proteolytic rate no stable fixed-point exists (blue curve in inset).

To test these ideas we model G clustering in the bacterial membrane with a stochastic Ising lattice-gas in twodimensions (see e.g. [27]). The dynamics are (conserved) particle-exchange, subject to a Metropolis acceptance criterion with a reduced interaction energy J ≡ J/(k B T ). We supplement these dynamics with a dimensionless proteolysis rate α (per monomer per timestep), where (to minimize finite-size effects) monomers are removed and replaced at random positions in the system. Using a typical monomer size of ∆x = 5nm and diffusivity of D = 70000nm 2 /s = δx 2 /(4δt) we have ∆t ≈ 10 -4 s [27], and the proteolytic lifetime is τ = ∆t/α. In units of ∆x, we use a linear lattice size L = 400 and check that finite-size effects are not significant. A typical bacterium is larger in every direction, with L ≈ 2000 [27].

III. RESULTS

Figure 2 illustrates the steady-state results of the stochastic simulations. For a high proteolytic rate (α = 10 -5 , or protein lifetime τ ≈ 10s) no large domains are seen. At an intermediate proteolytic rate (α = 5 × 10 -8 , or τ ≈ 2000s), stable domains with a characteristic size are seen. At a low proteolytic rate (α = 10 -9 , or τ ≈ 10 5 s) proteolysis is insufficient to stabilize the largest domain against coarsening to the limits of the system size. For the remainder of the paper, and for intermediate proteolytic rates, we quantify the peak size N 0 as well as the full-width-half-maximum W of the steadystate distribution of domain sizes, and explore how they vary with proteolytic rate, total membrane density of G, and reduced interaction J.

The growth dynamics of a given cluster of size N is given by the difference of incoming and outgoing flux, Ṅ = Ṅ+ (N ) -Ṅ-(N ). This gives us the transition probabilities of monomer addition Γ + (N ) ≡ Ṅ+ /( Ṅ-+ Ṅ+ ) and subtraction Γ -(N ) ≡ Ṅ-/( Ṅ-+ Ṅ+ ). If P (N ) is the resulting steady state probability distribution of clusters of size N , then the detailed balance condition of the transition probabilities is Γ + (N -1)P (N -1) = Γ -(N )P (N ). Approximating P (N ) as a continuous distribution of cluster sizes (where dP/dN ≃ P (N ) -P (N -1)), then in steady state we have

P (N ) ≃ P (1) exp N 1 1 - Γ -(n) Γ + (n -1) dn , (2)
where we choose P (1) such that m P (m) = 1.

We have measured the steady-state evaporation and condensation rates N ± (N ), and used the resulting transition probabilities Γ ± (N ) in Eqn. Intermediate proteolytic rates stabilize a peaked steady-state cluster size distribution with multiple large clusters, as illustrated by the corresponding snapshot of the system. The inset is a stability diagram from Eqn. 1, illustrating the effects of low (green, top curve, with only an unstable fixed point), medium (red, center curve, with an additional stable fixed point as indicated), and high (blue, bottom curve, with no fixed points) proteolytic rates on the domain size dynamics dN/dt vs. domain size N .

measured P (N ). We find reasonable agreement (see e.g. the dotted lines in Fig. 2 or the crosses in Fig. 3), showing that the dispersity of domain sizes around the stable fixed point N 0 (where , the distribution becomes peaked at N 0 = 0. This is consistent with Eqn. 1 and the stability plots in the inset of Fig. 2. The width decreases with α as the net stabilizing flux near the stable point increases. As shown by the inset of Fig. 3(a), the relative width W/N 0 nevertheless increases with α. The fraction f of G-monomers that are found in larger clusters (under the second peak in Fig. 2) slightly decreases with α but remains a considerable fraction (more than half) of the total, which indicates that proteolysis can be an efficient size-control mechanism. Experimentally, it appears that at least 20% of G-monomers are in ppili [28].

N + (N 0 ) = N -(N 0 )) is
In contrast, the steady-state peak cluster size N 0 , the width W near the peak, and the fraction of G in the peak clusters stay roughly constant over a large range of average G densities, ρ ≡ N g /L 2 , as shown in Fig. 3(b). We understand this as an effective (nonequilibrium) coexistence between a fixed density of G monomers and excess G in clusters -increasing ρ simply increases the number of clusters without significantly changing their sizedistribution. (When approximately one cluster is seen in the system, for smaller ρ, finite-size effects do appear -decreasing both N 0 and W . This is beginning to be apparent at the lowest ρ in Fig. 3(b). Strong finite-size effects are seen in the condensed-cluster fraction f , but simply arise from an approximately constant vapour density as ρ varies.) The lack of strong dependence of the cluster size N 0 on the total membrane density ρ is advantageous in terms of robust control of cluster size in the face of stochastic protein expression.

At a fixed proteolytic rate and expression level, the effects of varying J/(k B T ) (i.e. G-G interactions) are shown in Fig. 3(c). Thermal evaporation decreases with increasing J -leading both to an increasing f (see inset of Fig. 3(c)), a decreasing effective supersaturation, and a smaller stable N 0 . The fractional width W/N 0 of the peak of the cluster size distribution is narrower for weaker interactions, but at the same time a smaller fraction of G-monomers are in clusters.

IV. DISCUSSION

We find that intermediate levels of proteolysis controls the natural coarsening instability of condensed clusters, and leads to steady-state clusters with a relatively narrow distribution of sizes. For an intermediate proteolytic rate α = 5 × 10 -8 (turnover time τ ∼ 2000s), we obtain a cluster size N 0 ≈ 80 -remarkably close to the 85 G required to assemble a ppilus that spans the periplasmic space [9,10]. The fractional width is then about 30%, consistent with the lack of extracellular ppili observed under normal conditions. For cluster size-control to effectively control ppilus length, we predict that G-clusters are associated with the secreton base of the T2SS. Mutational variation of G-monomers to affect either their proteolytic susceptibility (α, see e.g. Fig. 3(a)) or G-G interactions ( J, see Fig. 3(c)) should also affect the distribution of ppilus lengths. We believe stoichiometric mechanisms, via GspK [11][12][13], control the ppilus length from individual associated G-clusters, while our proteolytic size control mechanism ensures that multiple Gclusters remain approximately equally sized.

We have shown that a novel mean-field fixed-point in the cluster size distribution arises from proteolysis, while the dispersity of cluster sizes around the fixed point arises from stochastic growth and shrinkage of clusters. The noise associated with proteolysis is intrinsically multiplicative, in that proteolysis only targets existing G proteins. Our lattice-gas model naturally implements both proteolysis and the thermal evaporation and condensation of clusters in a membrane.

The result is a monodisperse cluster-size distribution with a non-zero peak size, which qualitatively differs from some earlier work on lipid raft sizes in membranes that only found distributions peaked at N 0 ≈ 0 [25,26]. We believe this is due to the approximate evaporation/condensation dynamics [25] or the additive noise [26] used in those works. In contrast, earlier coarsegrained models of ternary mixtures with recycling [24] also applied to lipid nanodomains -did recover a nonzero peak size, though did not include recycling noise. Our microscopic two-component model is simpler, and the mean-field flow we present shows how the results are expected to be generic for proteolysis or recycling.

Proteolysis is not just for cellular cleanup. Targeted degradation can adjust timescales and levels of transcription or translation (see, e.g. [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF]). We have shown how it can also be used to control cluster sizes within the cell. Proteolysis contributes a new "evaporative" term that limits coarsening with a novel mean-field fixed point for the cluster size (inset of Fig. 2). The same mechanism will qualitatively apply whether the proteolysis targets all proteins in a cluster (as in Eqn. 1), appropriate for cytoplasmic or periplasmic proteases -or whether it targets the cluster periphery as might be appropriate for membrane associated proteases. In both cases the loss term in Eqn. 1 will grow with N , and so will lead to a stable fixed point at some N 0 . Proteolysis provides a new size-control mechanism to cells. We expect that proteolysis or analogous degradation terms, such as recycling, are widely exploited to achieve monodisperse steady-state clusters in other biological systems. abilities are dependent only on the pool of monomers. The total number of monomers is constant, so that N g = n + M i=1 l i , where l i are the lengths of the pili and n the number of monomers remaining in the monomeric pool. Since the growth rate is proportional to n, while the termination rate t is fixed, then the growth and termination probabilities are p + (n) = n/(n+t) and p -(n) = t/(n+t), respectively, where p + + p -= 1.

For a single pilus, the probability of achieving a maximum length l can be written as the product of the individual probabilities:

p 1 (l) = p + (N g ) × p + (N g -1) × • • • × p + (N g -(l -1)) × p -(N g -l). This leads to p 1 (l) = tΓ(N g + 1)Γ(N g + t -l) Γ(N g + t + 1)Γ(N g -l + 1) (3)
which is peaked around N p for t 1, as shown with M = 1 in Figure 4. A single pilus can achieve length-control by assembling most of the available pool of monomers before disassembly is triggered.

For the more biologically appropriate case of multiple pili per cell, we can use a mean-field approximation. In this case, t becomes M t, and N g becomes M N p -M l where l is the average length of the pili, and then p

M (l) = p + (M N p -M l) × p + (M N p -M l -1) × • • • × p + (M N p - M l -(l -1)) × p -(M N p -M l -l).
The p + are constant at O(1/M ), so that in the limit M ≫ 1 we recover an exponential distribution

p ∞ (l) = t + 1 N p + t e -t+1 Np+t l (4)
For the intermediate regime, with finite M > 1, we performed stochastic computer simulations -as illustrated in Fig. 4. The exponential limit is quickly approached for M 4. However, even for M = 2 the peak around N p seen for M = 1 is lost.

Appendix B -Analytical monotonicities of the width

For the proteolytic size-control mechanism described in the text, we detail here the formal derivations of the variation of the full width half maximum (FWHM) W of the non-zero peak in the size distribution. We consider the evaporative and condensation fluxes J + and J -, respectively, close to the stable fixed point size N = N 0 , where J + = J -= J 0 (see Fig. 5). We allow for linear dependence of the fluxes near the fixed point J + = a(N -N 0 ) + J 0

(5) J -= b(N -N 0 ) + J 0 (6) where the stability of the fixed point requires J -> J + for N > N 0 and J + > J -for N < N 0 , corresponding to the requirement that b < a.

A monotonically decreases with n, since Γ + (n -1)/Γ -(n) = (a(n -1) + J 0 )/(b(n -1) + J 0 ) < 1 where b > a. The FWHM condition, A(W/2) = 1/2, allows us to determine how W must respond to changes in J 0 and ∆ ≡ b -a.

Varying J 0 at FWHM we have (10) where c/2 < b and a < b so that d log A/dJ 0 > 0 and A monotonically increases with J 0 . We conclude that W increases with increasing J 0 .

d log A dJ 0 = 1 cW/4 + J 0 - 1 bW/2 + J 0 + W/2-1 n=1 1 na + J 0 - 1 nb + J 0
In the same manner,

d log A d∆ = W W ∆ + 2aW + 4J 0 - W W ∆ + aW + 2J 0 - W/2-1 n=1 n n∆ + na + J 0 (11)
is negative, so that A monotonically decreases with ∆ ≡ b -a. We conclude that W decreases with increasing b -a.

FIG. 2 :

 2 FIG.2:(Color online) Cluster size probability distribution, P (N ), vs. number of monomers in a cluster, N , for low (α = 10 -9 , green, right), medium (α = 5 × 10 -8 , red, center), and high (α = 10 -5 , blue, left) proteolytic rates (per protein per timestep, where ∆t ≈ 10 -4 s). Other parameters are J = 1.67 and ρ = 2 × 10 -3 . The dotted lines represent the reconstruction of the probability distribution via Eqn. 2. Intermediate proteolytic rates stabilize a peaked steady-state cluster size distribution with multiple large clusters, as illustrated by the corresponding snapshot of the system. The inset is a stability diagram from Eqn. 1, illustrating the effects of low (green, top curve, with only an unstable fixed point), medium (red, center curve, with an additional stable fixed point as indicated), and high (blue, bottom curve, with no fixed points) proteolytic rates on the domain size dynamics dN/dt vs. domain size N .

  due to stochastic fluctuations driven by the finite fluxes. Intuitively, and as proven in Appendix B, the stochastic width W of the size-distribution increases either with increasing absolute flux at the fixed point (i.e. |N ± |) or with decreasing net (stabilizing) flux (|N + -N -|) near the stable fixed point. Fig. 3(a) shows the peak position and the width of the steady-state cluster-size distribution vs. proteolytic rate α. Increasing α decreases the stable domain size until, for α 4 × 10 -6

FIG. 3 :

 3 FIG. 3: (Color online) (a) Dependence of the steady-state most-likely cluster size N0 (red circles) and the full-width-halfmaximum W (green squares) of the cluster-size distribution vs. the proteolytic rate α. Crosses (×) indicate corresponding values derived from fluxes near N0 via Eqn. 2. The inset shows the relative width W/N0 (blue squares) and the fraction of G-proteins in large domains (purple circles) vs. α. Error bars indicate statistical errors, which are often smaller than point sizes. Other parameters are J = 1.67 and ρ = 2 × 10 -3 . (b) The top figure shows size N0 (red circles) and width W (green squares) vs. average membrane density ρ ≡ Ng/L 2 . The bottom figure shows the dependence of relative width W/N0 (blue squares) and the cluster fraction f (purple circles). Other parameters: J = 1.67 and α = 5 × 10 -8 . (c) Size N0 (red circles) and width W (green squares) vs the clustering interaction J ≡ J/(kBT ). The inset shows the relative width W/N0 (blue squares) and the fraction of G-proteins in large domains f (purple circles) vs. J. Other parameters: ρ = 2 × 10 -3 and α = 5 × 10 -8 .

FIG. 4 :

 4 FIG. 4: (Color online) Probability distribution pM (l) of the maximum length for various number of pili (M = 1, 2, 3, 4, 6) sharing the same common pool of monomers. We use t = 0.3 and Np = 40, though qualitatively similar results are seen with other parameters. The continuous lines are the result of the stochastic simulation. The red dashed lines represents the two exact limiting distributions p1 and p∞, overlaying the M = 1 data and close to the M = 6 data, respectively.

FIG. 5 :

 5 FIG.5: Figure illustrating the mean-field fixed-point cluster size N0 that occurs where the evaporative and condensative fluxes are equal (J-= J+ = J0). The fixed-point is stable when evaporation is stronger than condensation (J-> J+) for N > N0 -as illustrated. The demonstration of monotonicities in this appendix are for N ≈ N0, where a linear approximation for J± holds.

Acknowledgments

This work was supported by Natural Sciences and Engineering Research Council (NSERC), Canadian Institutes for Health Research (CIHR), and Atlantic Computational Excellence Network (ACENET); computational resources came from ACENET. J.D. was also supported by the Human Frontier Science Program. We thank Patrick McKelvey for stochastic simulations of the stoichiometric length-control mechanism. We acknowledge useful discussions with Olivera Francetic and Anthony Pugsley.

Appendix A -Stoichiometric length-control

We present here analytical limits as well as computational simulations for a stoichiometric model for lengthcontrol. This system consists of M polymerized pili sharing a common pool of N g = M N p monomers, where N p is the average number of monomers per pilus. Each pilus has a probability p + of growing and p -of terminating growth and completely disassembling, and these prob-


Cristóbal Bertoglio -Philippe

Moireau -Jean-Frederic Gerbeau

Cristóbal Bertoglio

Philippe Moireau

Jean-Frederic Gerbeau

Sequential parameter estimation for fluid-structure problems. Application to hemodynamics

Keywords: Observation, Modeling, and Inverse problem, fluid-structure interaction, reduced order Unscented Kalman Filter, arterial stiffness estimation

. It does not require any adjoint or tangent problems. In addition, it can easily be run in parallel, which is of great interest in fluid-structure problems where the computational cost of the forward simulation is already a challenge in itself. We illustrate our methodology with the estimation of the artery wall stiffness from the wall displacement measurements -as they could be extracted from medical imaging -in a three-dimensional idealized abdominal aortic aneurysm. We also show preliminary results about the estimation of the proximal Windkessel resistance, which is an important parameter for setting appropriate fluid boundary conditions.

Estimation de paramètres séquentielle pour des problèmes fluide-structure. Application à l'hémodynamique Résumé : Nous proposons une stratégie d'estimation de paramètres robuste et efficace pour des problèmes d'interaction fluide-structure. La méthode est basée sur un algorithme de filtrage limité à l'espace des paramètres, connue sous le nom de Filtre de Kalman réduit [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF]. Il ne nécessite aucun adjoint ou problème tangent. De plus, il peut facilement être exécuté en parallèle, ce qui est d'un grand intérêt pour les problèmes d'interaction fluide-structure, où le coût de calcul de la simulation directe est déjà un défi en soi. Nous illustrons notre méthodologie avec l'estimation de la rigidité des artères à partir de mesures de déplacements de la paroi -tels qu'ils pourraient être extraits de l'imagerie médicale -dans un anévrisme de l'aorte abdominale idéalisée en trois dimensions. Nous montrons aussi des résultats préliminaires pour l'estimation de la résistance proximale de Windkessel, qui est un paramètre important pour les conditions aux limites du fluide.

Mots-clés : Problèmes inverses, interaction fluide-structure, Unscented Kalman Filter, estimation de la raideur des artères

Introduction

During the last decade, many works have been devoted to fluid-structure interaction (FSI) in large arteries. Important progress has been made to improve the robustness and the efficiency of the coupling algorithms (see for example [START_REF] Fernández | Fluid structure interaction problems in haemodynamics, chapter 9[END_REF] and the references therein). But to make the FSI models useful in the clinical context, they have to be personalized, i.e. adapted to the patient. This requires not only to get the patient anatomy, but also to fit the model parameters using the available clinical measurements. The problem of reducing the uncertainties of a model using partial observations of the state variables is well-known in many fields of science and engineering and is referred to as data assimilation.

The purpose of this article is to present a data assimilation strategy for FSI problems. We apply it on an academical -but realistic -test case of blood flow to demonstrate its capability to recover the artery wall stiffness from the displacements of the wall.

Data assimilation of distributed mechanical systems can be based on a variational approach (see e.g. [START_REF] Banks | Estimation techniques for distributed parameter systems, volume 1 of Systems & Control: Foundations & Applications[END_REF]5]) by minimizing a least square criterion which includes a regularization term and the difference between the observations and the model prediction. One of the main difficulties of this approach lies in the iterative evaluation of the criterion, involving many solutions of the forward problem, and its gradient -typically adjoint-based -which usually requires a laborious implementation. In FSI for hemodynamics, this approach has been investigated for one-dimensional models in [START_REF] Lagrée | An inverse technique to deduce the elasticity of a large artery[END_REF][START_REF] Martin | Parameter identification for a one-dimensional blood flow model[END_REF], and for three-dimensional problems in [START_REF] D'elia | Applications of variational data assimilation in computational hemodynamics[END_REF] by modifying the minimization problem in order to avoid the resolution of the adjoint equations.

In this work, we consider a sequential approach, based on a generalization of the Kalman filter. With sequential algorithms, the model prediction is improved at every time instant by analyzing the discrepancy between the actual measurements and the model observation outputs [START_REF] Bensoussan | Filtrage optimal des systèmes linéaires[END_REF][START_REF] Simon | Optimal state estimation: Kalman, H infinity, and nonlinear approaches[END_REF]. In fact, in the fully linear case, it can be proved that the Kalman filter gives the same result as a variational approach based on a least squares criterion.

Even though full Kalman filter is not tractable for distributed systems, some effective sequential procedures for mechanical systems have been recently introduced: for the state estimation [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: Perspective in medical imaging[END_REF], and for the parameter identification [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF][START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF]. In particular, in [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF], the authors formulate a reduced-order version of the Unscented Kalman Filter (UKF) [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF] inspired from the Singular Evolutive Extended Kalman filter (SEIK) proposed in [START_REF] Pham | A singular evolutive extended kalman filter for data assimilation in oceanography[END_REF] and apply it for characterizing solid mechanics constitutive parameters. Here, we investigate the derivation of this filter in order to estimate the parameters of a FSI system. With respect to other approaches, this algorithm does not need adjoint or tangent operators and can be very easily parallelized. As a result, the total CPU time needed for the data assimilation is of the same order of magnitude as the CPU time needed for one forward simulation and the implementation does not need important modifications in existing solvers. Preliminary results of the present work have been presented in [4].

The rest of this paper is organized as follows. Section 2 introduces the FSI equations and the data assimilation problem. Section 3 is devoted to an introduction to data assimilation techniques, geared to readers with no background in this field. In particular, a short derivation of the standard Kalman filter and some extensions to nonlinear cases are presented. The purpose of these sections is to provide the reader with a comprehensive overview of the reduced-order Uncented Kalman filter whereas Section 4 points out the specific issues encountered in the derivation of the filtering strategy in the fluid-structure context. Finally, numerical investigations are presented in Section 5.

Problem setting

Fluid-Structure model

We consider the mechanical interaction between an incompressible fluid and an elastic structure. The fluid is described by the Navier-Stokes equations, in a moving domain Ω f (t) ⊂ R d , d = 2, 3, in an Arbitrary Lagrangian Eulerian (ALE) formulation, and the structure by the linear elasticity equations in

Ω s (t) ⊂ R d . The fluid-structure interface is denoted by Σ def = ∂Ω s ∩ ∂Ω f and ∂Ω f = Γ in ∪ Γ out ∪ Σ, ∂Ω s = Γ d ∪ Γ n ∪ Σ,
are given partitions of the fluid and solid boundaries, respectively (see Figure 1). The coupled FSI problem reads as follows: find the fluid velocity u(t

) : Ω f (t) → R d , the fluid pressure p(t) : Ω f (t) → R and the structure displacement y(t) : Ω s → R d such that • Fluid equations:                ρ f ∂u ∂t |b x + ρ f (u -w f ) • ∇u -2µ∇ • ǫ(u) + ∇p = 0, in Ω f (t), ∇ • u = 0, in Ω f (t), u = u in , on Γ in , σ(u, p)n = -pn, on Γ out . (1)
• Structure equations:

           ρ s ∂ t w -η s ∇ • σ s ( w) -∇ • σ s ( y) = 0, in Ω s , w = ∂ t y, in Ω s , y = 0, on Γ d , σ s ( y)n s = 0, on Γ n , (2)
• Coupling conditions:

             y f = Ext y | b Σ , w f = ∂ t y f , Ω f (t) = I b Ω f + y f A ( Ω f , t), u = w f , on Σ(t), σ s ( y) n s + J f σ(u, p) F -T f n f = 0, sur Σ , (3)
with F f the deformation gradient and J f = det F f . Note that the equalities u = w f on Σ(t) and y f = Ext y | b Σ could be limited to the normal components (whereas the continuity of all the components of the fluid and solid velocities, u and w, has of course to be enforced on Σ(t)). This problem has to be completed with appropriate initial conditions: velocity u 0 and domain displacement y 0 f for the fluid, initial velocity w 0 and displacement y 0 for the solid.

Ω s (t) Σ(t) Σ Ω s Ω f Ω f (t) Γ in A(•,t) Γ in Γ out Γ out Γ n Γ d
In the hemodynamics problems considered in the work, the outlet pressure p is obtained by solving an ordinary differential equation coupled to the threedimensional problem:

             C dπ dt + π -P V R d = Q π| t=0 = π 0 p = π + R p Q Q = Γ out u • n f , (4)
where the venous pressure P V , the distal resistance R d , the proximal resistance R p and the capacitance C are assumed to be given. This "zero-dimensional" equation is known as the three-element Windkessel model (see for example [START_REF] Formaggia | Cardiovascular Mathematics. Modeling and simulation of the circulatory system[END_REF]). Note that the pressure π belongs to the state variables of the problem. In realistic configurations, several outlets have to be considered. In such a case, there are as many systems (4) as outlets (see for example [START_REF] Moireau | External tissue support and fluid-structure simulation in blood flows[END_REF] for FSI in a patient-specific aorta).

Data assimilation problem

A physical system like blood flowing in a compliant artery can be observed through various measurement modalities: artery wall movements obtained from 4D medical imaging -computed tomography (CT), magnetic resonance imaging (MRI) or ultrasound (US) -cross section blood flow rates by Phase Contrast MRI or US, pressure in few points provided by a catheter, etc. These measurements are usually limited to a few locations in space, and sometimes come from different cardiac cycles (in this case they are resynchronized with the electrocardiogram). They are of course subject to noise and their postprocessing can introduce some further inaccuracies. Moreover, only a limited number of physical quantities can be simultaneously obtained.

The physical system of interest can also be modeled by equations like those presented in the previous section. In that case, many quantities -displacement, velocity, pressure, stress -are available in "all" locations and at "any" time instant. But of course, the model itself contains approximations -due for example to the modeling choices and the numerical solution -and relies on parameters that are not perfectly known.

The objective of data assimilation is to take advantage of both measurements and models. Measurements can be used to reduce the uncertainties of the model, and the model can be used to access some "hidden" physical quantities, for example the mechanical stress in the artery wall. Data assimilation can also be viewed as a way to reduce the measurements noise by means of a model that takes into account the underlying physical principles.

The real state variables of the physical system of interest (for example displacements, velocities, etc.) are denoted by x(t), and the parameters of the model by θ ∈ R p . From now, it is assumed that the problem has already been discretized in space and we denote the semi-discrete state X(t) ∈ R N . The state variable is solution of the dynamical system:

Ẋ = A(X, θ) , X(0) = X 0 . (5)
In this study X(t) typically denotes the finite element approximation of the displacements and the velocities in the fluid and the solid, and θ denotes for example the Young's modulus in different regions of the artery wall. So the orders of magnitude of N and p are typically a few hundreds of thousands and a dozen respectively. A discontinuous representation of the parameters is chosen here for the sake of simplicity but higher order discretization can also be considered [START_REF] Xi | Myocardial transversely isotropic material parameter estimation from in-silico measurements based on reduced-order unscented kalman filter[END_REF].

The measurements Z(t) are defined by an observation operator H(t) applied to the real state x(t). Assuming that an approximation X n of x(n∆t) and an approximation H n of H(n∆t) are available, the discrete measurement vector Z n ∈ R m -supposed to be known -is given by:

Z n = H n (X n ) + ζ Z n ,
where ζ Z n represents the noise of the measurement device as well as the inaccuracy resulting from the discretizations. These measurements are assumed to be available at every simulation time step. When this is not the case, they are given by a linear interpolation in time.

Data assimilation techniques usually consists of minimizing a cost function like:

J(X, θ) = T 0 Z -H(X) 2 W -1 dt + θ -θ0 2 (P θ 0 ) -1 + X(0) -X0 2 (P X 0 ) -1 ,
with X satisfying (5). In this expression, X0 and θ0 are given a priori values for the initial condition and parameters, and

• W -1 , • (P X 0 ) -1 and • (P θ 0 ) -1
denote some norms used to measure the observations, the state and the parameters, respectively. These norms allow to give a different weight to the different terms and therefore account for the "confidence" in the different quantities. From a statistical viewpoint, the "confidence" can be viewed as the inverse of covariance matrices (W , P X 0 and P θ 0 ), which explains the notation. Defining ζ X = X 0 -X0 , ζ θ = θ-θ0 , the optimization problem reads (keeping the same notation J for the cost function):

J(ζ X , ζ θ ) = T 0 Z -H(X) 2 W -1 dt + ζ θ 2 (P θ 0 ) -1 + ζ X 2 (P X 0 ) -1 . (6)

RR n°7657

This minimization problem can be addressed by many techniques that are classically divided in two groups: the variational and the sequential approaches. The variational approach consists of minimizing this cost function by an optimization algorithm that is usually based on the computation of its gradient, obtained by solving an adjoint model. The sequential approach, also known as filtering, modifies the forward dynamics with a correction term that takes into account the discrepancy between actual measurements and observations generated by the model:

Ẋ = A( X, θ) + K(Z -H( X)) (7)
where X is called estimator of X. The quantity Z -H( X) is known as the innovation, and the operator K depends on the method. For linear problems, the most famous sequential approach is the Kalman filter. On a given time interval [0, T ], and if all the operators are linear, the variational method and the Kalman filter algorithm turn out to give the same estimation at t = T . Concerning the computational complexity, whereas the variational method has to solve several forward and adjoint problems on the whole interval [0, T ], the estimation in the sequential algorithm is computed by solving only once the filtered dynamics [START_REF] Evensen | Data assimilation: the ensemble Kalman filter[END_REF]. However, the optimal operator K is determined by operations (multiplications, inversions, etc.) involving full matrices of the size of the state and the observations, which makes Kalman-based filters prohibitive for discrete problems derived from partial differential equations (PDEs).

Fortunately, several alternatives exist to reduce the computational cost of filtering approaches. For instance, in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] the authors have designed a Luenberger observer [START_REF] Luenberger | An introduction to observers[END_REF] for the solid mechanics state variables in order to obtain a state filter operator that is tractable. Then, the remaining uncertainties -mostly the parameters in their context -can be handled by the use of a reduced-order Kalman filter on the system corrected by the Luenberger filter [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF][START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF]. Our purpose in this paper is to restrict our investigation to that last part by deriving the reduced-order Uncented Kalman Filter (UKF) proposed in [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] for FSI problems. In this context, if the usual state variables arising from the PDEs discretization are assumed to be stable, the remaining uncertainties are the parameters. Nevertheless, in Section 5, it will be shown that the Windkessel pressure, which is a state variable, can easily and advantageously be added to the set of filtered quantities.

In the next section, we give an informal introduction to some sequential methods. We start with the least squares estimation in a static case, and we gradually enrich the method until the reduced UKF algorithm that will be later applied to three-dimensional FSI inverse problem.

Kalman filtering in a nutshell

The material presented in this section is standard in many fields, like weather forecasting or oceanography. Since it is less common in mechanics, we thought it might be useful to present it in order to make the article as self-contained as possible. For a more detailed review, we refer the reader to [START_REF] Kailath | Linear estimation[END_REF][START_REF] Simon | Optimal state estimation: Kalman, H infinity, and nonlinear approaches[END_REF] among many others.

Static linear case: least squares estimation

To start with, assume that there is no dynamics, that the parameters are perfectly known and that the observation operator is linear. Assume that a "guess" Xis available for the state, for example provided by a model, associated with a confidence matrix (P -) -1 . Assume also that an observation Z = HX + ζ Z is available. Since ζ Z is unknown, we also associate a confidence matrix W -1 with the observations.

A quantity accounting for Xand Z, with their respective levels of confidence, can simply be obtained by minimizing the quadratic cost function:

J( X) = 1 2 ( X -X-) T (P -) -1 ( X -X-) + 1 2 (Z -H X) T W -1 (Z -H X).
An elementary computation shows that the minimizer of J is given by

X+ = X-+ K (Z -H X-) , (8)
where

K = P + H T W -1 (9)
is the Kalman matrix and

P + = (P -) -1 + H T W -1 H -1 .
Applying the Woodbury matrix inversion lemma to P + , the Kalman matrix can also be written as

K = P XZ P Z -1 , (10)
with

P XZ = P -H T and P Z = W + HP -H T . (11)
This result can also be obtained with statistical arguments considering the unknown state X as a random variable. Assuming that X-= E(X) and that the error on the guess has a covariance matrix P -= E((X -X-)(X -X-) T ) (which is the inverse of the above "confidence matrix"). Then, looking for an X such that X = E(X) and such that the trace of the associated covariance matrix is minimized, it can be shown that the result X+ is again given by [START_REF] Fernández | Coupling schemes for incompressible fluid-structure interaction: implicit, semi-implicit and explicit[END_REF]. In what follows, we will adopt this statistical point of view.

It will be useful to note that the covariance of X+ is given by P + = (I -KH)P -or

P + = P --P XZ P Z -1 P XZ T . ( 12
)

Dynamic linear case: the Kalman filter

Now assume that the state variable satisfies a linear dynamics, without any model uncertainties:

Ẋ = AX + F,
and with an uncertain initial guess X(0) = X0 + ζ X ( X0 being given and ζ X unknown). After time discretization, the dynamics reads:

X n+1 = A n+1 X n + F n+1 . (13)
Suppose that the following observations are available:

Z n = H n X n + ζ Z n , (14)

RR n°7657

where ζ Z n includes the measurement noise and the discretization error. The linear Kalman filter [START_REF] Kalman | New results in linear filtering and prediction theory[END_REF] can simply be presented as follows. Assume that X+ n is known with a covariance P + n . First, the model is used to compute a prediction:

X- n+1 = A n+1 X+ n + F n+1 , (15)
By linearity of A n+1 , the covariance

P - n+1 = E((X n+1 -X- n+1 )(X n+1 -X- n+1 ) T ) is given by P - n+1 = A n+1 P + n A T n+1 . (16)
Then, repeating mutatis mutandis the least squares argument of the case without dynamics, a correction taking into account the observation is given by:

X+ n+1 = X- n+1 + K n+1 (Z n+1 -H n+1 X- n+1 ). (17)
with

K n+1 = P XZ n+1 (P Z n+1 ) -1 (18)
where

P XZ n+1 = P - n+1 H T n+1 and P Z n+1 = W n+1 + H n+1 P - n+1 H T n+1 .
As in the static case, the covariance associated with X+ n+1 is given by

P + n+1 = P - n+1 -P XZ n+1 P Z n+1 -1 P XZ n+1 T . (19)
Thus, we can iterate and perform the prediction step at the new time step n + 1.

It can be proved that these estimations of the state and covariances can also be obtained by minimizing the functional

J(ζ X ) = 1 2 ζ 2 (P - 0 ) -1 + 1 2 n+1 k=0 Z k -HX k 2 W -1 k ( 20
)
under the constraint [START_REF] Hoteit | A simplified reduced order Kalman filtering and application to altimetric data assimilation in Tropical Pacific[END_REF]. In other words, the Kalman filter estimator X+ n+1 is equal to X n+1 ( ζX ), where ζX is the minimizer of the variational problem on [0, t n+1 ].

Extensions to nonlinear cases

The Kalman filter presented above is only valid when the dynamics and the observation operator are linear. To address the fluid-structure problems considered in this work, it is necessary to extend the algorithm to nonlinear cases.

The most straightforward extension consists of applying steps ( 16), ( 18) and ( 19) with the tangent operators

A ′ n+1 def = A ′ n+1 ( X+ n ) and H ′ n+1 def = H ′ n+1 ( X- n+1
), and steps [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF] and [START_REF] Kailath | Linear estimation[END_REF] with the nonlinear operators A n+1 (•) and H n+1 (•):

               X- n+1 = A n+1 ( X+ n ) P - n+1 = A ′ n+1 P + n (A ′ ) T P XZ n+1 = P - n+1 (H ′ ) T , P Z n+1 = W n+1 + H ′ P - n+1 H ′ n+1 T X+ n+1 = X- n+1 + K n+1 (Z n+1 -H n+1 ( X- n+1 )) , K n+1 = P XZ n+1 (P Z n+1 ) -1 P + n+1 = P - n+1 -P XZ n+1 P Z n+1 -1 P XZ n+1 T (21)

RR n°7657

But the resulting algorithm, called Extended Kalman Filter (EKF), has two drawbacks: the computation of the tangent operators and the precision of the estimated values.

To understand the latter issue, it is interesting to consider how the mean value of a random variable is transported by a nonlinear operator A. Let X denote a random variable in R N , its mean value X = E(X) and its covariance P = E((X -X)(X -X) T ). A Taylor expansion around X gives:

Y def = A(X) = A( X)+A ′ ( X)•(X-X)+ 1 2 (X-X) T •A ′′ ( X)•(X-X)+O(|X-X| 3 )
Hence, the propagated mean is given by

E(Y ) = A( X) + 1 2 A ′′ ( X) : P + O(E(|X -X| 3 )) , (22)
and the propagated covariance by

E (Y -E(Y )) (Y -E(Y ) T = A ′ ( X) P A ′ ( X) T + O(E(|X -X| 3 )) . ( 23
)
We see in [START_REF] Moireau | Erratum of article "reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] that the simple propagation of the mean value by the nonlinear operator is only locally second order accurate (the term A ′′ ( X) : P being second order in |X -X|). This will be improved with the technique presented in section 3.4. Another way to approximate means and covariances consists of choosing a set of vectors X (i) , 1 ≤ i ≤ r, called particles, which are propagated with the operator A(•), i.e.,

Y (i) = A(X (i) ) , 1 ≤ i ≤ r .
Then, the mean and covariances are obtained by computing the empirical statistics of the propagated particles Y (i) . Of course, to be effective, the particles have to be carefully chosen. In this work, we use this type of mean and covariance approximations in the context of the Unscented Kalman Filter [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF], [START_REF] Julier | A new method for the nonlinear transformation of means and covariances in filters and estimators[END_REF], which is based on deterministic particles. It will be briefly presented in the next section.

The UKF filter presents a lot of similarities with the Ensemble Kalman Filter [START_REF] Evensen | Data assimilation: the ensemble Kalman filter[END_REF] which also uses a moderate number of particles -even if this number is often a larger than for UKF -to compute empirical statistics. Moreover nonlinear particle filters consider a random generation of the particles X (i) according to a given probability density function. These filters usually involve a very large number of particles, which yields a very large number of solutions of the forward problem. We refer for example to [START_REF] Simon | Kalman filtering with state constraints: A survey of linear and nonlinear algorithms[END_REF].

Unscented transforms and application to Kalman filtering

The Unscented Kalman Filter (UKF) algorithm is based on the unscented transformation. Its principle is to approximate the nonlinear propagation of the mean and the covariance of a random vector via the propagation of well-chosen deterministic particles.

The basic idea can easily be explained in one dimension. Let X be a random variable in R (so N = 1), let X denote its mean value E(X) and σ = Var(X) its standard deviation. Two particles are defined: X(1) = X + σ, X(2) = Xσ.

These particles are constructed by adding to the known state X two values σ and -σ called sigma-points. The empirical mean value of the particles, namely ( X(1) + X(2) )/2, is by construction equals to the mean value X. In order to evaluate the empirical mean value of the particles propagated by the nonlinear operator, let us consider the Taylor expansions:

A( X(1) ) = A( X) + σA ′ ( X) + σ 2 2 A ′′ ( X) + O(| X(1) -X| 3 ) , A( X(2) ) = A( X) -σA ′ ( X) + σ 2 2 A ′′ ( X) + O(|X (2) -X| 3 ) .
The empirical mean is by construction

Ŷ def = 1 2 A(X (1) ) + A(X (2) ) = A( X)+ σ 2 2 A ′′ ( X)+O(| X(1) -X| 3 )+O(| X(2) -X| 3 ).
(24) Comparing ( 22) and ( 24), we see that the empirical mean of the propagated particle is an approximation of E(A(X)) better than A( X) since it includes by construction the second order term σ 2 A ′′ ( X)/2. Note that the covariance are approximated with the same order of accuracy in the two methods.

This idea can be generalized in N -dimensions by defining appropriate sigmapoints and their respective weights -see [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] for a comprehensive review. The natural generalization of the one-dimensional standard deviation is the square root of the covariance matrix. Its computation, for example by a Cholesky factorization, is the most expensive part of the UKF algorithm, but usually negligible compared to the propagation step for the problems considered in this work.

The empirical mean and covariance are defined by:

X def = E α X ( * ) def = 1≤i≤r α i X (i) , Cov α X ( * ) , Z ( * ) def = 1≤i≤r α i (X (i) -X)(Z (i) -Z) T , (25)
where the weights α i ∈ R + are given and depend on the choice of the sigmapoints. In the present work, we use the ones proposed in [START_REF] Julier | The spherical simplex unscented transformation[END_REF], only based on r = N + 1 particles (see Appendix A).

It remains to explain how to use the unscented transformation in the filtering algorithm. The idea is to replace the formulae [START_REF] Fernández | A projection semiimplicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF] used in the original Kalman algorithm by the evaluations of the empirical covariance and mean value of the particles. More precisely, P XZ is replaced by Cov α (Y ( * ) , Z ( * ) ), with

Y (i) = A(X (i) ) and Z (i) = H(Y (i) ), and P Z is replaced by W + Cov α (Z ( * ) , Z ( * ) ).
To summarize, here is the UKF algorithm. Assume that the sigma-points I (i) , 1 ≤ i ≤ r, are given, as well as X+ 0 and P + 0 . For n ≥ 0, we have a 3 steps recursive algorithm:

• Sampling. Generation of the particles:

C n = P + n (Cholesky factorization) X(i)+ n = X+ n + C T n I (i) , 1 ≤ i ≤ r (26)
RR n°7657

• Prediction. Resolution of one time step of the model for each particle (can naturally be done in parallel) and computation of the empirical mean value and covariance:

     X(i)- n+1 = A n+1 ( X(i)+ n ) , 1 ≤ i ≤ r X- n+1 = E α ( X( * )- n+1 ) P - n+1 = Cov α ( X( * )- n+1 , X( * )- n+1 ) (27)
• Correction.

                   Z (i) n+1 = H( X(i)- n+1 ) , 1 ≤ i ≤ r P Z n+1 = Cov α (Z ( * ) n+1 , Z ( * ) n+1 ) + W n+1 P XZ n+1 = Cov α ( X( * )- n+1 , Ẑ( * )- n+1 ) K n+1 = P XZ n+1 (P Z n+1 ) -1 X+ n+1 = X- n+1 + K n+1 (Z n+1 -E α (Z ( * ) n+1
))

P + n+1 = P - n+1 -P XZ n+1 P Z n+1 -1 P XZ n+1 T (28)

Factorized formulation of the UKF

In this section, we show how the construction of P + n+1 , P XZ n+1 and P Z n+1 is modified when the covariance matrix P - n+1 can be factorized in the form LU -1 L T , with L ∈ M N,p and U ∈ M p,p . This will be useful to derive the reduced-order UKF algorithm in the next section.

We consider a family of sigma-points

I (i) ∈ R p , i = 1 . . . r. It is proved in [23, Sect. 2.2.1]
that the empirical covariance of the particles X (i) and Z (i) (see [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspective in medical imaging[END_REF]) can be expressed as

Cov α (X ( * ) , Z ( * ) ) = [X ( * ) ]D α [I ( * ) ] T [I ( * ) ]D α [I ( * ) ] T -1 [I ( * ) ]D α [Z ( * ) ] T (29)
where [X ( * ) ] ∈ M N,r is the matrix whose column i corresponds to the particle

X (i) , 1 ≤ i ≤ r, and where D α def = diag [α i ] ∈ M r,r
is the diagonal matrix containing the weights of formulae [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspective in medical imaging[END_REF].

The covariance matrix P - n+1 is by definition Cov α ( X( * )-n+1 , X( * )-n+1 ). Thus, according to (29), we rewrite it as

P - n+1 = L X n+1 P -1 α (L X n+1 ) T , (30)
with

L X n+1 def = [ X( * )- n+1 ]D α [I ( * ) ] T and P α def = [I ( * ) ]D α [I ( * ) ] T . (31)
Similarly, applying again (29) to P Z n+1 and P XZ n+1 (see their definition in ( 28)):

P Z n+1 = W n+1 + L Z n+1 P -1 α (L Z n+1 ) T , P XZ n+1 = L X n+1 P -1 α (L Z n+1 ) T , (32)
with

L Z n+1 def = [Z ( * ) n+1 ]D α [I ( * ) ] T . (33)
Thus, the a posteriori covariance P + n+1 can be written as:

P + n+1 = P - n+1 -P XZ n+1 P Z n+1 -1 P XZ n+1 T = L X n+1 P -1 α -P -1 α L Z n+1 T W n+1 + L Z n+1 P -1 α L Z n+1 T -1 L Z n+1 P -1 α L X n+1 T = L X n+1 U -1 n+1 L X n+1 T , ( 34
) RR n°7657
where U n+1 is defined, applying the Woodbury inversion lemma, by

U n+1 def = P α + L Z n+1 T W -1 n+1 L Z n+1 . (35)
Next, applying again the Woodbury inversion lemma to P Z n+1 , we obtain:

P Z n+1 -1 = W -1 n+1 -W -1 n+1 L Z n+1 U -1 n+1 L Z n+1 T W -1 n+1 . (36)
Thus, we deduce the factorized form of the Kalman-like matrix K n+1 = P XZ n+1 (P Z n+1 ) -1 :

K n+1 = L X n+1 P -1 α L Z n+1 T W -1 n+1 -W -1 n+1 L Z n+1 U -1 n+1 L Z n+1 T W -1 n+1 = L X n+1 U -1 n+1 L Z n+1 T W -1 n+1 . ( 37
)
The key point here is to be able to perform operations with P + n+1 and K n+1 using only L X n+1 , L Z n+1 and P α , and not directly P - n+1 . This property proves to be extremely useful when the rank of P - n+1 is much smaller than the size of X, as will be shown in the next section.

Remark 1

According to (9), K = P + H T W -1 . Thus, using the expression (34) of P + , we obtain K = L X U -1 L X T H T W -1 . By comparison with (37), we see that L Z formally corresponds to HL X , the "observed state sensitivity". Matrix L X (L θ ) -1 corresponds to the sensitivity of the state X with respect to the parameters θ. The factorization formula [START_REF] Simon | Kalman filtering with state constraints: A survey of linear and nonlinear algorithms[END_REF] is instrumental in identifying these analogies. We refer to [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] for more details.

Reduced-order UKF for parameter estimation

State and parameter estimation In principle, the UKF algorithm can easily be generalized to those cases when the uncertainties also affect the parameters. To do so, we consider an augmented state X = (X, θ) whose size is

N = N + p. At t = 0, X (0) = X 0 + ζ X = (X 0 , θ 0 ) + (ζ X , ζ θ
). The extended dynamical system reads:

X n+1 def = (X n+1 , θ n+1 ) = A n+1 (X n ) , (38)
with the initial condition X 0 = (X 0 , θ 0 ). When the parameters θ are constant in time, A n+1 is simply defined by

A n+1 (X n ) = (A n+1 (X n , θ n ), θ n ).
Then, even with a small number of parameters, classical Kalman-like filter cannot be applied to this augmented state because of the initial state size. However, if we assume that the uncertainties only affect the part θ of vector X , we will show in the next paragraph how to use a reduced order filter formulation which limits the computation of the filter operator to a subspace of small dimension -here the parameter space -and therefore can be computed.

Reduced-order UKF We now have all the material to present the reducedorder UKF algorithm that is used in this paper to estimate some parameters of a fluid-structure interaction system. As explained above, any Kalman-like algorithm can be used to estimate the state and the parameter by simply considering an augmented dynamical system RR n°7657 in the variable X = (X, θ). Nevertheless, as already mentioned, this approach is intractable with large system as those considered here. To circumvent this difficulty, the choice made in this study is to assume that the uncertainty is limited to the p parameters θ. Thus, the covariance matrix has a rank p and we can take advantage from the factorized formulation of Section 3.5.

Suppose that we do not have any uncertainties in the initial condition of the state, i.e., ζ X = 0. Then the initial covariance can be factorized:

P + 0 = 0 0 0 Cov(ζ θ ) = L X 0 L θ 0 U -1 0 (L X 0 ) T (L θ 0 ) T ,
where

U 0 def = Cov(ζ θ ) -1 is a (small) p × p matrix, L X 0 def = 0, L θ 0 def = ✶. Hence,
and this is a key point, we only need p + 1 sigma-points in R p for the sampling step in the UKF algorithm:

       C 0 = U -1 0 (Cholesky factorization) X(i)+ 0 = X+ 0 + L X 0 C T 0 I (i) = X+ 0 , 1 ≤ i ≤ p + 1 θ(i)+ 0 = θ+ 0 + L θ 0 C T 0 I (i) , 1 ≤ i ≤ p + 1 .
Reduced order estimation algorithm Consider the simplex sigma-points

I (i) , 1 ≤ i ≤ p + 1, their weights collected in the matrix D α def = diag [α i ] ∈ M p+1,p+1
. Then, for given values of X0 , θ0 and Cov(ζ θ ), perform • Initialization: initialize the sensitivities as

L θ 0 = ✶, L X 0 = 0, U 0 = Cov(ζ θ ) -1 (39a)
• Sampling: generate the particles by

     C n = U -1 n (Cholesky factorization) X(i)+ n = X+ n + L X n C T n I (i) , 1 ≤ i ≤ p + 1 θ(i)+ n = θ+ n + L θ n C T n I (i) , 1 ≤ i ≤ p + 1 (39b)
• Prediction: propagate the particles with the dynamical system

     ( X(i)- n+1 , θ(i)- n+1 ) = A n+1 ( X(i)+ n , θ(i)+ n ) X- n+1 = E α ( X( * )- n+1 ) θ- n+1 = E α ( θ( * )- n+1 ) (39c)
• Correction: use the innovation

Γ (i) n+1 = Z n+1 -H( X(i)- n+1 ) (1 ≤ i ≤ p + 1)
to correct the predicted state and parameters

                   L X n+1 = [ X( * )- n+1 ]D α [I ( * ) ] T ∈ M N,p L θ n+1 = [ θ( * )- n+1 ]D α [I ( * ) ] T ∈ M p L Z n+1 = [Γ ( * ) n+1 ]D α [I ( * ) ] T U n+1 = P α + (L Z n+1 ) T W -1 n+1 L Z n+1 ∈ M p X+ n+1 = X- n+1 + L X n+1 U -1 n+1 L Z n+1 T W -1 n+1 E α (Γ ( * ) n+1 ) θ+ n+1 = θ- n+1 + L θ n+1 U -1 n+1 L Z n+1 T W -1 n+1 E α (Γ ( * ) n+1 ) (39d)
RR n°7657

Remark 2

We have implicitly used some specific properties of the simplex sigma-points in the derivation of the algorithm as fully described in [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF]. For more general sigma-points, the algorithm is more intricate. We refer to [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] and in particular its erratum [START_REF] Moireau | Erratum of article "reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF].

Of course in practice, the state X also contains uncertainties. Although this problem will not be addressed in this paper, let us just mention that a possible solution consists of using the reduced order UKF in combination with another filtering technique, like a Luenberger observer [START_REF] Luenberger | An introduction to observers[END_REF]. This approach consists of adding a feedback term in the time-continous dynamics

Ẋ = A( X, θ) + K L (Z -H( X)) ,
where K L has to be designed for each physics and each type of observations in order to make the estimation error X -X decrease rapidly. We refer the reader to [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] where this approach has been proposed and applied to solid mechanics. For a study about Luenberger observers in fluid-structure problems, we refer to [3].

4 Reduced Order filtering for FSI problems

Fluid-Structure algorithm

Many procedures are known to couple the fluid (1) and the solid (2) problems (see, e.g., [START_REF] Fernández | Coupling schemes for incompressible fluid-structure interaction: implicit, semi-implicit and explicit[END_REF][START_REF] Fernández | Fluid structure interaction problems in haemodynamics, chapter 9[END_REF] and the references therein). In this work, we choose a semi-implicit partitioned algorithm. The scheme is called "partitioned" because it keeps independent the fluid and structure solvers and it is called "semiimplicit" because only a part of the fluid problem needs subiterations with the structure. A Newmark scheme is used in the solid and a first order Chorin-Temam projection method in the fluid. The coupling scheme treats explicitly the ALE-advection-diffusion step and implicitly the projection step, following the approach proposed in [START_REF] Fernández | A projection algorithm for fluid-structure interaction problems with strong added-mass effect[END_REF][START_REF] Fernández | A projection semiimplicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF].

Considering the variables ( ũn , y n f , y n , w n , a n ) which correspond to the fullydiscrete fluid velocity field, fluid domain displacement, solid displacement, solid velocity and solid acceleration, respectively, the complete fluid-structure algorithm reads: Given the initial condition ( ũ0 , y 0 f , y 0 , w 0 , a 0 ), find ( ũn , y n f , y n , w n , a n ), for n > 0, as follows:

1. For n = 0 initialize the pressure as:

-∆p 0 = - ρ f τ ∇ • ũ0 in Ω 0 f , ∇p 0 • n = 0 on Γ in , p 0 = p(t 0 ) on Γ out , ∇p 0 • n = - ρ f τ w -1 2 -ũ0 • n on Σ 0 , , (40)
with w -1 2 = w 0τ a 0 /2 and Ω 0 f = Ω f . 2. ALE-Advection-diffusion step (explicit coupling with the solid):

RR n°7657

• Update fluid domain:

y n+1 f = Ext( y n | b Σ ), Ω n+1 f = (I b Ω f + y n+1 f )( Ω f ), w n+1 f = ∂ τ y n+1 f in Ω n+1 f . (41)
• Fluid viscous sub-step:

               ρ f ũn+1 -ũn τ b x + ρ f ( ũn -w n+1 f ) • ∇ ũn+1 -2µ∇ • ǫ( ũn+1 ) = -∇p n in Ω n+1 f , ũn+1 = u in (t n+1 ) on Γ in , 2µǫ( ũn+1 )n f = 0 on Γ out , ũn+1 = w n-1 2 on Σ n+1 . ( 42
)
3. Projection step (implicit coupling with the solid): Initialize w

n+ 1 2 0 = w n-1 2 , and iterate (k > 0) until | y n+1 k -y n+1 k-1 | 2 + | w n+ 1 2 k -w n+ 1 2 k-1 | 2 < ε , on Σ n+1 .
with ε a given tolerance.

• Fluid projection sub-step:

               -∆p n+1 k = - ρ f τ ∇ • ũn+1 in Ω n+1 f , ∇p n+1 k • n = 0 on Γ in , p n+1 k = p(t n+1 ) on Γ out , ∇p n+1 k • n = - ρ f τ w n+ 1 2 k-1 -ũn+1 • n on Σ n+1 . (43)
• Solid:

                     ρ s ∂ τ w n+1 k -∇ • σ s ( y n+ 1 2 k ) = 0 in Ω s , ∂ τ y n+1 k = w n+ 1 2 k in Ω s , y n+1 k = 0 on Γ d , σ s ( y n+ 1 2 k )n s = 0 on Γ n , σ s ( y n+ 1 2 k )n s = -J n+1 f σ f ( ũn+1 , p n+1 k )(F f n+1 ) -T n on Σ. ( 44
)
4. Set

y n+1 = y n+1 k , w n+1 = w n+1 k , a n+1 = w n+1 -w n τ , n = n + 1
and go to Step 1.

In general, a matrix free Newton method is used to solve the interface equation resulting from the coupling between (43) and (44). In this work, since only linear solids are considered, the interface problem is simply solved with a matrix free GMRES algorithm. The fluid and the solid subproblems are dicretized in space with the P 1 finite element. The advection in the fluid is stabilized with a streamline diffusion term.

Algorithmic aspects of FSI filtering

Implementation Beside its efficiency, a very appealing feature of the proposed approach is the simplicity of its implementation. In algorithm (39a)-(39d), steps (39a), (39b) and (39d) are totally independent of the problem and can be implemented in an external software. Since this part is generic, it is very simple to change the number and the kind of parameters that are used in the filter.

Step (39c) corresponds to one time step of the dynamical system with a given initial condition. It therefore only requires to be able to "restart" the solvers from any state.

Particle initialization Let us give some more details about the "restart" of the fluid-structure algorithm, since this is an important aspect of the estimation procedure.

In our case, the discrete state X n corresponds to the degrees of freedom of ( ũn , y n f , y n-1 2 , w n-1 2 , a n ). At each time step, and for each particle separately, the following steps are performed:

• In the solid:

-Recover y n and w n by

y n = y n-1 2 + τ 2 w n-1 2 , w n = w n-1 2 + τ 2 a n .
• In the fluid:

-Reconstruct Ω n f = (I b Ω f + y n f )( Ω f ) -Reconstruct the pressure p n by solving                -∆p n = - τ ρ f ∇ • ũn in Ω n f , ∇p n • n = 0 on Γ in , p n = p(t n ) on Γ out , ∇p n • n = - ρ f τ w n-1 2 -ũn • n on Σ n . , (45)
Then perform steps 2 and 3 of the semi-implicit algorithm to get the new state by

X n+1 = ( ũn+1 , y n+1 f , y n+ 1 2 , w n+ 1 2 , a n+1 ).
Total vs. incremental ALE In many ALE solvers, the domain is updated in an incremental way, from Ω n to Ω n+1 , i.e.,

w n+1 f = Ext( w n-1 2 | Σ n ), Ω n+1 f = Ω n f + τ w n f .
This approach is not suitable in the present framework since it would results in incorrect fluid domains in each particle. This is the reason why we update the domain from Ω f using the Lagrangian displacement y n f . For very large displacements, this may complicate the computation of the mesh deformation (using for example a nonlinear pseudo-elasticity problem). In the applications presented in this work, a simple linear harmonic extension proved to be sufficient.

Windkessel boundary conditions

In blood flow simulations, the pressure p on Γ out is given by the ordinary differential equation ( 4). Its discrete form reads (assuming P V = 0):

             π n+1 = R d τ R d C + τ π n + R d C R d C + τ Qn , p(t n+1 ) = R p Qn+1 + π n+1 , Qn = Γ out ũn • n , (46)
where R p , R d and C are given constants. As already mentioned, when this kind of models is used for the fluid boundary conditions, π n has to be included in the state, i.e., X n = ( ũn , y n f , π n , y n-1 2 , w n-1 2 , a n ). When the system is "restarted", equation (46) 2 is used to define the value p(t n ) in (45) 3 .

Parameter range constraints For physical reason, the parameters are usually restricted to a subset of R p . For instance, Young's modulus, densities and viscosities have to be positive. This constraint can simply be enforced by reparametrizing the physical parameters in such a way that the estimation can indeed be done in the whole space R p . In this work, the physical parameters (e.g. the Young's modulus) are written as 2 θi , and the estimation is performed on θ i . Note that this reparametrization modifies the statistical meaning of the covariance. For example, when the covariance on θ is equal to identity, the Young's modulus 2 θ has the same probability to be twice or half as the initial value.

Choice of the covariance W n Due to the analogy with variational methods (i.e., comparing for instance [START_REF] Luenberger | An introduction to observers[END_REF] and ( 6)), we use W -1 n = γM Γm , with M Γm the L 2 mass matrix of the observation region Γ m , or its diagonal-lumped version. Hence, we propose to choose

γ = β w -1 (47)
where β is a positive scalar value that has to be set and

w -1 = τ obs σ -2 T ref H ref , (48)
where

T ref is a fixed reference time, H ref ≈ |Γ m |
, and τ obs is the time sampling of the observations, and assuming that ζ Z ∼ N (0, σ 2 I). These definitions is motivated by keeping dimensionless the tuning parameters associated with the measures and the a priori knowledge of the parameters. Note that the measure term "tends" to the continuous space-time norm when spatial and temporal resolution get higher.

5 Numerical examples

Forward model description

Consider the idealized abdominal aortic aneurysm (AAA) shown in Figure 2. Its length is 22.95 cm, its minimal diameter 1.7 cm (tubular part), its maximal diameter 4.98 cm (aneurysm) and its wall thickness 0.2 cm. The solid volume is divided into 5 zones (along the axis), associated with different Young's moduli: 2 θi M P a, with θ 1 = 1.5, θ 2 = 0.5, θ 3 = -1, θ 4 = 0.2, θ 5 = 1. Its Poisson ratio is 0.46, its density 1.2 g/cm 3 and viscosity η s = 10 -3 s. The fluid viscosity is 0.035 P o and its density 1 g/cm 3 . Figure 2: Idealized AAA: solid subregions (blue, green, red, cyan and magenta, respectively), fluid velocity field, pressure at the outlet and displacements at the aneurysm.

The boundary conditions are set as follows. The inlet and the outlet parts of the solid are clamped. In the fluid, a parabolic velocity profile is enforced at the inlet, with a realistic inflow (peak velocity ≈ 96 cm/s). At the fluid outlet, the three-dimensional model is coupled to a three-element Windkessel model (4) with R p = 700, R d = 5 • 10 3 and C = 2 • 10 -4 (all three in CGS unit system) and P V = 0. The time step is τ = 0.002 s.

At t = 0, the pressure is constant and equal to 80 mmHg, whereas all the other state variables are zero. During the whole simulation, the stress received by the structure is corrected by the initial one. Doing so, the solid only "feels" the difference with the diastolic phase. This is a simple way to account for the prestress in linear elastodynamics. We refer to [START_REF] Moireau | External tissue support and fluid-structure simulation in blood flows[END_REF], and references therein, for a discussion about prestress computation in a more general framework.

Estimation of the Young's modulus

Estimation setting We apply the reduced order UKF to estimate the Young's moduli, E i = 2 θi M P a, starting from θi = 0, i = 1, . . . , 5. The measurements are the displacements on all nodes of the fluid-structure interface Σ. In other words, the observation operator H consists of a m × n y matrix (whose entries are only 0 or 1), where n y is the number of displacement degrees of freedom of the whole solid and m is the number of displacement degrees of freedom at the interface Σ.

The initial parameter covariance is assumed to be

Cov(ζ θ ) = α✶ ( 49
)
where α is a given positive parameter. The noised signal is given by

Z n+1 = H Y d n+1/2 + σ ζ ,
where Y d n+1/2 is a synthetic solid displacement, ζ ∼ N (0, I), and σ = 0.2 mm, i.e., 10% of the maximal displacement at the solid wall.The measurements covariance is computed as indicated in Section 4.2. The scalar gain w -1 is obtained by formulae (48) setting T ref = 0.8 s, H ref = 76 cm 2 ≈ | Σ|. Moreover, we resample the noised measurements Z n+1 in time with τ obs = 10τ = 0.02s, and we re-interpolate linearly so that the filter can be applied at every time step of the simulation. The perfect signals and the noise (without time resampling for the sake of clarity) are shown in Figure 3 for representative points of the five regions of the vessel. Note that the signal-to-noise ratio decreases when going from the center to the ends of the AAA.

First estimation results

In Figure 4, we present two estimation results for E i = 2 θi obtained for two different covariances α = 4 and α = 9 for θ i in (49). We remind that the higher α the lower the confidence in the initial values of the parameters. The constant β in (47) is for the moment fixed to 1, which corresponds to a rather low confidence in the measurements. Several comments are in order. First, comparing the two graphs of Figure 4, we observe that the estimation algorithm has more freedom to adapt the parameters for large values of α. Second, the result is very good in region 3 but seems poor in the other regions. This is of course a consequence of the low SNR and the low value of β. Nevertheless, even in those cases, the algorithm is able to distinguished regions with high and low stiffness. This can be considered as a promising result for applications like detection of arteriosclerosis plaques from very noisy measurements. Sensitivity with respect to β The results presented in Figure 4 can be improved by increasing the value of β, which means increasing our confidence in the measurements. Figure 5 shows the behavior of the estimated parameter at the end of the cardiac cycle when β varies.

As expected, the sensitivity to β is higher for large a priori covariances Cov(ζ θ ) = α✶, and the estimated value is much more sensitive to α in region 1, which has the poorest SNR. But we also observe that for β reasonably large (about 10), the values of E are correctly estimated in all the other regions. Hence, at least for this example, we can conclude that, when the noise is not too RR n°7657 high, the estimation does not strongly depend on the user-defined parameters α and β. Estimation of the covariance The results presented so far corresponds to the estimated mean value of the parameters. An important feature of Kalman filtering is to also provide the estimated covariances, which is a valuable information about the confidence we can have in the results.

In the reduced UKF framework, the a posteriori covariance of the parameters is given by

P θ n+1 = L θ n+1 U -1 n+1 L θ n+1 T .
Figure 6 shows the results for the parameters mean θn and mean plus/minus the standard deviation θn ± diag(P θ n ). For the sake of clarity, the results are presented here in terms of θ instead of 2 θ , and only for the worst and best regions (1 and 3 respectively). We observe that the a posteriori variance is much higher in the region where the parameter is poorly estimated. We are therefore informed by the algorithm that the results are less reliable in that case, even when the correct value is unknown. This point is particularly important when dealing with real data. 5). The thick continuos lines represent the mean value, the thin continous lines the mean plus/minus standard deviation, and the dashed lines the reference value as before.

Estimation of the Windkessel's proximal resistance

As indicated in section 4.2, this estimation strategy does not involve many changes in the original software and can be implemented in a generic way outside the specific solvers. It is therefore quite simple to change or add parameters to estimate. In order to illustrate this versatility, we estimated the Windkessel's proximal resistance R p from the same measurements as before, i.e. synthetic noisy wall displacements. For this preliminary test, all the other parameters were supposed to be known. More complex situations will be considered in future works. As for the Young's modulus, we reparametrize R p as R p = 500•2 θ , and we perform the estimation on θ. The result shown in Figure 7 was obtained with Cov(ζ θ ) = α✶, with α = 1. The value of β = 100 was chosen so that the final estimation result is not sensible to it anymore. Note that can perfectly recover the reference value of R p .

Error in the initial condition and filtering Windkessel's pressure

Up to now, it has been assumed that the initial state is perfectly known. This is of course not the case in practice, and we would like to illustrate in this section the effect of an error in the initial condition. We will also show that the results can be improved by considering the Windkessel's pressure as a part of the filtered variables in the reduced UKF algorithm (we remind that up to now the filtered variables were restricted to the parameters).

Estimation with inexact initial condition

To generate the inexact initial condition, the forward model is first run for one cardiac cycle with an homogeneous Young's modulus corresponding to θ i = 0. The initial condition X0 used for the estimation corresponds to the state obtained at the end of the cycle (t = 0.8 s). Compared to the exact solution, the pressure difference is 2.4 mmHg and the displacement difference is 0.16 mm in the AAA center (see Figure 8). The velocity and displacement differences in the fluid and the solid are negligible.

Figure 9 (left) shows the parameter estimation results obtained with α = 4, β = 100 (the value of β is chosen so that we are in the flat zone shown in Figure 5) . Comparing with Figure 4, we see how the perturbation of the initial state deteriorates the estimation results. The pressure is lower due to the error in the initial condition (Figure 10, right, dashed cyan curve), which affects the parameter estimation.

As already mentioned, a first solution to address this problem is to filter the state with a physical Luenberger feedback, like in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] for elastodynamics and [3] for FSI. The complexity of Luenberger filters remains moderate when the system dimension increases since it is based on regular mechanical stabilization terms.

In the next section we propose a simpler solution, less general, but apparently quite efficient for blood flow simulations.

Filtering the Windkessel's pressure So far, the reduced order UKF algorithm has only been applied to the parameters. Now we propose to also apply it to the Windkessel's pressure π which is a state variable of the problem. Contrary to the other state variables -the fluid and solid velocities and displacements, which depend on the spatial discretizationπ n has by definition a limited size: just one real value in the experiment considered here, or a few real values in a case with several outlets. Therefore, it is straightforward to include it in the parameter space. From a practical viewpoint, this only consists of considering π n as a (time dependent) parameter in the above algorithm. If a new processor can be affected to this new "parameter", the computational cost remains almost unaffected. From the physical point of view, the fluid pressure is mainly governed by the Windkessel's pressure. Hence, improving the knowledge of this quantity would have a global impact on the estimation of the whole FSI system. Figure 9 (right) shows the estimation results in the case with an initial condition error (see the previous paragraph). We observe that the result is significantly improved.

Figure 10 shows the displacements and the pressure for different cases: reference (the solution to recover), not filtered, estimated with a perfectly known initial condition, estimated with error in the initial condition with and without filtering the Windkessel's pressure. Note that in all the filtered cases, the displacements are properly recovered. But we observe a substantial improvement in the pressure estimation when the Windkessel's pressure is filtered. This good estimation of the general amount of stress in the system, combined with the displacement accuracy, implies that the Young's modulus estimation is better than in the situation when π n is not filtered.

Conclusions

We have derived a procedure to estimate uncertain physical parameters in systems involving the mechanical interaction of a viscous incompressible fluid and an elastic structure. The method is based on the reduced-order Unscented Kalman Filter, introduced in [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF]. The algorithm does not need any tangent or adjoint problems and can easily be run in parallel, using as many processors as the number of parameters plus one. Doing so, the computational time needed to solve the inverse problem is of the same order as the time needed by a forward simulation using one processor. Some numerical results have been presented for noised synthetic data corresponding to an idealized geometry of aneurysm. We have estimated the artery Young's modulus in 5 regions using measurements of the wall displacement. We have first considered a problem with a perfectly known initial condition, and we have investigated the sensitivity of the results to the noise and the a priori parameter covariance matrices. Then, we have perturbed the initial condition of the system in order to illustrate its impact on the parameter estimation. Finally, we have shown that considering the Windkessel's pressure as a parameter is an easy way to improve the estimation performance, without significantly increasing the computational cost, and without any additional measurements.

To illustrate the versatility of the method, we have also presented a preliminary result about Windkessel resistance estimation. Future works will include the estimation of other parameters of the problem. Other kinds of measurements will also be considered, in particular in the fluid, like those obtained by catheters or phase-contrast MRI, which should improve the global observability -or identifiability -of the system.

A Simplex sigma-points

In this work we used the so-called simplex sigma-points I (i) , r = p + 1 , that correspond to vectors of zero mean and √ p covariance in R p . Together with the weights α i , they are computed recursively as (see [START_REF] Pham | Filtres de Kalman singuliers évolutifs pour l'assimilation de données en océanographie[END_REF][START_REF] Hoteit | A simplified reduced order Kalman filtering and application to altimetric data assimilation in Tropical Pacific[END_REF])

[I * 1 ] = - 1 √ 2α , 1 √ 2α , α = 1 p + 1 = α i ∀i ,

Figure 1 :

 1 Figure 1: Reference (left) and current domain (right) in the ALE formulation.

Figure 3 :

 3 Figure 3: Noise compared to the typical wall displacements in the five regions. Time responses (left), and signal-to-noise ratio (SNR: signal mean divided by standard deviation) (right).

Figure 4 :

 4 Figure 4: Results for the estimation of the Young's modulus with the reduced UKF algorithm for α = 4 (left) and α = 9 (right) for β = 1. The dashed lines correspond to the correct values. The color are the same as in Figure 3

Figure 5 :

 5 Figure 5: Results for the estimation of the Young's modulus (at t = 0.8) with the reduced UKF algorithm for α = 4, 9, for each section and different values of β. The relation color-region is the same as in Figure 3. The dashed lines represents the reference values.

Figure 6 :

 6 Figure6: Results for θn for Regions 1 (blue) and 3 (red) with α = 4, β = 100 (left) and α = 9, β = 20 (right) (β chosen in both cases to be in the flat region of figure5). The thick continuos lines represent the mean value, the thin continous lines the mean plus/minus standard deviation, and the dashed lines the reference value as before.

Figure 7 :

 7 Figure 7: Results for the estimation of the proximal Windkessel resistance Rp with β = 100.

Figure 8 :

 8 Figure8: Response of the forward model with θ i = 0, ∀ i, for displacements at the AAA center (left) and proximal pressure p at Γ out (right) . The dashed lines show the gap between t = 0 s and t = 0.8 s.

Figure 9 :Figure 10 :

 910 Figure 9: Parameter estimation in presence of an error in the initial condition with α = 4, β = 100. Only 5 stiffness parameters included in the parameter space for the reduced UKF (left), and and effect of inclusion the Windkessel's pressure (right).

Acknowledgement-

  The authors wish to thank Dominique Chapelle and Miguel Fernández for fruitful discussions. The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement nÂř 224495 (euHeart project).

, 2 ≤

 2 d ≤ p + 1 .

RR n°7657