Emmanuel Maitre

Thomas Milcent

Georges-Henri Cottet

Annie Raoult

Yves Usson

Applications of level set methods in computational biophysics

  L'archive ouverte pluridisciplinaire

Introduction

Biophysics and biomechanics are two fields where Fluid-Structure interactions play an important role, both from the modeling and computing points of view. In many 3D applications, flow and solid models coexist with biochemical systems. For such problems it is desirable to have at hand computing tools wich readily couple models of different nature (typically Eulerian for fluids, Lagrangian for solids), easily enforce continuity conditions at the interface and enable to handle reaction-diffusion systems.

In a series of papers [START_REF] Cottet | A level-set formulation of immersed boundary methods for fluidstructure interaction problems[END_REF][START_REF] Cottet | A level-set method for fluid-structure interactions with immersed surfaces[END_REF][START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF] such models were derived to compute the interaction of 3D incompressible fluids with elastic membranes or bodies. These models are based on Eulerian formulations of elasticity and rely on the use of level set functions, both to capture the fluid-solid interfaces and to measure the elastic stresses. Interface conditions are implicitly enforced through the elastic forces acting on the flow equations. They can be seen as alternative to more conventional ALE methods where Eulerian and Lagrangian formulations of the fluid and the solid are coupled through explicit enforcement of interface conditions.

The goal of the present paper is to present applications of this method to two types of problems arising in biophysics. The first problem is the computation of equilibrium shapes of biological vesicles. In this case the fluid structure model is a dynamical model for shape optimization, in the spirit of [START_REF] Biben | An advected field approach to the dynamics of fluid interfaces[END_REF], in contrast with more classical geometric approaches [START_REF] Simon | Existence of surfaces minimizing the Willmore functional[END_REF][START_REF] Du | A phase field formulation of the Willmore problem[END_REF]. Elastic stresses and immersion of the vesicle in an incompressible fluid are used to enforce constraints of constant area and volume.

In the second problem we are concerned with numerical simulations of spontaneous cardiomoyocyte contractions. In that case, the model couples an incompressible anisotropic medium with a reactiondiffusion system for the calcium concentrations. This coupling is through a calcium dependent active stress in the elastic medium. Our approach differs from that in [START_REF] Pustoc'h | Finite element modelling of the calcium-induced contraction cardiomyocytes based on time-lapse videomicroscopy[END_REF][START_REF] Okada | Three-dimensional simulation of calcium waves and contraction in cardiomyocytes using the finite element method[END_REF] by the fact that no remeshing of the structure is needed during the cell deformation. This leads to significant computational savings.

An outline of the paper is as follows. In Section 2 we focus on the problem of equilibrium shapes for biological vesicles. We present our level set formulation for the shape optimization and show numerical results for biological vesicles. In Section 3 we turn to the elastic deformation of a cardiomyocyte. We recall the level set Eulerian formulation derived in [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF] for a transverse isotropic elastic body. We couple this model with the reaction diffusion model [START_REF] Goldbeter | Minimal model for signal-induced Ca 2+ oscillations and for their frequency encoding through protein phosphorylation[END_REF] and show numerical results illustrating the association of cell contraction with calcium waves. Section 4 is devoted to some concluding remarks.

Equilibrium shapes of 3D phospholipidic vesicles

Phospholipidic vesicles are routinely considered as physical models in particular for red blood cells. Their membrane is a bilayer made of a fixed amount of molecules. As a result, it only responds to change of area or breakup. Taking into account the hydrodynamics is necessary in order to be able to study the behavior of these 3D cells when they are immersed in a flowing fluid. Phase-field models have been developed and used in 2D simulations, but 3D simulation of vesicles dynamics in shear flow is still a challenging problem. As a first step in this direction we consider here the problem of finding equilibrium shapes of these vesicles, constrained to have a fixed volume and area. An important parameter is the volume ratio

η = 3V (4π) 1/2 A 3/2 (1)
where V is the volume and A is the area, which measures the ratio between the volume of the cell and the volume of the sphere having the same area. Our computational approach mimics the underlying biophysical dynamics, in the sense that we assume that the cell is moving in an incompressible fluid and is subject to a very stiff elastic stress localized on the membrane. The curvature energy that the vesicle is supposed to minimize is used to derive an external force driving it towards its equilibrium.

Level-set formulation

In this section we give a level set formulation for the curvature driven dynamics of elastic membranes immersed in an incompressible fluid. Consider a domain Ω of R 3 containing some incompressible fluid into which a vesicle is immersed. This vesicle is considered to be an elastic surface. One way to describe the motion of this surface is to introduce a function ϕ whose zero level-set is the surface [START_REF] Osher | Level set methods and Dynamic Implicit Surfaces[END_REF]. Given a signed distance function ϕ 0 such that the initial interface is given by Γ 0 = x ∈ Ω, ϕ 0 (x) = 0 , the problem of localizing the structure is reduced to an advection of the function ϕ by the fluid velocity u. The velocity is solution to a Navier-Stokes system with a singular source term which accounts for the elastic forces acting on the fluid. As observed in [START_REF] Cottet | A level-set formulation of immersed boundary methods for fluidstructure interaction problems[END_REF], when u is incompressible, the change of area of {ϕ = 0} is recorded in |∇ϕ|. Following [START_REF] Cottet | A level-set method for fluid-structure interactions with immersed surfaces[END_REF], this makes it possible to express the area energy in terms of ϕ alone:

E a [ϕ] = Ω E(|∇ϕ|) 1 ε ζ( ϕ ε )dx
and the associated force is given by

f a [ϕ] = P ∇ϕ ⊥ ∇[E (|∇ϕ|)] -E (|∇ϕ|)κ(ϕ) ∇ϕ |∇ϕ| |∇ϕ| 1 ε ζ( ϕ ε ). (2)
In this formula ζ is a cut-off function classically used in level-set methods to spread the singular force on the mesh, P ∇ϕ ⊥ is the orthogonal projector on ∇ϕ ⊥ and κ(ϕ) = div ∇ϕ |∇ϕ| is the mean curvature. The constitutive law r → E (r) describes the response of the membrane to a change of area. As we already pointed out, in the present application the membrane is nearly inextensible, and we choose E (r) = λ(r -1) for large values of λ. This somehow corresponds to penalizing the change of area.

In order to account for curvature effects, we next introduce the following energy:

E c [ϕ] = Ω G(κ(ϕ))|∇ϕ| 1 ε ζ( ϕ ε )dx
where |∇ϕ| 1 ε ζ( ϕ ε )dx is an approximation of the surface measure. A common choice is G(r) = 1 2 r 2 but terms identification is made easier in the following by keeping a general function G. The strategy to compute the curvature force is to take the time derivative of the energy and identify it with the power of the generated force:

dE c dt = - Ω f c • udx.
Let us first compute the differential of E c :

dE c [ϕ](δ) = Ω G (κ(ϕ)) div ∇δ |∇ϕ| - ∇ϕ • ∇δ |∇ϕ| 3 |∇ϕ| 1 ε ζ( ϕ ε )dx + Ω G(κ(ϕ)) ∇ϕ • ∇δ |∇ϕ| 1 ε ζ( ϕ ε ) + G(κ(ϕ))|∇ϕ| 1 ε 2 ζ ( ϕ ε )δdx.
The two terms in the second integral may be combined: upon integrating the first one by parts, one obtains

- Ω G(κ(ϕ))κ(ϕ) 1 ε ζ( ϕ ε )δ + G(κ(ϕ)) ∇ϕ |∇ϕ| 1 ε 2 ζ ( ϕ ε )∇ϕδ + ∇G(κ(ϕ)) • ∇ϕ |∇ϕ| 1 ε ζ( ϕ ε )δ.
Therefore

dE c [ϕ](δ) = Ω G (κ(ϕ)) div P ∇ϕ ⊥ (∇δ) |∇ϕ| |∇ϕ| 1 ε ζ( ϕ ε )-G(κ(ϕ))κ(ϕ) 1 ε ζ( ϕ ε )δ-∇G(κ(ϕ))• ∇ϕ |∇ϕ| 1 ε ζ( ϕ ε )δdx
which from the expression of the mean curvature κ(ϕ) also reads

dE c [ϕ](δ) = Ω G (κ(ϕ)) div P ∇ϕ ⊥ (∇δ) |∇ϕ| |∇ϕ| 1 ε ζ( ϕ ε ) -div(G(κ(ϕ)) ∇ϕ |∇ϕ| ) 1 ε ζ( ϕ ε )δdx.
Since P ∇ϕ ⊥ (∇δ) • ∇ϕ = 0 the first term may be integrated by parts to give

- Ω ∇ |∇ϕ|G (κ(ϕ)) • P ∇ϕ ⊥ (∇δ) 1 |∇ϕ| 1 ε ζ( ϕ ε ) = - Ω P ∇ϕ ⊥ |∇ϕ|∇G (κ(ϕ)) • ∇δ |∇ϕ| 1 ε ζ( ϕ ε )dx
where the symmetry of the projector on ∇ϕ ⊥ has been used. Integrating once more by parts we find

dE c [ϕ](δ) = Ω div -G(κ(ϕ)) ∇ϕ |∇ϕ| + 1 |∇ϕ| P ∇ϕ ⊥ ∇[|∇ϕ|G (κ(ϕ))] 1 ε ζ( ϕ ε )δdx.
Let us now compute the time derivative of the energy. Using the advection equation on ϕ,

d dt E c [ϕ] = dE c [ϕ](ϕ t ) = dE c [ϕ](-u • ∇ϕ) = - Ω f c (x, t) • udx (3)
which by identification gives:

f c [ϕ] = div -G(κ(ϕ)) ∇ϕ |∇ϕ| + 1 |∇ϕ| P ∇ϕ ⊥ ∇[|∇ϕ|G (κ(ϕ))] 1 ε ζ( ϕ ε )∇ϕ. (4)
A more general derivation of curvature driven level set models for shape optimization will be given in [START_REF] Maitre | A comparison between level-set and classical shape optimization for curvature dependent functionals[END_REF].

The two forces ( 2) and ( 4) are finally inserted as forcing terms in the Navier-Stokes equations, leading to the following model: given and initial velocity field u 0 and an initial interface ϕ 0 , find (u, ϕ) solution to

           ρ(ϕ)(u t + u • ∇u) -div(µ(ϕ)D(u)) + ∇p = f a [ϕ] + f c [ϕ] on Ω×]0, T [ div u = 0 on Ω×]0, T [ ϕ t + u • ∇ϕ = 0 on Ω×]0, T [ u = u 0 ϕ = ϕ 0 on Ω × {0}
where, in the case of homogeneous interior and exterior fluids we have

µ(ϕ) = µ 1 H( ϕ ε ) + µ 2 (1 -H( ϕ ε )) and ρ(ϕ) = ρ 1 H( ϕ ε ) + ρ 2 (1 -H( ϕ ε )) + ν 1 ε ζ( ϕ ε )
where ν is the membrane surface density (in its reference state) and

H(r) = r -∞ ζ(s)ds.
The boundary condition to be enforced on ∂Ω plays a marginal role. For simplicity we in general choose homogeneous Dirichlet boundary conditions.

Note that for this model the following energy equality holds for all t ∈ [0, T ]:

1 2 Ω ρ(ϕ)u 2 dx + Ω E(|∇ϕ|) 1 ε ζ( ϕ ε ) dx + Ω G(κ(ϕ))|∇ϕ| 1 ε ζ( ϕ ε )dx + 1 2 t 0 Ω µ(ϕ)D(u) 2 dx dt = 1 2 Ω ρ(ϕ 0 )u 2 0 dx + Ω E(|∇ϕ 0 |) 1 ε ζ( ϕ 0 ε ) dx + Ω G(κ(ϕ 0 ))|∇ϕ 0 | 1 ε ζ( ϕ 0 ε )dx
which shows that the spreading of elastic and curvature forces inherent to the level-set method does not introduce any energy dissipation. We also remark that the resolution of the full problem fluid/membrane, while not mandatory to obtain equilibrium shapes (but mandatory to study the dynamical behavior of vesicles in shear flow), brings some advantages from the viewpoint of volume conservation: indeed we will use a projection method which will ensure this conservation at the discrete level. Note that in order to solve the minimization problem without any fluid, it is necessary to add a volume constraint which is usually enforced through a Lagrange multiplier approach. This may result in a loss of accuracy for volume conservation.

Numerical results

The numerical results presented here show two typical situations of optimal shapes for 3D vesicles. The final shape depends on the volume ratio coefficient r defined in [START_REF] Brown | Accurate Projection Methods for the Incompressible Navier-Stokes Equations[END_REF]. The first test case corresponds to r = 0.8, giving a peanut minimizing shape (figure 1, top pictures). The second one to r = 0.586 looks like a real red blood cell (figure 1, bottom pictures). In each series of pictures we have represented a sequence of shapes from the initialization stage (left pictures) to the steady-state optimal shapes (right pictures). These shapes qualitatively agree with those observed for corresponding η values ( [START_REF] Seifert | Shape transformations of vesicles: phase diagram for spontaneous curvature and bilayer coupling model[END_REF]).

Ongoing works deal with the extension of the present method to simulate the dynamical behavior of 3D vesicles in shear flows.

The numerical resolution of Navier-Stokes equations is performed using a finite differences solver (projection method) on a MAC mesh [START_REF] Brown | Accurate Projection Methods for the Incompressible Navier-Stokes Equations[END_REF] of size 128 3 . In order to ensure volume conservation which is crucial in this problem, the level-set function is advected using a fifth order WENO scheme. Since the level-set function is used through its gradient to compute the stretching, we do not perform the usual redistancing operation on the level-set function [START_REF] Osher | Level set methods and Dynamic Implicit Surfaces[END_REF]. Instead, we use the renormalization ϕ |∇ϕ| to measure the distance to interface. This approach was proved in [START_REF] Cottet | A level-set method for fluid-structure interactions with immersed surfaces[END_REF] to be efficient from the point of view of both volume conservation and interface force calculations.

Eulerian three-dimensional modelisation and simulation of cardiomyocyte contraction induced by calcium waves

In a recent article, Okada et al. [START_REF] Okada | Three-dimensional simulation of calcium waves and contraction in cardiomyocytes using the finite element method[END_REF] investigated the mechanism of calcium wave propagation in connection with cardiomyocyte contraction. They developped a 3D simulator using the model of [START_REF] Subramanian | Underlying mechanisms of symetric calcium wave propagation in rat ventricular myocytes[END_REF] for the Ca 2+ dynamics and relying on the Negroni and Lascano's contraction model [START_REF] Negroni | A cardiac muscle model relating sarcomere dynamics to calcium kinetics[END_REF] which couples Ca 2+ concentration with force generation. For the elastic part an isotropic Saint Venant-Kirchhoff hyperelastic model was assumed and myofibrills, Z-lines, sarcolemma, cytoskeleton and cytoplasm were represented by various finite element families. In our paper we adopt a similar approach in an Eulerian framework: following [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF] we use a level set approach of the fluid-structure coupling that occurs between the surrounding fluid and the cardiomyocyte, considering these two as a unique incompressible continuous medium. The microscopic internal structure of the cardiomyocyte is not described: the passive property of the myocyte is given by nonlinear elasticity, with a transverse isotropy assumption accounting for the topology of the sarcolemma. The calcium dynamics is coupled through an active stress law given by Stuyvers et al. [START_REF] Stuyvers | Effect of stimulation rate, sarcomere length and Ca 2+ on force generation by mouse cardiac muscle[END_REF], as described in Tracqui et al. [START_REF] Pustoc'h | Finite element modelling of the calcium-induced contraction cardiomyocytes based on time-lapse videomicroscopy[END_REF]. While our model does not pretend to reproduce the internal structure as precisely as in [START_REF] Okada | Three-dimensional simulation of calcium waves and contraction in cardiomyocytes using the finite element method[END_REF], it is more realistic in some respects in the elastic part. In particular it is worth noticing that the Saint Venant-Kirchhoff constitutive law considered in [START_REF] Okada | Three-dimensional simulation of calcium waves and contraction in cardiomyocytes using the finite element method[END_REF] should not be used for the large deformations observed in myocytes.

Description of the model

In this section we rely on the level set framework developed in [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF] for anisotropic elastic bodies in interaction with incompressible fluids and couple this model with a differential system for the Calcium concentration responsible for the active stress.

Eulerian formulation of the passive continuous medium

The cardiomyocyte is immersed in a fluid which lies in a bounded domain Ω ⊂ R 3 . We denote by u the divergence-free velocity field of the whole continuous medium, assumed to be C 1 and to vanish on the boundary ∂Ω.

(H) u ∈ C 1 (Ω × [0, T ]) and u = 0 on ∂Ω × [0, T ]
The interface between the cardiomyocyte and the fluid is captured by a level set function ψ 0 . Note that ψ 0 needs not be a signed distance function. In our calculations it was obtained from experimental data, specifically by confocal miscroscopy [START_REF] Usson | Method for the study of the threedimensional orientation of the nuclei of myocardial cells in fetal human heart by means of confocal scanning laser microscopy[END_REF]. Let us introduce the characteristics of the vector field u. We denote by s → X(s; x, t) ∈ R 3 the solution of the differential system ∂X ∂s = u(X, s)

with "initial" condition X(t; x, t) = x. Classically, under the assumptions (H), the map x → X(s; x, t) is a C 1 diffeomorphism from Ω to Ω. Since u is incompressible, one has div u = 0 and thus the Jacobian of X, denoted by J, is equal to 1. Following [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF], to compute X in an Eulerian fashion, we use the following transport equation satisfied by X as a funtion of t, x :

X t (s; x, t) + u(x, t) • ∇X(s; x, t) = 0 (5)
still with the same initial condition on t = s. Once X is computed at time t (for s = 0), several quantities can be easily obtained. The cardiomyocyte boundary position is given by the zero level set of ψ(x, t) = ψ 0 (X(0; x, t)). The left Cauchy-Green tensor is given by ([2], p. 15, [START_REF]Ciarlet Mathematical elasticity[END_REF], p. 43) B = F F T where F (x, t) = (∇X)(t; X(0; x, t), 0) = (∇X) -1 (0; x, t)

the last equality being obtained by differentiation of X(t; X(0, x, t), 0) = x. For an elastic material whose response is isotropic at point ξ = X(0; x, t), the Cauchy stress tensor at x is given by ([2], p. 50, [START_REF]Ciarlet Mathematical elasticity[END_REF], p. 115) T (x) = T D (X(0; x, t), B(x, t)).

If this material is incompressible, then this constitutive equation becomes ( [START_REF] Gurtin | Topics in finite elasticity[END_REF], p. 45, or [START_REF]Ciarlet Mathematical elasticity[END_REF], p. 259 after applying a Piola transform):

-p(x)I + T D (X(0; x, t), B(x, t))

where I is the identity. If we describe the surrounding fluid as Newtonian, the stress tensor in this part of the continuous medium is given by -p(x)I + µD(u) where µ stands for the viscosity and D(u) = 1 2 (∇u + ∇u T ). The conservation of momentum may thus be written as

ρ(u t + u • ∇u) -div σ + ∇p = f
where σ = σ S χ {ψ<0} + σ F χ {ψ>0} and σ S = T D (X(0; x, t), B(x, t)), σ F = µD(u). Of course the momentum equation has to be understood in the sense of distributions. For computational purposes, a regularization of χ using the level set function ψ must be introduced. If we denote as above by r → H( r ε ) an approximation of the Heaviside function, the regularized stress

σ ε = σ S (1 -H( ψ ε )) + σ F H( ψ ε )
varies smoothly across the interface ψ = 0 and the momentum equation may be understood in the classical sense.

Transverse isotropy

In the myocardic tissue, cells have a rod-like shape and are assembled along fibres and bundled by collagen. Thus anisotropy is expected. More precisely, following [START_REF] Mourad | Description topologique de l'architecture fibreuse et modlisation mcanique du myocarde[END_REF], p. 80, the myocardic tissue can be considered as transverse isotropic [START_REF] Ogden | Nonlinear elasticity, anisoptropy, material staility and residual stresses in soft tissue, Biomechanics of Soft Tissue in Cardiovascular Systems[END_REF]. Let τ be the preferred direction for the cardiomyocyte, i.e. its long axis at time t = 0. Such a material can be characterized with a strain energy which depends on F and τ ⊗ τ . The stress tensor has then the following general expression

σ S = -pI + 2α 1 B + 2α 2 (tr(B)B -B 2 ) + 2α 4 F τ ⊗ F τ + 2α 5 (F τ ⊗ BF τ + BF τ ⊗ F τ ) (6)
where α i is the derivative of the strain energy with respect to the invariant number i, with

I 1 = tr(B), I 2 = 1 2 [tr(B) 2 -tr(B 2 )], I 4 = |F τ | 2 , I 5 = (BF τ ) • (F τ ). (7)
While I 1 2

4 stands for the fibre elongation (preferred direction), I 1 2

5 carries also an information on the elongation in the direction normal to the unpreferred directions. For the sake of simplicity we will restrict ourselves in this paper to the case where α 5 = 0. As the quantity computed by resolution of ( 5) is (x, t) → X(0; x, t), we have to express the stress law in terms of its components X i (0; x, t), i = 1, 2, 3. As F = (∇X) -1 (0; x, t), and det F = 1 by incompressibility, this is easily done:

F (x, t) = cof ∇X T =   X ,x 2 × X ,x 3 X ,x 3 × X ,x 1 X ,x 1 × X ,x 2   = ∇X 2 × ∇X 3 ∇X 3 × ∇X 1 ∇X 1 × ∇X 2
where X ,x i × X ,x j are row vectors, and ∇X i × ∇X j column vectors. The components of B = F F T are thus obtained from two-by-two scalar products of the X ,x i × X ,x j . For the invariants I 1 and I 2 involved in [START_REF] Du | A phase field formulation of the Willmore problem[END_REF], after some elementary computations, there holds

I 1 = tr(B) = |X ,x 2 × X ,x 3 | 2 + |X ,x 3 × X ,x 1 | 2 + |X ,x 1 × X ,x 2 | 2 = | cof ∇X| 2 tr(B 2 ) = tr(B) 2 -2(|X ,x 1 | 2 + |X ,x 2 | 2 + |X ,x 3 | 2 ), I 2 = |X ,x 1 | 2 + |X ,x 2 | 2 + |X ,x 3 | 2 = |∇X| 2 (8)
After some tedious yet elementary algebra, using (X ,x 1 × X ,x 2 ) • X ,x 3 = 1, (tr B)B -B 2 has the following simple expression:

(tr B)B -B 2 =   |X ,x 2 | 2 + |X ,x 3 | 2 -X ,x 1 • X ,x 2 -X ,x 1 • X ,x 3 -X ,x 1 • X ,x 2 |X ,x 1 | 2 + |X ,x 3 | 2 -X ,x 2 • X ,x 3 -X ,x 1 • X ,x 3 -X ,x 2 • X ,x 3 |X ,x 1 | 2 + |X ,x 2 | 2   and α 1 , α 2 , α 4 are functions of | cof ∇X| 2 , |∇X| 2 , | cof ∇X T τ | 2 .

Active contraction and final model

We now come to the coupling of the elastic properties with the biochemistry taking place inside the cell. For the active behavior of the cardiomyocyte, we follow [START_REF] Pustoc'h | Finite element modelling of the calcium-induced contraction cardiomyocytes based on time-lapse videomicroscopy[END_REF] where the active stress is added to σ S . With our notations it corresponds to adding in (6) the term

T 0 γ(Z(x, t))
to α 4 , where Z is the intracellular Ca 2+ concentration, and γ is the following Hill function [START_REF] Stuyvers | Effect of stimulation rate, sarcomere length and Ca 2+ on force generation by mouse cardiac muscle[END_REF]:

γ(Z) = Z n H Z n H 50 + Z n H
The constant T 0 was fixed to 5.5kP a such that predicted and experimental amplitudes match. The calcium dynamics is given by the following reaction-diffusion system [START_REF] Subramanian | Underlying mechanisms of symetric calcium wave propagation in rat ventricular myocytes[END_REF][START_REF] Goldbeter | Minimal model for signal-induced Ca 2+ oscillations and for their frequency encoding through protein phosphorylation[END_REF]:

∂Y ∂t = ν 2 (Z) -ν 3 (Y, Z) -k f Y (9)
∂Z ∂t = ν 0 + ν 1 β -ν 2 (Z) + ν 3 (Y, Z) + k f Y -kZ + ∇ • (D∇Z) (10)
Z(r, t 0 ) = Z 0 ; Y (r, t 0 ) = Y 0 (11)
where the diffusion tensor is diagonal. As experimentally observed by [START_REF] Subramanian | Underlying mechanisms of symetric calcium wave propagation in rat ventricular myocytes[END_REF], we considered a ratio

D 33 D 11 = D 22
D 11 = 0.5 between diffusion along the sarcolemma direction and the transverse directions. The calcium fluxes ν 2 and ν 3 are given by the following Michaelis-Menten functions:

ν 2 = V M 2 Z n K n 2 + Z n (12)
ν 3 = V M 3 Y m K m R + Y m Z p K p A + Z p (13)
We refer to [START_REF] Pustoc'h | Finite element modelling of the calcium-induced contraction cardiomyocytes based on time-lapse videomicroscopy[END_REF] for precise biological meaning of constants and function appearing in ( 9)-( 13), and to table 2 for the values used in the forthcoming numerical simulations.

The final model consists of the following set of equations:

ρ(u t + u • ∇u) -div σ ε + ∇p = f, div u = 0 X t + u • ∇X = 0, F = cof ∇X T , B = F F T σ ε = σ S (1 -H( ψ ε )) + σ F H( ψ ε )
with σ S = 2α 1 B + 2α 2 (tr(B)B -B 2 ) + 2(α 4 + T 0 γ(Z(x, t)))F τ ⊗ F τ and σ F = µD(u), coupled with the system [START_REF] Gurtin | An introduction to continuum mechanics[END_REF][START_REF] Gurtin | Topics in finite elasticity[END_REF][START_REF] Maitre | A comparison between level-set and classical shape optimization for curvature dependent functionals[END_REF].

Parameter Value Unit ν 0 0.45 µM.s -1 k 2.2 s -1 ν 1 4 µM.s -1 V M 2 65 µM.s -1 V M 3 500 µM.s -1 K 2 1.2 µM K A 0.92 µM K R 3.5 µM Y 0 0.

Numerical results

For the Navier-Stokes and advection equations we used the finite-difference method described in Section 2.2. The reaction-diffusion system was solved on the same mesh using a seven-point discretisation of the (diagonal) diffusive terms. We considered two cases. In the first case the calcium concentration is homogeneous in the cell. In the second case the calcium concentration is initialized randomly and produces after some transient period a coherent wave propagating across the cardiomyocyte. In both cases the cardiomyocyte geometry was acquired from data obtained in [START_REF] Usson | Method for the study of the threedimensional orientation of the nuclei of myocardial cells in fetal human heart by means of confocal scanning laser microscopy[END_REF] using confocal microscopy (see figure 2). In all our calculation we used a grid of 128 3 points in a box surrounding the cell (the cell itself represented about 50% of the computational box). However, to obtain a better quality visualisation we ran a better resolved transport equation for the level set using velocity values interpolated from the low resolution results. Calcium concentration where represented on the fine mesh by interpolation from the lower resolution calculations.

Uniform contraction

In this first test case (figure 3.2.1), the coefficient β which controls the source of calcium is set constant in space. This results in an homogeneous calcium release and in a uniform contraction along the whole body of the cardiomyocyte. Note that this leads to a nonlinear elasticity problem with large displacements. The cardiomyocyte is discretized on a 384 3 grid for the advection equation, while the fluid-structure equations are solved on a 128 3 grid.

Coupling with a calcium wave

In this case, the coefficient β is a spatially localized function. This triggers a calcium wave which is starting on the left front of the cell and is propagating towards the other end (see figure 3.2.2). Note that the propagation is faster along the fiber direction, due to the higher diffusion coefficient in the FIGURE 3: Uniform contraction of a cardiomyocyte resulting from an homogeneous calcium release reaction-diffusion CICR system (see table). In that case the calcium peaks come with a deformation in a plane transverse to the principal axis of the cell. The efficiency, in terms of contraction along the principal axis, is clearly much lower than in the previous case. Note that for clarity, in this experiment as in the previous one, the fluid is not represented.

Conclusion

We have presented level set methods based on Eulerian representation of elasticity to deal with fluidstructure interactions. Two biophysical examples were provided to illustrate our method in this field, for vesicles shape optimization and cardiomyocyte contraction. These examples focus on interaction of a fluid with an immersed membrane, or an anisotropic elastic material. Both involve elasticity with large displacements, which would be very time consuming to deal with in the classical ALE method.

FIGURE 1 :

 1 FIGURE 1: Shape optimization for equilibrium shapes of biological shapes. Top pictures: η = 0.8, bottom pictures: η = 0.586. Initialization to steady-state from left to right. Computations made by Thomas Milcent.

FIGURE 2 :

 2 FIGURE 2: Geometry of a cardiomyocyte obtained by confocal imagery [25] as used to initialize the level set function ψ

FIGURE 4 :

 4 FIGURE 4: Calcium wave propagating into a cardiomyocyte. The color codes for the calcium concentration. Total simulation time: 1s

Table 2 :

 2 Parameters values used in numerical simulations

		1	µM
	Z 0	10	µM
	Z 50	2.5	µM
	β	0.05	-
	D 11	300	µm 2 s -1
	D 22	150	µm 2 s -1
	D 33	150	µm 2 s -1
	T 0	5.5	kP a

Acknowledgments

This work was partially supported by the French Ministry of Education though ACI program NIM (ACI MOCEMY contract # 04 5 290).


Mohamed Bouchakour

Guillaume Jeannic

Florent Autrusseau

JND MASK ADAPTATION FOR WAVELET DOMAIN WATERMARKING

One of the most challenging issues for watermarkers is to tune the strength of their watermark embedding. The strength is usually an α parameter which is increased until a reasonable trade-off between invisibility and robustness is achieved. The watermarking community needs efficient Just Noticeable Difference (JND) masks to optimally embed the watermarks. The Fourier transform is particularly adapted to the Human Visual System (HVS) modeling. In this work, we evaluate the usability of the JND mask in the wavelet domain. The use of the mask in the DWT domain involves some approximations. We will see here that the HVS decomposition and the wavelet decomposition do not perfectly fit altogether. The efficiency of the so obtained mask is tested both in terms of invisibility and robustness.

INTRODUCTION

It is widely admitted that among the different requirements needed in watermarking applications, the robustness and invisibility are very important. The watermark's robustness is inversely proportional to the invisibility. Thus, in watermarking context, optimising the invisibility versus robustness trade-off is crucial. Perceptual modeling is very important in watermarking context, it is in fact crucial to embed the watermark in the perceptually significant image areas, otherwise, the watermark wouldn't resist to lossy compression. Of course, as many studies have been conducted on the design of perceptual models for image compression techniques, several DCT or DWT based perceptual masks can be found in the watermarking literature [START_REF] Huang | Perceptual Watermarking Using a Wavelet Visible Difference Predictor[END_REF][START_REF] Chih-Wei | Exploring effective coefficients in transform-domain perceptual watermarking[END_REF]. However, most of the DCT perceptual masks are simply based on quantization matrices, and do not take into account more complex processes, such as masking effect. One of the most advanced DCT/DWT perceptual watermarking technique was proposed by Podilchuk and Zeng [START_REF] Podilchuk | Image-adaptive watermarking using visual models[END_REF]. The authors designed an image adaptive watermarking algorithm exploiting Watson's works for the perceptual masks implementation. Bartolini et al. have This work is partially supported by the ANR-Safimage project studied several perceptual masks in [START_REF] Bartolini | Mask building for perceptual hiding frequency embedded watermarks[END_REF]. The first one exploited a complex multiple channels HVS model, the second was based on the local variance computation, and the last one was made from heuristic considerations. They used a bank of filters to extract an appropriate frequency range, and a Sobel filtering detecting the image edges. The authors claimed that the masks based on heuristic consideration presented better detection results than HVS based masks. The watermark weighting has always been a very challenging problem. Authors in [START_REF] Voloshynovskiy | A Stochastic approach to content adaptive digital image watermarking[END_REF] introduce the Noise Visibility Function (NVF), a content adaptive watermarking embedding scheme was designed for noise-like watermarking embedding. However very few complex HVS models are used for watermarking purpose. Evidently, an important drawback of HVS models lie in their complexity. Providing to the watermarking community an efficient (HVS based) low complexity JND mask remains an open challenge.

The goal of this paper is to adapt a previously designed JND mask in the wavelet domain. Effectively, we recently proposed in [START_REF] Autrusseau | A robust image watermarking technique based on quantization noise visibility thresholds[END_REF] a JND mask based on quantization noise visibility thresholds, and an adaptation of the mask in a Fourier domain embedding technique. As we will see in section 2 there are some incompatibilities between the Fourier splitting of the human visual system model and the Fourier representation of the wavelet sub-bands. However, as explained later, the quantization thresholds being quite similar for neighboring visual sub-bands, in this work, we spread the watermark in one DWT sub-band, and verify that most of the Fourier representation of the watermark is included within at most two neighbouring visual bands. We hereby adapt the previous embedding technique operating in the Fourier domain into the wavelets. We evaluate the efficiency of the mask regarding both invisibility and robustness. The Stirmark benchmark is used to evaluate the robustness of the proposed algorithm and PSNR, wPSNR and C4 are used to assess the quality of the marked images.

This paper is structured as follows: Section 2 presents the HVS model that we use as well as the JND mask and its adaptation to Wavelets. In section 3 we present the adaptation into wavelets and the embedding technique. Finally, section 4 gives experimental results for both in-visibility and robustness.

HUMAN VISUAL SYSTEM MODEL

Based on psychophysics experiments conducted in our lab, we have derived a Perceptual Channel Decomposition (PCD). The filters of the PCD are similar to the cortex filters developed by Watson. Interested reader may refer to [START_REF] Autrusseau | A robust image watermarking technique based on quantization noise visibility thresholds[END_REF] for further details on the HVS model and JND masks. As pointed out in [START_REF] Zeng | An overview of the visual optimization tools in JPEG[END_REF], there are some incompatibilities between the HVS models decompositions and Wavelet sub-bands. Effectively, in the Fourier spectrum, sub-bands at 30 • and 150 • are treated separately along the visual pathways, and processed by different cells in the visual cortex, whereas, in the DWT domain this information is grouped into the same HH sub-band. As a consequence, for example, modeling the self-masking effect in these sub-bands is limited. One coefficient in the transform domain might represent different signals at 30 • and 150 • which do not mask each other. The PCD, defined in cycle/degree, in the frequency domain is given Figure 1-a. pixel location in the i th crown and the j th angular channel is defined as the ratio between the luminance of the reconstructed (i, j) sub-band at the specified pixel location and the mean luminance of all the below radial sub-bands. This gives the equation 1

C i,j (m, n)= L i,j (m, n) L i (m, n) , ( 1
)
where i represents the i th radial channel and L i (m, n) is the local mean luminance at the (m, n) position (i.e.

spatial representation of all Fourier frequencies below the considered visual sub-band). We use the local contrast definition to determine the allowable watermark strength. Previous studies, conducted on the perceptual decomposition 1-a, determined invisible quantization conditions. Each perceptual sub-band was independently quantized and the image quality was assessed by a set of observers. The optimal quantization step (∆C), which does not visually affect the image, has been introduced. The formula for computing the ∆C i,j values can be found in [START_REF] Autrusseau | A robust image watermarking technique based on quantization noise visibility thresholds[END_REF], and

∆L i,j (m, n)=∆C i,j × L i (m, n) . ( 2
)
∆L i,j (m, n) represent the JND mask computed for each (i, j) sub-band. Finally a perceptual weighting coefficient K i,j is computed from the watermark's spatial domain representation and the visual mask(eq. 3):

K i,j = argmin m,n ∆L i,j (m, n) W S (m, n) , ( 3
)
where W S (m, n) is the watermark's spatial representation before weighting process by factor K i,j for each (m, n) spatial position. It is very important to notice that the JND masks proposed by this technique are suitable for specific frequency contents i.e. for a chosen sub-band, the frequency content of the embedded watermark should ideally be totally restrained in the same frequency content. Actually, the frequency representation of the watermark should be fully included within a PCD sub-band, whereas the corresponding JND mask is entirely made of the lower frequency disk. Restraining the watermark within a PCD sub-band is indeed easy in the Fourier domain [START_REF] Autrusseau | A robust image watermarking technique based on quantization noise visibility thresholds[END_REF], however, as previously emphasized, the wavelets sub-bands do not overlap entirely in the PCD sub-bands (Figure 1). One solution to best match the PCD, would be to embed only in the intersection of both sub-bands. A watermarking technique was designed in the DWT domain, in order to confirm the mask's efficiency to others embedding domains. A three stages wavelet transform (9/7 filters) was applied on the input image, a noise-like watermark was embedded independently in the HH 2 or LH 1 . The detection technique remains the same as the one previously presented in [START_REF] Autrusseau | A robust image watermarking technique based on quantization noise visibility thresholds[END_REF]: cross-correlation is computed between the stored watermark and the extracted DWT coefficients.

EXPERIMENTAL RESULTS

Figure 2-a shows the spatial representation of the weighted watermark along with its Fourier representation (2-b), where the PCD is superimposed to explicitly show that most of the watermark's energy is contained within the appropriate visual sub-bands. It is important to notice on Figure 2-a that the proposed JND mask do not allow a strength adaptation in high activity areas, the watermark's strength is not increased in the image edges or textures. The JND mask being composed of weighted low frequencies (eq. 2) it provides a global watermark weighting parameter. Figure 2-c shows the variance of the watermark for each PCD sub-band (computed in the Fourier domain). This plot confirms that the frequency representation of the DWT watermark is indeed maintained in the visual sub-bands (peaks at positions 6, 7 and 11, see sub-bands numbering in Figure 2-a), and thus, the ∆C value of sub-band (IV, 1) can be used here for the mask implementation. As previously explained in section 3, the normalized cross-correlation is computed between the wavelet sub-band of the marked image, and the wavelet representation of the weighted watermark. Stirmark attacks were used to assess the robustness of the DWT embedding technique. Figure 3-b represents the cross correlation max value (Y-axis) as a function of 40 selected Stirmark attacks (X-axis). Unlike the Fourier domain embedding technique [START_REF] Autrusseau | A robust image watermarking technique based on quantization noise visibility thresholds[END_REF], the watermark is more widely spread into the frequency domain, and thus, the weighting parameter (K i,j in eq. 3) is sensibly lower and so are the correlation peaks (Y-axis in Figure 3-b). However, false positive and false negatives have been computed and the optimal detection threshold is set to 0.075 (dashed vertical line in Figure 3-a). The detection rate was found to be 60% for watermarks in sub-band HL 1 (dashed lines) and 80% for watermarks in HH 2 (solid lines). Effectively, as less coefficients are marked in HH 2 , the strength is increased, and so is the robustness. The main goal of this work is not to propose a full DWT domain watermarking technique, but rather to adapt a spatial JND based on complex HVS properties into the wavelet domain and ensure the usability of the mask. Table 2 summarizes the quality assessment of the proposed embedding technique compared to previous works operating in the Fourier domain [START_REF] Autrusseau | A robust image watermarking technique based on quantization noise visibility thresholds[END_REF]. PSNR, wPSNR, SSIM and C4 were used (refer to [START_REF] Marini | Evaluation of standard watermarking techniques[END_REF] for details and performances of the quality metrics). We can notice on Table 2 that the quality requirements are fulfilled, all metrics present very good results. Most quality measures appeared to be better than the ones in [START_REF] Autrusseau | A robust image watermarking technique based on quantization noise visibility thresholds[END_REF], where subjective experiments were used to validate the objective scores. This slight improvement in objective quality measures is due to a wider spreading of the watermark in the Fourier frequencies. However, the wider is the watermark, the lower might be its strength.

Thus the robustness performances in the DWT domain might be slightly lower than in [START_REF] Autrusseau | A robust image watermarking technique based on quantization noise visibility thresholds[END_REF]. For SSIM and C4, the closer the measure is to one, the best is the quality. Table 2 shows results for watermarked image in the HL 1 and HH 2 sub-bands. For comparison purpose, the bottom value gives results for the watermarking technique operating in the Fourier domain (8 × 8 watermark in sub-band (IV, 1)).

CONCLUSION

Although invisibility is a strong requirements for digital image watermarking, very few JND masks take into account advanced HVS features. Several perceptual masks based on heuristic considerations are proposed, but in order to make the embedding possible in the whole input image, including the smooth areas, complex HVS models must be considered. We proposed here an adaptation of a spatially defined JND mask into the wavelet domain. The proposed mask model have been evaluated in terms of both invisibility and robustness capabilities and showed good results regarding both requirements.

Figure 1 -Fig. 1 .

 11 Fig. 1. Superimposition of perceptual channel decomposition and wavelet tiling of the spectrum.

Fig. 2 .

 2 Fig. 2. (a) Spatial representation of a HL 1 watermark along with (b) its Fourier representation and (c) variance within each of the 17 PCD sub-bands.

  Cross-correlation peaks for 40 attacks.

Fig. 3 .

 3 Fig. 3. Detection performances.

Table 1 .

 1 Table 1 below gives the values provided by observers during subjective experiments. Experimental ∆C i,j for every PCD sub-band.

			radial selectivity
	angular selectivity	LF	III	IV	V
	1	0.5	0.0034	0.0066	0.026
	2		0.004	0.010	0.04
	3		0.0034	0.010	0.04
	4		0.004	0.0066	0.026
	5			0.010	0.04
	6			0.010	0.04

3. PERCEPTUAL WATERMARKING

As previously explained, we can define the visibility of quantization noise in visual channel content for complex signals. We will now exploit this property in a watermarking context for the strength determination process. Since this model operates in a psychovisual space, the input image has to be converted into luminances (in Cd/m 2 ) depending on the display. The monitor's "gamma function" is used to transform the digital grey level values N (m, n) into the photometric quantity known as luminance: L(m, n)=L min + L max ×

Table 2 .

 2 Quality assessment within two distinct levels.

	Image		PSNR	SSIM	wPSNR	C4
		HL1	49.3	0.991	50.7	0.978
	lena	HH2	52.8	0.996	54.1	0.974
		from [2]	48.7	0.964	49.9	0.974
		HL1	53.2	0.996	54.9	0.949
	boats	HH2	52.8	0.995	54.5	0.951
		from [2]	45.2	0.980	46.9	0.985
		HL1	51.2	0.997	53.2	0.946
	goldhill	HH2	52.7	0.998	54.6	0.931
		from [2]	47.0	0.993	48.7	0.974
		HL1	53.2	0.996	54.9	0.986
	baboon	HH2	52.8	0.999	57.7	0.977
		from [2]	42.8	0.992	47.5	0.978
		HL1	53.1	0.996	54.3	0.910
	kodie	HH2	52.8	0.996	54.1	0.904
		from [2]	51.8	0.987	53.4	0.980


P-Y Henry

L Bousset

P Sourrouille

P Jarne

Partial selfing, ecological disturbance and reproductive assurance in an invasive freshwater snail

Keywords: demography, density, mating system, Physa acuta, reproductive assurance, self-fertilization

Although reproductive assurance (RA) might play a central role in the evolution of the selfing rate, this hypothesis has never been seriously investigated in an hermaphroditic animal. We studied the mating system of the freshwater snail Physa acuta in which the availability of mating partners might be highly variable, because this species is an efficient colonizer occupying unstable habitats. A total of 11 populations differing in ecological disturbance regime (water level, openness) and snail densities were monitored over 2 years. The outcrossing rate was estimated in ca 10 families per population using microsatellite markers and the progenyarray approach. Components of fecundity and survival were recorded for each progeny. Predominant outcrossing (tm=0.94) was detected, with a few individuals (4%) purely selfing. The outcrossing rate did not explain among-family variation in fitness components. None of the predictions formulated under the RA hypothesis were verified: (i) selfing was related neither to disturbed habitats, nor to temporal density fluctuations, (ii) it was positively related to population density, (iii) it co-occurred with multiple paternity, and (iv) it did not induce delayed reproduction. Explanations for these negative results are discussed in light of other arguments supporting the RA hypothesis in P. acuta, as well as alternative theories explaining the occurrence of partial selfing, as either a genetically fixed or plastic trait.

Introduction

The capacity to self-fertilize is obviously advantageous: offspring can be produced even when mates are unavailable. Reproductive assurance (RA) has therefore been a popular hypothesis for explaining the evolution of the selfing rate (see [START_REF] Jarne | The evolution of the selfing rate in functionally hermaphrodite plants and animals[END_REF]. Further investigations (eg [START_REF] Lande | The evolution of self-fertilization and inbreeding depression in plants. I. Genetic models[END_REF] rationalized the understanding of the genetical processes driving the evolution of the selfing rate, converting the problem mainly into a trade-off between the automatic advantage (50%) of gene transmission under selfing and the lowered fitness of selfed relative to outcrossed progeny, that is, inbreeding depression. The influence of environmental constraints has since then been incorporated into several theoretical models that account for population structure and dynamics (eg [START_REF] Ronfort | A stochastic model of selection on selfing rates in structured populations[END_REF][START_REF] Pannell | Effects of population size and metapopulation dynamics on a mating-system polymorphism[END_REF][START_REF] Cheptou | The evolution of self-fertilization in density-regulated populations[END_REF]. Theoretical studies focusing specifically on the role of the reproductive assurance hypothesis (RAH) showed that delayed selfing is always advantageous when it occurs as RA [START_REF] Lloyd | Self-and cross-fertilization in plants. II. The selection of self-fertilization[END_REF]Tsitrone et al, 2003a), even in a metapopulational context [START_REF] Pannell | Baker's law revisited: reproductive assurance in a metapopulation[END_REF]; but see [START_REF] Pannell | Effects of population size and metapopulation dynamics on a mating-system polymorphism[END_REF].

The wide acceptance of the central role of RA in the evolution of selfing (eg [START_REF] Lloyd | Demographic factors and mating patterns in Angiosperms[END_REF][START_REF] Schoen | How does self-pollination evolve? Inferences from floral ecology and molecular genetic variation[END_REF][START_REF] Barrett | The evolution of plant sexual diversity[END_REF] basically relies on empirical evidence. Major arguments are the over-representation of self-compatible species among long-distance colonizing species (reviewed in [START_REF] Pannell | Baker's law revisited: reproductive assurance in a metapopulation[END_REF], the repeated evolution of selfing populations or species from outcrossing ancestors in constraining habitats (eg [START_REF] Dole | Reproductive assurance mechanisms in three taxa of the Mimulus guttatus complex (Scrophulariaceae)[END_REF], and the occurrence of autogamous characters in populations experiencing chronically low density (eg [START_REF] Dole | Reproductive assurance mechanisms in three taxa of the Mimulus guttatus complex (Scrophulariaceae)[END_REF][START_REF] Städler | Population biology, genetic structure, and mating system parameters in freshwater snails[END_REF]. More direct arguments are derived from studies on plant pollination showing that selfing correlates with pollen-vector unreliability and/or low plant population density (eg [START_REF] Lloyd | Demographic factors and mating patterns in Angiosperms[END_REF][START_REF] Watkins | Outcrossing rates as related to plant density in Phlox drummondii[END_REF][START_REF] Fausto | Reproductive assurance and the evolutionary ecology of selffertilization in Clarkia xantiana (Onagraceae)[END_REF]. However, the evolutionary significance of the RAH has been recently questioned (eg [START_REF] Barrett | The evolution of plant sexual diversity[END_REF], the advantage of selfing being outweighed by seed discounting and inbreeding depression [START_REF] Herlihy | Genetic cost of reproductive assurance in a self-fertilizing plant[END_REF]. Moreover, although widely studied in plants, RA has hardly been investigated in hermaphroditic animals.

We engaged in an ecological and genetical study aiming at evaluating the role of RA in a hermaphroditic freshwater snail, Physa acuta (Gastropoda: Pulmonata; [START_REF] Henry | Dynamique des populations et évolution du système de reproduction. Etude chez le gastéropode hermaphrodite Physa acuta (Pulmoné)[END_REF]Tsitrone et al, 2003b; Comment citer ce document : Henry, P.-Y., Bousset, L., Sourrouille, P., Jarne, P. (2005). Partial selfing, ecological disturbance and reproductive assurance in an invasive freshwater snail. Heredity, 95, 428-436. , DOI : 10.1038/sj.hdy.6800731

This file is the post-print of: Henry et al. Heredity (2005) 95:428-436 Doi: 10.1038/sj.hdy.6800731 2 [START_REF] Bousset | Population biology of the invasive freshwater snail Physa acuta approached through genetic markers, ecological characterization and demography[END_REF]. Freshwater snails are good candidates for such a work [START_REF] Jarne | The evolution of the selfing rate in functionally hermaphrodite plants and animals[END_REF][START_REF] Städler | Population biology, genetic structure, and mating system parameters in freshwater snails[END_REF]. P. acuta predominantly outcrosses, probably because of strong inbreeding depression [START_REF] Jarne | The influence of self-fertilization and grouping on fitness attributes in the freshwater snail Physa acuta: population and individual inbreeding depression[END_REF][START_REF] Henry | Environment-dependent inbreeding depression in a hermaphroditic freshwater snail[END_REF]Tsitrone et al, 2003b). Residual selfing in outcrossing animal species might constitute an adaptive strategy of RA (Tsitrone et al, 2003a), and there are some ecological reasons for thinking that it might be so in P. acuta. First, this species inhabits a wide range of more or less temporary habitats varying in openness (eg ponds versus rivers). As a consequence, population size is likely to vary dramatically through time because of drying out events or floods (reviewed in [START_REF] Städler | Population biology, genetic structure, and mating system parameters in freshwater snails[END_REF]. Second, P. acuta is the commonest species in isolated habitats (ie ponds) of the area studied here (reviewed in [START_REF] Henry | Dynamique des populations et évolution du système de reproduction. Etude chez le gastéropode hermaphrodite Physa acuta (Pulmoné)[END_REF], in accordance with its invasive status [START_REF] Dillon | Populations of the European freshwater pulmonate Physa acuta are not reproductively isolated from American Physa heterostropha or Physa integra[END_REF]Bousset and Jarne, unpublished data). As colonization has been shown to be associated with founding events in freshwater snails [START_REF] Charbonnel | A spatial and temporal approach to microevolutionary forces affecting population biology in the freshwater snail Biomphalaria pfeifferi[END_REF][START_REF] Meunier | Field and experimental evidence of preferential selfing in the freshwater mollusc Lymnaea truncatula (Gastropoda, Pulmonata)[END_REF], this opens opportunities for selfing as a RA strategy. In this context, Tsitrone et al (2003a) have shown theoretically that a key-factor, if RA is to play a role in the evolution of selfing, is the waiting timethe time that an isolated individual should wait before selfing when compared to an outcrossing individual. Its occurrence has been experimentally documented, and this character proved to be heritable (Tsitrone et al, 2003b;Escobar, David and Jarne, unpublished data). This direct approach suggests that the RAH can explain residual selfing rates in outcrossing species.

The RAH is here tested through an indirect, correlative approach, aiming at determining whether population and individual variation in selfing rate are affected by environmental, demographic and individual factors as predicted under the RAH. In all, 11 sites harbouring populations of P. acuta were studied representing four habitat types characterized by their water regime (temporary versus permanent) and openness (ponds versus rivers). Population density and water level were monitored monthly for 2 years. The outcrossing rate was estimated in ca 10 families per population using microsatellite markers and the progenyarray approach [START_REF] Ritland | Extension of models for the estimation of mating systems using n independent loci[END_REF].

Predictions from the RAH were formulated under two models for the genetic basis of selfing. First, selfing was assumed to be a genetically fixed trait, with natural selection acting on the selfing rate. The selfing rate would be adapted to the past history of mate availability. The RAH predicted that mate availability and the outcrossing rate covary across ecological disturbance regimes (water regime and openness), both being lower in temporary than in permanent habitats, and in ponds than in rivers. Second, selfing was assumed to be a genetically plastic trait, with natural selection acting on the adjustment of the selfing rate to environmental conditions (which is possibly the case in P. acuta; Tsitrone et al, 2003b). The selfing rate would thus be adapted to the instantaneous mate availability. Predictions from the RAH were: the outcrossing rate increases with increasing density, multiple paternity is positively associated with the outcrossing rate, and the time to first egg laying of isolated individuals under laboratory conditions is negatively associated with the outcrossing rate.

Materials and methods

Populations studied and progeny arrays

A total of 11 populations of P. acuta were monitored monthly from July 2000 to 2002 in the Montpellier area, France. Ecological details of these populations are given in Table 1 (see also [START_REF] Henry | Dynamique des populations et évolution du système de reproduction. Etude chez le gastéropode hermaphrodite Physa acuta (Pulmoné)[END_REF]. Analysis of population genetic structure is presented in [START_REF] Bousset | Population biology of the invasive freshwater snail Physa acuta approached through genetic markers, ecological characterization and demography[END_REF]. Populations were distinguished according to habitat openness (rivers=open versus ponds=closed) and water regime (permanent versus temporary). Temporary sites dried out at least once during the 2-year survey. The snail density was visually scored (ca 0, 1, 5, 10, 25, 50 or 100 ind/m 2 ) at each visit. When snails were spatially aggregated, we reported the range of density (dsampl in Table 1), and the range mean was used in statistical analyses. The ecological variables are habitat openness (Open versus Closed), water regime (Permanent versus Temporary), density at sampling (dsampl in individual/m 2 ; intervals are given when P. acuta were spatially aggregated), harmonic mean of population densities (dharm), and frequency of occurrence of population densities of one ind/m 2 or less (dfreq). A population covariate is the proportion of fertile individuals (pfert). The individual covariates are the size of G0 individuals (mm), number of eggs per capsule (Negg), egg-laying delay (tlaying in days), incubating time (tincub in days), and survival 7 days after first hatching (S7 d). Nfam, NG0 and NG1 are the number of families, of G0 and G1 individuals genotyped. Fis is the inbreeding coefficient, and tFis the population outcrossing rates based on Fis. tm and ts are the multilocus and single-locus outcrossing rates, and tm-ts the amount of apparent selfing accounted for by biparental inbreeding. rs and rp are the correlations of selfing among families and of outcrossed paternity within families. Stars indicate estimates significantly higher than 0 (Fis, tm-ts, rs) or lower than 1 (tm, ts, rp).

In July 2000, 20 mature individuals (denoted G0) were collected in each population. They were measured (length) and isolated in 75 ml plastic boxes less than 6 h after field collection. They were maintained in the laboratory at 25°C under a 12:12 photoperiod regime and fed ad libitum with boiled lettuce. Water was changed every 3 days. Plastic boxes were checked daily for egg capsules, until the first capsule was laid, over a period of 21 days. The egg-laying delay (tlaying; mean: 5.1±3.9 days) was calculated as the time elapsed between isolation and first egg laying. The fraction of fertile individuals (pfert) was the fraction of individuals that laid eggs over these 21 days. When several capsules were laid by a given snail, only the oldest capsule, based on embryo development, was analysed. Its number of eggs (Negg) was recorded. Hatching was checked daily until the first juvenile (denoted G1) was detected, from which we derived the minimum incubation time (tincub), that is, the time elapsed between egg laying and first hatching. Survival (S7 d) was defined as the number of surviving G1 7 days after first hatching divided by the number of eggs per capsule.

A G0 mother and its G1 offspring constituted a family. Since P. acuta stores sperm after copulation (eg [START_REF] Wethington | Selfing, outcrossing, and mixed mating in the freshwater snail Physa heterostropha: lifetime fitness and inbreeding depression[END_REF], preferentially outcrosses [START_REF] Jarne | The influence of self-fertilization and grouping on fitness attributes in the freshwater snail Physa acuta: population and individual inbreeding depression[END_REF]Tsitrone et al, 2003b;[START_REF] Bousset | Population biology of the invasive freshwater snail Physa acuta approached through genetic markers, ecological characterization and demography[END_REF]) and G0 individuals were placed in optimal conditions for egg laying, we assumed that the estimated outcrossing rate reflects that under natural conditions.

Estimating the outcrossing rate G0 and G1 individuals were stored at -80°C and in 95% alcohol, respectively, prior to DNA extraction. DNA extraction, PCR amplification and allele size scoring were performed according to [START_REF] Sourrouille | Microsatellite variation in the freshwater snail Physa acuta[END_REF]. The seven loci used are described in [START_REF] Monsutti | Dinucleotide microsatellite loci reveal a high selfing rate in the freshwater snail Physa acuta[END_REF] and [START_REF] Sourrouille | Microsatellite variation in the freshwater snail Physa acuta[END_REF]. Loci were coamplified by PCR (Pac1, Pac2 and Pac5; Pac4 and Pac7; AF108762 and AF108764; [START_REF] Bousset | Population biology of the invasive freshwater snail Physa acuta approached through genetic markers, ecological characterization and demography[END_REF] and products were run simultaneously on ABI Prism Genetic Analysers (310 or 3100). PCR amplification and allele size scoring were performed simultaneously for the G0 mother and its progeny. In all, 10 families per population and seven juveniles per family were genotyped. When less than seven juveniles were available, at least one family was added in order to genotype the same overall number of juveniles per population (Table 1). Families and juveniles were randomly chosen. Mendelian segregation of alleles at the loci used was checked elsewhere [START_REF] Bousset | Population biology of the invasive freshwater snail Physa acuta approached through genetic markers, ecological characterization and demography[END_REF]. When null alleles were detected (loci Pac4, Pac7), data from these families and loci were omitted from the analysis. Estimates of allelic frequencies per locus and per population This file is the post-print of: Henry et al. Heredity (2005) 95:428-436 Doi: 10.1038/sj.hdy.6800731 4 are given in Web Table 1, available online. The mean number of polymorphic loci per population ranged between 3 and 6, and the mean number of alleles per loci and population varied between 1.8 and 4.2 [START_REF] Bousset | Population biology of the invasive freshwater snail Physa acuta approached through genetic markers, ecological characterization and demography[END_REF]).

The mating system parameters were estimated using maximum likelihood procedures (software MLTR 2.4; [START_REF] Ritland | Extension of models for the estimation of mating systems using n independent loci[END_REF]. Both population-level multilocus (tm) and single-locus (ts) outcrossing rates were estimated. Family-level multilocus outcrossing rates are noted tmf. The difference between tm and ts gives a minimum estimate of apparent selfing due to biparental inbreeding [START_REF] Ritland | Extension of models for the estimation of mating systems using n independent loci[END_REF]. Correlations of selfing (rs) document the among-family variation in selfing, with a non-null correlation indicating significant heterogeneity [START_REF] Ritland | Extension of models for the estimation of mating systems using n independent loci[END_REF]. Correlations of paternity (rp) represent the proportions of full-sibs (versus halfsibs) among outcrossed offspring [START_REF] Sun | Mating system of yellow starthistle (Centaurea solstitialis), a successful colonizer in North America[END_REF][START_REF] Ritland | Extension of models for the estimation of mating systems using n independent loci[END_REF] with rp significantly lower than 1 indicating multiple paternity. The Expectation-Maximization method and the 'parent chosen in probability' option were used. SD of estimates were calculated based on 500 bootstraps assuming Normal distributions. The resampling unit was, respectively, the family and the individual within family for population-and familylevel estimates. SD for averages over all population or family estimates were computed with the δ-method [START_REF] Mood | Introduction to the Theory of Statistics[END_REF]. Fis, tm-ts and rs were considered as significantly higher than 0, and tm, ts and rp as lower than 1 when more than 95% of the 500 resampled estimates were, respectively, higher than 0 or lower than 1. Assuming that selfing is the sole source of inbreeding, we calculated the selfing rate corresponding to Fis values according to: tFis=(1-Fis)/(1+Fis) [START_REF] Pollak | On the theory of partially inbreeding finite populations. I. Partial selfing[END_REF].

Selfing rates and fitness traits

Since juveniles were genotyped 7 days after hatching, tm and tmf estimates are expected to be biased positively, because of inbreeding depression [START_REF] Jarne | The evolution of the selfing rate in functionally hermaphrodite plants and animals[END_REF]. This is particularly likely since inbreeding depression is strong in juvenile P. acuta [START_REF] Wethington | Selfing, outcrossing, and mixed mating in the freshwater snail Physa heterostropha: lifetime fitness and inbreeding depression[END_REF][START_REF] Jarne | The influence of self-fertilization and grouping on fitness attributes in the freshwater snail Physa acuta: population and individual inbreeding depression[END_REF][START_REF] Henry | Environment-dependent inbreeding depression in a hermaphroditic freshwater snail[END_REF]Tsitrone et al, 2003b). The tm value should therefore explain population-level variation in the proportion of fertile G0 individuals (pfert), and tmf should explain family-level variation in the number of eggs per capsule (Negg), incubating time (tincub) and survival at 7 days (S7 d).

Tests of the RAH

Predictions from the RAH depend on the assumed genetic basis of the selfing rate (see Introduction; Table 2).

(i) The selfing rate could be a genetically fixed trait (predictions prefixed by F) evolving under the selective pressure of temporal fluctuation in mate availability. Mate availability was expected to differ with habitat. Episodes of low density should be less frequent in permanent than in temporary habitats, because population bottlenecks associated with drying-out events are expected to be more frequent in temporary habitats [START_REF] Charbonnel | A spatial and temporal approach to microevolutionary forces affecting population biology in the freshwater snail Biomphalaria pfeifferi[END_REF][START_REF] Trouvé | Evolutionary implications of a high selfing rate in the freshwater snail Lymnaea truncatula[END_REF][START_REF] Bousset | Population biology of the invasive freshwater snail Physa acuta approached through genetic markers, ecological characterization and demography[END_REF]. They should also be more frequent in closed than in open habitats, because immigration more easily refills populations after a population bottleneck in open than in closed habitats [START_REF] Bousset | Population biology of the invasive freshwater snail Physa acuta approached through genetic markers, ecological characterization and demography[END_REF]. The probability of low density episodes was estimated based on the harmonic mean density (dharm) and the proportion of visits when density was one individual per m 2 or less (dfreq) over our 2-year survey. This threshold was chosen because we were technically unable to estimate density precisely at lower values. In response to these expected variations of mate availability, the outcrossing rate tm was predicted to be higher in permanent than in temporary habitats (F1), in open than in closed habitats (F2), positively correlated to dharm (F3), and negatively correlated to dfreq (F4).

(ii) The selfing rate could be a plastic trait adapted to realise the best compromise between inbreeding depression and instantaneous mate availability (predictions prefixed by P). (P1) At the population level, the probability of outcrossing should be positively correlated with partner availability, estimated as local density during genetic sampling (dsampl). (P2) Under delayed selfing, partner availability should not only affect the fraction of selfed ovules, but also the number of outcrossing partners. The probability of outcrossing should thus be negatively correlated to the proportion of half-sibs among outcrossed offspring, rp. (P3) Previous work has shown the existence of a waiting time in P. acuta, defined as the difference in age at first egg-laying between outcrossing and selfing virgin individuals (Tsitrone et al, 2003a(Tsitrone et al, , 2003b)). By analogy, the individual egg-laying delay (tlaying) was expected to negatively correlate with the outcrossing rate.

Statistical analyses

The probability of outcrossing was analysed using logistic regression models, followed by likelihood ratio tests (LRT). To convert the normal ML estimates of outcrossing rate tm into a binomial variable, that is, the probability of being an outcrossed offspring, the number of genotyped progeny per population NG1 was multiplied by tm and 1-tm (Table 1). This produced the estimated number of outcrossed and selfed juveniles per population. Using logistic regressions, with tm as the dependent variable and the reconstructed number of juveniles as sample size, accounted for the binomial variance of tm estimates. Note that this is not the case under analyses of variance (ANOVA) with classical transformation (eg arcsin) of tm which, however, provided similar results (results not shown). The explanatory variables were habitat openness, water regime and density variables (dharm, dfreq and dsampl). The effects of openness, water regime and their interaction were tested in a single model, and that of density variables in separate models. The dependent variables dfreq, pfert and S7 d were also analysed using logistic regression models: dfreq was modelled in relation to openness, water regime and their interaction, and pfert in relation to tm. Familylevel Ln(tlaying+1), Ln(Negg), and Ln(tincub) were analysed using ANOVAs. The explanatory variables for S7 d, Ln(tlaying+1), Ln(Negg) and Ln(tincub) were the log-transformed individual size, tlaying+1, Negg, tincub and tmf (see Table 3 for more details on each model). The homogeneity of dharm among habitats was tested using a Kruskal-Wallis test, followed by Mann-Whitney tests of the predicted relationships with habitat openness and water regime. Logistic regression and ANOVA models were built with JUMP 3.2.1 (S.A.S. Institute Inc.), and nonparametric tests with Statistica 5.5 (Statsoft Inc.).

Table 3: Analysis of fitness components and outcrossing rate in P. acuta

The table reports the F-values and their associated probability (P). The dependent variables (first row) are the logtransformed egg-laying delay (tlaying; since individuals laid eggs within 1 day, the dependent variable was Ln(tlaying+1)), number of eggs per capsule (Negg), minimum incubating time (tincub) and survival at 7 days (S7 d). The first three variables were analysed with ANOVA models and S7 d with a logistic regression model. The explanatory variables were the logtransformed individual size, Ln(Negg), Ln(tlaying+1) and the family-level multilocus outcrossing rate tmf. This file is the post-print of: Henry et al. Heredity (2005) 95:428-436 Doi: 10.1038/sj.hdy.6800731 6

Results

Mating system characterization P. acuta preferentially outcrossed in all populations. Population-level multilocus outcrossing rates ranged from 0.731 to 0.998 (average: 0.936±0.019). The selfing rate significantly differed from 0 in three populations out of 11 (Lam12, Vio7, Vio1; Table 1). However, the correlation of selfing (rs) among families indicated that outcrossing rates significantly varied among families in 10 populations out of 11 (Table 1). The family-level outcrossing rates ranged from 0.001 to 1 (average: 0.885±0.002). Partial selfing was detected in at least one family per population (except in Mos6), and purely selfed families occurred in three populations only (Lam12, Vio7 and Vio1). tm-ts ranged from -0.042 to 0.12 (average: 0.041±0.011) with significant biparental inbreeding occurring in five populations (significant tm-ts in Table 1). Fis differed from zero in four populations out of 11 (range: 0.088-0.210). However, no obvious relationship could be established with inbreeding since Fis was significant in Vio2 and Vio8 where neither selfing, nor biparental inbreeding was detected. Population-level correlations of paternity (rp) ranged from 0.028 to 0.629 (Table 1), indicating multiple paternities within progenies in all populations but Vio2. The family-level average was 0.260±0.613 (range: 0.001-1) corresponding to a mean number of 3.8 fathers per family (ie 1/rp; [START_REF] Sun | Mating system of yellow starthistle (Centaurea solstitialis), a successful colonizer in North America[END_REF].

Selfing rates and fitness traits

The logistic regression model for S7 d did not fit the data (χ 2 =576.3, d.f.=91, P<10 -4 ). Deviances were thus divided by a variance inflation factor ĉ (deviance/d.f.=6.333) before testing effects with LRTs. No effect of the outcrossing rate on any of the four reproductive parameters studied was detected among families (Table 3). This also holds for the proportion of fertile individuals (LRT=0.189, d.f.=1, P=0.664). None of the individual traits expected to confound the effect of the mating system significantly explained variation in fitness components (Table 3).

Tests of the RAH

Mate availability depended on the disturbance regime. The harmonic mean density (dharm) differed among the four habitats (H=7.833, d.f.=3, N=11, P=0.050) with slightly higher values in closed temporary than in other sites (Table 1). When habitat types were tested separately, there was no detectable effect for either habitat openness (U=6, Z=-1.643, P=0.100), or for water regime (U=13, Z=0.365, P=0.715). The probability of occurrence of densities of one individual or less per m 2 (dfreq) varied among habitats (interaction Openness*Water regime, LRT=29.818, d.f.=1, P<10 -4 ). However, it did not follow the predicted pattern: low densities were more frequent in open temporary (45%, N=42 visits) than in closed permanent (36%, N=72), closed temporary (10%, N=62) and open permanent sites (8%, N=40).

None of the predictions of the RAH (Table 2) were supported by our results. Although variation in outcrossing probability could be explained by habitat characteristics (interaction Openness*Water regime, Table 4), they did not fit with variation in dfreq, at variance with predictions from the RAH: the overall outcrossing rate was lower in closed permanent (0.86) than in open temporary (0.92), closed temporary (0.95) and open permanent (0.98). Moreover, this variation was no more significant when tested in the presence of dsampl (results not shown; F1 and F2 not fulfilled). dharm and dfreq did not explain the among-population variation in outcrossing probability (Table 4; F3 and F4 not fulfilled).

Table 4: Likelihood-ratio tests (LRT) for ecological correlates of the population-level outcrossing probability

The dependent variable was the reconstructed binomial data corresponding to the estimated multilocus outcrossing rate and sample size per population (see Materials and methods). The ecological covariates were openness (closed versus open), water regime (permanent versus temporary), the harmonic mean density (dharm), the frequency of occurrence of densities ≤1 ind m -2 (dfreq) and the density at sampling (dsampl). All effects were tested against the null model.

Comment citer ce document : Henry, P.-Y., Bousset, L., Sourrouille, P., Jarne, P. (2005). Partial selfing, ecological disturbance and reproductive assurance in an invasive freshwater snail. Heredity, 95, 428-436. , DOI : 10.1038/sj.hdy.6800731

This file is the post-print of: Henry et al. Heredity (2005) 95:428-436 Doi: 10.1038/sj.hdy.6800731 7

The population-level probability of outcrossing was negatively dependent on the density at sampling (Table 4; Figure 1), contrary to the predicted positive correlation (P1). This effect remained strongly significant when using a more conservative statistical approach under which the sample size per population was reconstructed such as to be proportional to the SD of tmp (results not shown). At the population level, multiple paternities occurred in the three populations in which significant selfing was detected (Table 1). At the family level, there was no relationship between tmf and rp (r=0.088, d.f.=98, P=0.385; P2 not fulfilled). The egg-laying delay did not depend on the family-level multilocus outcrossing rate (Table 3;prediction P3), although essentially selfing individuals tended to exhibit long egg-laying delays.

Figure 1

Relationship between population-level multilocus outcrossing rate (tm±SD) and population density at sampling (dsampl) in P. acuta (N=11 populations).

Discussion

The mating system of P. acuta P. acuta reproduced primarily through outcrossing, with a population-level average outcrossing rate of 0.936. Over all families, the average outcrossing rate was 0.885. These results are in agreement with those in a natural population from the same area (0.921 in [START_REF] Jarne | The influence of self-fertilization and grouping on fitness attributes in the freshwater snail Physa acuta: population and individual inbreeding depression[END_REF]. [START_REF] Wethington | Selfing, outcrossing, and mixed mating in the freshwater snail Physa heterostropha: lifetime fitness and inbreeding depression[END_REF] found a rather lower value (0.385), but this was under laboratory conditions with North-American individuals. Our estimates based on Fis (average tFis : 0.89) were close to those derived from progeny arrays.

Examining fitness components showed that the variation in the proportion of fertile individuals, number of eggs per capsule, incubation time, and offspring survival were not explained by variation in outcrossing rate. This might be interpreted as evidence of limited inbreeding depression. However, this would be surprising since strong inbreeding depression has been consistently detected in P. acuta [START_REF] Wethington | Selfing, outcrossing, and mixed mating in the freshwater snail Physa heterostropha: lifetime fitness and inbreeding depression[END_REF][START_REF] Jarne | The influence of self-fertilization and grouping on fitness attributes in the freshwater snail Physa acuta: population and individual inbreeding depression[END_REF][START_REF] Henry | Environment-dependent inbreeding depression in a hermaphroditic freshwater snail[END_REF]. Moreover, in the partially selfing population Lam12, Fis did not differ from zero, suggesting the occurrence of inbreeding depression between the juvenile and adult stages. Alternatively, this might be due to imprecision in family-level estimates of tmf which is not accounted for in our analyses. Discounting the estimation error on tmf blurs the acceptance of the null hypothesis. To circumvent this general problem, a comprehensive statistical framework is required that would allow comparative analyses of mating system parameters, and in which the contribution of individual and population covariates could be directly built in the Maximum Likelihood procedure and tested through comparisons of nested models.

Tests of the RAH

The RAH has been repeatedly invoked to explain the maintenance of selfing in sessile or poorly mobile organisms. Its role was suggested in freshwater snails because these organisms often inhabit unpredictable temporary habitats (eg [START_REF] Jarne | The evolution of the selfing rate in functionally hermaphrodite plants and animals[END_REF][START_REF] Städler | Population biology, genetic structure, and mating system parameters in freshwater snails[END_REF][START_REF] Charbonnel | A spatial and temporal approach to microevolutionary forces affecting population biology in the freshwater snail Biomphalaria pfeifferi[END_REF]. Tsitrone et al (2003a) This file is the post-print of: Henry et al. Heredity (2005) 95:428-436 Doi: 10.1038/sj.hdy.6800731 8 animal species (see Introduction). Their model has been supported by experimental work conducted in P. acuta that showed that the waiting time is heritable (Tsitrone et al, 2003b). There is therefore evidence from a direct approach that the RAH explains residual selfing in highly outcrossing species.

The RAH was here evaluated using an indirect, correlative approach. None of the tests performed (Table 2) support the hypothesis that P. acuta self-fertilizes in response to a low probability of mate encounter. At the population level, the ecological constraints on the mating system related to habitat, via population density, were not those expected: temporary or closed habitats were no more likely to go regularly through low-density periods than permanent or open habitats (the opposite trend was indeed found). A conservative interpretation of our results is that the outcrossing rate was not related to those environmental variables that might influence the probability of mate encounter, that is, water regime, habitat openness, temporal fluctuations of density, and occurrence of low-density events. The density of G0 individuals at sampling did explain a significant fraction of the among-population variation in outcrossing rate. However, the relationship was opposite to the prediction: increasing population density was associated with decreasing outcrossing rate. Another surprising result under the RAH is that multiple paternity and the selfing rate were positively associated: why did individuals self their own ovules when mate density was high enough to ensure multiple paternities in some offspring (see also Cheptou et al, 2002)? On the whole, the indirect approach developed here did not provide any evidence in favour of the RAH as a major determinant of the selfing rate. This is compatible with results from another indirect approach in which the probability of mate encounter was experimentally manipulated [START_REF] Henry | Dynamique des populations et évolution du système de reproduction. Etude chez le gastéropode hermaphrodite Physa acuta (Pulmoné)[END_REF].

It must be emphasized that our predictions on the probability of mate encounter were based on shortterm ecological trends and population dynamics (over a 2-year survey). This did not capture the occurrence of much rarer events of low mate encounter probability, that is, due to abnormal disturbance events (eg, prolonged drying out) or very rare founding events. Those events might have indeed been the ecological factors maintaining selfing as a RA mechanism. Founding events very likely occur in P. acuta since this species displays high dispersal abilities (reviewed in [START_REF] Henry | Dynamique des populations et évolution du système de reproduction. Etude chez le gastéropode hermaphrodite Physa acuta (Pulmoné)[END_REF] and is the most successful worldwide invasive freshwater snail [START_REF] Dillon | Populations of the European freshwater pulmonate Physa acuta are not reproductively isolated from American Physa heterostropha or Physa integra[END_REF]Bousset and Jarne, unpublished data). Such founding events might be scarce in the studied area, since neither demographic, nor genetic signatures were detected [START_REF] Henry | Dynamique des populations et évolution du système de reproduction. Etude chez le gastéropode hermaphrodite Physa acuta (Pulmoné)[END_REF][START_REF] Bousset | Population biology of the invasive freshwater snail Physa acuta approached through genetic markers, ecological characterization and demography[END_REF]). The mating system would rapidly revert to preferential outcrossing when mates become available, because of inbreeding depression [START_REF] Stephenson | Interrelationships among inbreeding depression, plasticity in the self-incompatibility system, and the breeding system of Campanula rapunculoides L. (Campanulaceae)[END_REF][START_REF] Pannell | Effects of population size and metapopulation dynamics on a mating-system polymorphism[END_REF], and signals of RA would no more be detectable.

Very few studies in animals have evaluated the RAH, and most failed. In predominantly selfing snails, [START_REF] Trouvé | Evolutionary implications of a high selfing rate in the freshwater snail Lymnaea truncatula[END_REF] and [START_REF] Charbonnel | The influence of genetic factors and population dynamics on the mating system of the hermaphroditic freshwater snail Biomphalaria pfeifferi[END_REF] concluded that the among-population variation in selfing rate is inconsistent with this hypothesis. The only supporting results come from comparative analyses showing that selfing gastropods are better colonizers [START_REF] Selander | The genetic structure of populations as illustrated by molluscs[END_REF][START_REF] Bengtsson | Do pioneers have r-selected traits? Life history patterns among colonizing terrestrial gastropods[END_REF]. Less controversial arguments come from studies in plants, where negative relationships between self-pollination (and morphological or incompatibility-related traits) and density (of plants or of pollination vectors) have been repeatedly documented (eg Rankin et al, 2002 and references therein; but see eg Cheptou et al, 2002 for opposite results). However, it is worth noting that most results are derived from correlative observations of fitness components (the indirect approach above). Moreover, the only satisfactory quantitative test dismissed the evolutionary importance of RA because the advantage of autogamy was outweighed by seed discounting and inbreeding depression [START_REF] Herlihy | Genetic cost of reproductive assurance in a self-fertilizing plant[END_REF]. On the whole, the significance of selfing as a RA strategy, particularly in animals, remains to be demonstrated.

So, why do P. acuta and other snail species partially self-fertilize?

Although several studies have dealt with the reproductive biology of P. acuta, the maintenance of some selfing remains unexplained [START_REF] Henry | Dynamique des populations et évolution du système de reproduction. Etude chez le gastéropode hermaphrodite Physa acuta (Pulmoné)[END_REF]; but see Tsitrone et al, 2003b). We briefly review several theories, which are alternatives to the RAH, which may be consistent with the empirical evidence from P. acuta. These explanations can be extended to other freshwater snail species with a mating system similar to that of P. acuta (see [START_REF] Henry | Environment-dependent inbreeding depression in a hermaphroditic freshwater snail[END_REF] for a list). A first view assumes that partial selfing in P. acuta is the expression of a mixed-mating system, that is, a polymorphism with genetically fixed selfing rates. Three evolutionary scenarios can be envisaged. (i) The selfing rate greatly varies among families [START_REF] Jarne | The influence of self-fertilization and grouping on fitness attributes in the freshwater snail Physa acuta: population and individual inbreeding depression[END_REF]this study;see Trouvé et al, 2003 This file is the post-print of: Henry et al. Heredity (2005) 95:428-436 Doi: 10.1038/sj.hdy.6800731 9 inbreeding depression [START_REF] Jarne | The influence of self-fertilization and grouping on fitness attributes in the freshwater snail Physa acuta: population and individual inbreeding depression[END_REF][START_REF] Henry | Environment-dependent inbreeding depression in a hermaphroditic freshwater snail[END_REF]Tsitrone et al, 2003b). The required conditions for a coevolution between loci determining the selfing rate and fitness driving to mixed mating therefore seem to be met [START_REF] Uyenoyama | Ecological and genetic factors directing the evolution of self-fertilization[END_REF]. However, the observation that no genetic covariation between selfing propensity and inbreeding depression was found in P. acuta (Escobar, David and Jarne, unpublished data) weakens this hypothesis. (ii) Biparental inbreeding, due to restricted dispersal in association with densitydependent recruitment [START_REF] Uyenoyama | Inbreeding and the cost of meiosis: the evolution of selfing in populations practicing biparental inbreeding[END_REF][START_REF] Ronfort | A stochastic model of selection on selfing rates in structured populations[END_REF], could also be relevant and was detected in five populations out of 11 in this study. However, significant selfing rates were not systematically associated with biparental inbreeding, and density-dependent recruitment, although likely, remains to be investigated in P. acuta. (iii) Mixed mating can evolve in structured populations (eg [START_REF] Pannell | Effects of population size and metapopulation dynamics on a mating-system polymorphism[END_REF]. The highly fragmented, potentially unstable populations of P. acuta, in association with preferential outcrossing, high dispersal abilities and high inbreeding depression, nicely fit the situation modelled by [START_REF] Pannell | Effects of population size and metapopulation dynamics on a mating-system polymorphism[END_REF]: in the case of extinction-colonization, preferential outcrossing should have a selective advantage over preferential selfing thanks to a higher number of viable propagules.

A second view envisages the selfing rate as a plastic trait, and plasticity itself as the target of natural selection, adjusted to environmental conditions [START_REF] Stephenson | Interrelationships among inbreeding depression, plasticity in the self-incompatibility system, and the breeding system of Campanula rapunculoides L. (Campanulaceae)[END_REF][START_REF] Pannell | Effects of population size and metapopulation dynamics on a mating-system polymorphism[END_REF]. Tsitrone et al (2003a) demonstrated that delaying selfing according to mate availability increases both the probability of mate encounter and future reproduction through resource reallocation towards growth. A side prediction is that, if the waiting time before selfing is positively related to survival, the propensity to self-fertilize should be higher when survival is low. If high population density represents a high mortality risk, this would explain why the outcrossing rate was here negatively related to population density. Also, it could explain withinindividual temporal variation, like age dependency, in selfing rate [START_REF] Watkins | Outcrossing rates as related to plant density in Phlox drummondii[END_REF][START_REF] Stephenson | Interrelationships among inbreeding depression, plasticity in the self-incompatibility system, and the breeding system of Campanula rapunculoides L. (Campanulaceae)[END_REF][START_REF] Rankin | Mating system instability in Schiedea menziesii (Caryophyllaceae)[END_REF], which remains to be documented in hermaphroditic animals.

Comment citer ce document : Henry, P.-Y., Bousset, L., Sourrouille, P., Jarne, P. (2005). Partial selfing, ecological disturbance and reproductive assurance in an invasive freshwater snail. Heredity, 95, 428-436. , DOI : 10.1038/sj.hdy.6800731

  Comment citer ce document : Henry, P.-Y., Bousset, L., Sourrouille, P.,Jarne, P. (2005). Partial selfing, ecological disturbance and reproductive assurance in an invasive freshwater snail. Heredity, 95, 428-436. , DOI : 10.1038/sj.hdy.6800731

Table 2 :

 2 Comment citer ce document : Henry, P.-Y., Bousset, L., Sourrouille, P., Jarne, P. (2005). Partial selfing, ecological disturbance and reproductive assurance in an invasive freshwater snail. Heredity, 95, 428-436. , DOI : 10.1038/sj.hdy.6800731This file is the post-print of:Henry et al. Heredity (2005) 95:428-436 Doi: 10.1038/sj.hdy.6800731 5 Seven predictions (numbered as in text) relating the outcrossing rate (t) to ecological characteristics, and fitness traits, as derived from the reproductive assurance hypothesis Predictions were derived under either fixed (noted F) or plastic (noted P) genetic determination of selfing. Habitats were characterized as permanent or temporary, and as open or closed. The outcrossing rate (and therefore the prediction) is either at the population, or at the family level. dsampl: density at sampling; dharm: harmonic mean of population densities; dfreq: frequency of occurrence of population densities of one ind/m 2 or less; rp: correlation of outcrossed paternity within families; tlaying: time elapsed between isolation and first egg laying.

  Interaction terms are denoted by a star.indicates that the effect was not included in the model. Comment citer ce document : Henry, P.-Y., Bousset, L., Sourrouille, P., Jarne, P. (2005). Partial selfing, ecological disturbance and reproductive assurance in an invasive freshwater snail. Heredity, 95, 428-436. , DOI : 10.1038/sj.hdy.6800731

  clarified the conditions under which RA is likely to play a role in mainly outcrossing Comment citer ce document : Henry, P.-Y., Bousset, L., Sourrouille, P., Jarne, P. (2005). Partial selfing, ecological disturbance and reproductive assurance in an invasive freshwater snail. Heredity, 95, 428-436. , DOI : 10.1038/sj.hdy.6800731

  and references therein for results in other snail species), as do selfing propensity and Comment citer ce document : Henry, P.-Y., Bousset, L., Sourrouille, P., Jarne, P. (2005). Partial selfing, ecological disturbance and reproductive assurance in an invasive freshwater snail. Heredity, 95, 428-436. , DOI : 10.1038/sj.hdy.6800731

Table 1 :

 1 Ecological characterization, population and individual covariates, and mating system parameters estimates (with SD in parentheses) for 11 populations of Physa acuta

	Comment citer ce document :
	Henry, P.-Y., Bousset, L., Sourrouille, P., Jarne, P. (2005). Partial selfing, ecological
	disturbance and reproductive assurance in an invasive freshwater snail. Heredity, 95, 428-436. , DOI
	: 10.1038/sj.hdy.6800731

Acknowledgements

We are very grateful to N Charbonnel, P-O Cheptou, R Williams and two anonymous reviewers for comments on the manuscript, to T Lenormand for statistical advices, and to the Service Commun des Marqueurs Génétiques at CEFE, especially C Debain, for technical support with the genetic analyses. This project was funded by grants to PJ from the Centre National de la Recherche Scientifique (Programme Jeune Equipe), the Ministère de l'Ecologie et du Développement Durable (Programme InvaBio) and the Bureau des Ressources Génétiques. PYH benefitted from a fellowship from the Ministère de l'Education Nationale, de la Recherche et de la Technologie, and LB from an ASC fellowship from the Institut National de la Recherche Agronomique.

Supplementary information

Web Table 1 Estimated allelic frequencies (output from the MLTR software) at seven microsatellite loci in the 11 populations of Physa acuta analysed. Values at loci Pac 4 and Pac7 were obtained once null homozygotes were discarded. Frequencies in Mos6 are raw data from G0 individuals because the small sample size prevented computing estimations. Alleles are designated by their size in base pairs.