O Sandre

C Ménager
email: menager@ccr.jussieu.fr

J Prost

V Cabuil

J-C Bacri

A Cebers

Shape transition of giant liposomes induced by an anisotropic spontaneous curvature

Keywords: numbers: 75.50 Mm, 87.16 Dg, 82.65 Dp, 82.70 Dd, 83.50 Pk

We explore how a magnetic field breaks the symmetry of an initially spherical giant liposome filled with a magnetic colloid. The condition of rotational symmetry along the field axis leads either to a prolate or to an oblate ellipsoid. We demonstrate that an electrostatic interaction between the nanoparticles and the membrane triggers the shape transition.

Phospholipidic bilayers are the basic structural element of most biological membranes. In water these bilayers can form closed vesicles also called liposomes. Due to this analogy with cells, liposomes are commonly used in biophysical studies. In spite of the obvious simplification it implies, artificial liposomes are used to mimic some biological phenomena in well defined, low component environments [1]. For example, analogies have been found between shapes of red blood cells (discocytes, stomatocytes, echinocytes, spherical ghosts) and shapes of synthetic liposomes. Most shapes have a spherical topology and an axis of rotational symmetry: spheres, prolate ellipsoids, dumbells, pear shapes, oblate ellipsoids, discocytes, stomacytes. Other shapes have a symmetry axis of finite order (starfishes), no symmetry axis (3 axed ellipsoids), or another topology (torus). The transformations between those shapes can be triggered by an external parameter: temperature [2,3], gravity [4], a pH gradient [5,6] or a photochemical reaction [7]. The area-differential-elasticity model (ADE) enables the attribution of all morphologies to domains in a single phase diagram [8]. The two control parameters are: (i) the spontaneous curvature c 0 which reflects the asymmetry between the two sides of the bilayer ; (ii) the 'reduced volume' v which quantifies the ratio of volume to area (high v for an inflated vesicle, low v for a flaccid one). It was shown in [9] that a vesicle can move from one equilibrium shape to another by subtle changes of c 0 , i.e. by perfusion of different sugars (in osmotic equilibrium). In reference [10] the authors theoretically consider asymmetric bilayers in contact with colloids, but they do not include any angular anisotropy.

In this article, we describe an anisotropic spontaneous curvature due to magnetic nanoparticles placed in the interior of giant liposomes. An aqueous solution of those particles is an ionic ferrofluid [11]. The liposomes filled with such a ferrofluid are called 'magnetoliposomes'. Here, we observe magnetoliposome shapes under a magnetic field of low intensity (400 Oe). In some cases liposomes elongate in the direction of the magnetic field (Fig. 1(a)), while in others they are compressed at their poles (Fig. 1(b)). In the following, we describe the relevant experimental parameters determining the shape: the concentration of magnetic nanoparticles inside the liposome (C p ), the salt concentration of the medium (C s ), the radius of the liposome (R 0 ) when spherical (with no field) and its axes a and b under deformation. Qualitatively, oblate deformations occur for small, weakly magnetized liposomes, at low ionic strengths. On the other hand, prolate deformations are observed for large, highly magnetized liposomes in salty media. In the last section we propose a model to explain this transition between prolate and oblate shapes.

Preparation of the magnetoliposomes

The phospholipid constituting the membrane is 1-2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) purchased as powder from Avanti Polar Lipids. The process used to encapsulate the ferrofluid inside the liposomes with maximal efficiency differs slightly from usual ones: the first step is to pre-hydrate the phospholipid film with the colloidal solution of ferrofluid, and then swell it with pure water. A small amount of perfectly dry powder (around 1mg) of DOPC is mixed with 10µL of the aqueous dispersion (pH 7) of magnetic nanoparticles (γFe 2 O 3 ) and spread and sheared with a glove finger on a glass support (Petri dish) to obtain an oily orange film. This film is presumably a lamellar phase swelled with charged particles. Immediately following the shearing, 1 mL of distilled water is poured onto the film of fat to start the spontaneous swelling of liposomes. Samples are placed in a water bath at 45°C and observed after half an hour with an optical microscope (Leica 40x, NA 0.65). Pictures from a CCD camera are digitized with a frame-grabber (LG-3, Scion Corp., Frederick MD). Most of the liposomes prepared this way are quasi-spherical with diameters ranging from 10 to 100 µm. Their interiors are orange-colored and their membranes exhibit thermal fluctuations.

Measurement methods

There is heterogeneity in the colloid entrapment yields of different liposomes. First, the film pre-hydrated with the colloid is mixed manually and hence the particles are not dispersed homogeneously. In addition, liposomes in the same sample can also differ in their 'histories': their membranes close at different times during the swelling process, which competes with the simple diffusion of the ferrofluid by water (without entrapment in a closed membrane). Thus, the concentration C p of magnetic nanoparticles encapsulated in a given liposome must be measured independently. Therefore we make magnetophoresis experiments. Briefly, the method consists in applying a controlled magnetic field gradient. The intensity of the field gradient is 100 Oe cm -1 , and the geometry of the field lines is described in detail in [12,[START_REF] Sandre | Cebers in Giant Vesicles[END_REF]. Because magnetoliposomes have a higher magnetic susceptibility χ than the surrounding liquid, they move towards increasing field intensity. Their velocity is constant and corresponds to the balance of the magnetic force and the drag force exerted by the outer fluid. For ellipsoidal shapes the two forces have analytical expressions, which allows a precise measurement of χ.

= χ

. The magnetization at saturation of γFe 2 O 3 , m s =360 Oe, and the average volume of a nanoparticle V p are known from former studies [START_REF] Bacri | [END_REF].

Maghemite particles have diameters of the order of 10 nm. Their surface is coated by citrate ligands so that they bear negative charges at pH 7: this is essential for the colloidal stability of the ferrofluid. With no precaution, the concentration of tri-sodium citrate in ferrofluids is rather high (C s 0 ≈90 mM) because of the non adsorbed citrate species in solution in equilibrium with the citrate ligands. This can be decreased by dialysis through a cellulose membrane (Spectra/Por MWCO 12-14000, ROTH, France). With this method, the salt concentration C s 0 is reduced to 3 to 80 mM as obtained from the measured conductivity of the ferrofluid. In the same way as the nanoparticles, the electrolyte is also diluted in the vesicles compared to the initial ferrofluid. We make the assumption that the sodium citrate salt and the nanoparticles have identical entrapment yields:

0 0 p p s s C C C C =
. By contrast, if the sodium and citrate ions were not entrapped as efficiently as the particles (e.g. due to their faster diffusion before the bilayers close), then the electrolyte would be homogeneously diluted when adding pure water. For instance, starting with a ferrofluid of C s 0 =80 mM would lead to magnetoliposomes in a liquid bath with a salt concentration C s =0.8 mM. The suspension of magnetoliposomes would then have a moderate electrical conductivity (≈1mS cm -1 predicted). This is not observed experimentally, as the measured conductivity is less than 40 µS cm -1 . Therefore we conclude that the electrolyte is indeed co-encapsulated with the magnetic colloid in the liposome interiors.

Dimensions of magnetoliposomes are measured by optical bright field microscopy using image analysis software (Image, NIH). It is well known that the projected area of giant liposomes is smaller than the true surface of the membrane [15], the difference being absorbed by thermal undulations. The effect of a moderate tension is to flatten the fluctuations. Quasi-spherical vesicles can change their shape at constant volume by using the excess area hidden in the fluctuating membrane. In our case, tension originates from magnetic polarization forces acting on the magnetoliposomes. The shape of a liposome submitted to a magnetic field can be described by an ellipsoid with axial symmetry around the field direction. A slight deviation from the rotational symmetry is observed by confocal microscopy (Fig. 2). Giant liposomes filled with a ferrofluid have indeed a higher mass density than the suspending bath:

-3 3 4 cm g 10 5 10

- - × - ≈ ρ ∆
. Thus they are somehow flattened by gravity 0 g  in the z direction. In the next section we develop a simple theory that ignores gravity by considering only axisymmetric ellipsoids. In a standard approximation, shapes are predicted in the zero gravity limit, and gravitational corrections are considered afterward [8,16]. The dimensionless number b K R g g 4 0 0 ρ ∆ = has been introduced in [4] to predict deformations of vesicles due to gravity (K b is the membrane bending modulus). Magnetoliposomes have g values between 10 -3 and 1. This is why the largest ones are not perfectly axisymmetric. Nevertheless, we can clearly distinguish between prolate and oblate shapes, either elongated or flattened in the y direction (parallel to the applied magnetic field). We also point out here the discrepancy of our case with spontaneous transitions between a prolate and an oblate shape reported in [16]. In the latter case, the axis of rotational symmetry also flickers from parallel to perpendicular to the bottom surface of the observation chamber. In the shape transition triggered by a magnetic field, both prolate and oblate ellipsoids keep their axis parallel to In summary, the following experiments have been performed. First, we synthesize different ferrofluids with fixed particle concentration C p 0 , and varying ionic concentrations C s 0 . Then we encapsulate these ferrofluids in giant vesicles, thereafter characterized by optical microscopy. We measure their radius R 0 in the quasi-spherical state. Then we use a home built setup of magnetophoresis to simultaneously measure the ellipsoidal deformation under a magnetic field (value of e 2 and its sign), and the velocity under the field gradient (perpendicular to the field direction). As explained above, with these measurements we can calculate the concentrations C p and C s in individual liposomes.

Theory

The overall shape of a liposome subjected to a field H 0 must minimize the total free energy, which is the sum of the magnetic energy E m , the bending energy E b, and the surface energy of the vesicle ( )

τ π S R -4 0 2
, where τ is the membrane tension, and S the surface area of the vesicle. The idea is to calculate the free energy of an initially spherical liposome deformed at constant enclosed volume into an axisymmetric ellipsoid. In the expansion of the free energy as a function of the ellipticity, the terms of order higher than e 2 are neglected. The surface energy thus disappears since it varies as e 4 .

The magnetic term is proportional to the liposome volume. It is related to the shape through the where χ is the magnetic susceptibility of the encapsulated liquid. Using a standard expansion of D as a function of e [START_REF] Chikazumi | Physics of Magnetism[END_REF], one obtains the approximate expression:

( )

2 3 0 2 0 45 4 constant e R H E m πχ - ≈ (2)
This magnetic energy is decreased by positive values of e 2 , and thus it favors a prolate deformation.

In the following we show that the prolate/oblate transition under magnetic field is driven by a competition between the magnetic energy E m and the bending energy E b . The strong influence of the ionic strength on the deformation suggests that electrostatic interactions are involved between the charged nanoparticles and the slightly charged membrane. We make here the assumption that some charges exist on the DOPC phospholipidic bilayers. This assumption is supported by previous studies which report: i) the existence of a weak negative ζpotential of a few mV for vesicles made of natural egg lecithin [START_REF] Carrion | [END_REF] or synthetic phosphatidylcholine lipids [19]; and ii) the evidence of electrostatic repulsion between pure phosphatidylcholine bilayers, measured with a Surface Force Apparatus [20].

The behavior of charged membranes has been studied extensively in theoretical papers [21,22,23]. The electrostatic contributions to the bending modulus K b and to the spontaneous curvature c 0 have been considered in [24] in the Debye-Hückel approximation of the electrostatic potential and in the case of electrically uncoupled leaflets. It appears that electrostatic repulsion produces an enhancement of the bending modulus compared to the uncharged case: surface charges make the bilayer stiffer. The spontaneous curvature c 0 characterizes the bilayer asymmetry [25]. It has been derived in reference [24] for a bilayer bearing different surface charge densities on its two sides. In that case the salt concentrations were identical inside and outside the vesicle, and hence the Debye lengths as well. Here we develop another model that fits better to the case of magnetoliposomes:

asymmetry comes from two different screening lengths κ i -1 and κ o -1
, inside and outside respectively.

In a suspension of magnetoliposomes, the outer phase is a usual salt solution, while the inner phase is a complex colloidal system combining the salt and the charged magnetic nanoparticles. In order to obtain a simple solution of this problem, we describe this inner phase as a mixture of electrolytes: -the particles (concentration C p ), identified as magnetic polyions with an effective charge Z -the other ions (charge Z i for type i) having the same concentration C s i as in the suspending bath. In the outer phase, electrostatic interactions are screened by the ions (in practice sodium and citrate, but no particles) above a standard isotropic Debye length

2 1 2 2 1 4           = ∑ - i i i w o e Z C kT π e κ
, where e w is the dielectric constant of water (e w =78.5). Inside the liposomes, the decay length of electrostatic interactions is denoted κ i -1

and takes into account all the charged species, salt ions and nanoparticles. The modifications on the bending modulus and on the spontaneous curvature are estimated using the same approach as in [24]. At first we calculate the electrostatic energy E el of a simple body -an infinite cylinder with radius R -then we expand it as a function of R. Finally we identify the terms with the standard expression of the bending energy, supposing that these expressions remain the same for a body of general geometry:

2 0 2 3 2 3 2 2 2 1 2 1 8 3 8 3 2 1 2 1 1 1 4 2 1 2       - ≡         + + - + + ≈ c R K R R R R R E el b o i o i i o w el κ κ κ κ κ κ e πσ π (3)
The elastic constants are derived in the assumption of a slight difference between κ i and κ o . The bending modulus is then Here σ denotes the surface charge density. It is inferred that σ is the same on both monolayers, because the pH is buffered inside and outside by the citrate salt. The asymmetry between the two sides of the membrane arises from different electrolyte compositions inside and outside the vesicle. This small difference has a negligible effect on the modulus, but it is dominant for the spontaneous curvature.

In the model of a liposome filled with a ferrofluid, we make an additional approximation. We suppose that the difference between κ i and κ o is due only to the presence of nanoparticles inside, and their absence outside. In other words, we neglect the difference of salt concentrations inside and outside, responsible for a classical isotropic c 0 . The model focuses on the field dependence of the inner Debye length κ i -1

. The concentration profile of nanoparticles near the inner leaflet of the membrane is indeed sensible to  H . Thus, a magnetic term has to be added in the linearized Poisson-Boltzmann equation. Due to that the calculation of the electrostatic energy of the membrane in the presence of the magnetic field turns out to be rather complicated problem since simultaneous solution of the electrostatic and magnetic field equations is involved. As result at the curvature expansion of the electrostatic energy of the membrane its anisotropy arising due to the magnetic interactions appears. This electrostatic energy curvature expansion is carried out in the forthcoming publication [27]. Its results confirm the conclusions of much simpler model described below. The simple model is based on the assumption of the locally flat membrane and the reduction of the solution of the coupled electrostatic and magnetostatic field equations to 1D approximation. Obtained in such way screening constants are substituted in the expressions of the bending elasticity constant and spontaneous curvature leading to the anisotropic curvature elasticity energy of the membrane. Further consideration is carried out in the frame of the simple model. We begin by solving the magnetostatic problem in a flat membrane geometry, with a magnetic field at an angle α with the normal to the plane of the membrane. Let x be the coordinate in the direction normal to the membrane.

Because the vesicle radius is three orders of magnitude larger than the particles sizes, we can neglect the membrane curvature and the finite size of the nanoparticles. Thus we denote the nanoparticle concentration C p (x) at a distance x from the membrane, the x component of the field H x (x), and the electric potential ψ(x). Far away from the membrane those variables tend respectively towards C p , H 0 cos(α) and zero. The solution of Maxwell

equation 0 4 = + ) M H ( div   π is: p p p p x ) x ( C C ) cos( H ) x ( H χ π χ π α 4 1 4 1 0 + + = (5)
The net flux of nanoparticles in the x direction is then calculated using their mobility u under an osmotic, an electrophoretic and a magnetophoretic force:

x H H ) x ( uC x ) x ( C ukT x Z ) x ( uC J x x p p p p x ∂ ∂ χ ∂ ∂ ∂ ∂ψ + - - = (6)
At equilibrium, J x =0 and: compares the two pressures acting on the membrane: the magnetic pressure calculated at the liposome poles ( 2 2 0 2 πχ H ) and the osmotic pressure exerted by the nanoparticles (C p kT). Numerical values of λ range from 10 -3 to 5×10 -2 , as seen in Table 1. Looking back to the electrostatic potential ψ(x), an anisotropic decay factor can be defined, including the screening by all the charged species:

( ) ( ) x Z ) x ( C x ) x ( C ) x ( C kT ) ( cos H C ) x ( C kT p p p p p p p p ∂ ∂ψ ∂ ∂ χ π α χ π χ π - =         + + +
2 1 2 1 1       + + = ) ( cos d ) ( o i α λ κ α κ (8) It depends through d Z C Z C p i s i i = ∑ 2 2
on the ionic strength of the dispersion and through λ on its magnetic characteristics. The anisotropic part of κ i is negligible as long as the screening is dominated by the salt ions, i.e.

for high ionic strengths (small d). In the opposite case, the angular dependence of the distribution of particles induces an angular dependence of

) ( i α κ
, and hence of c 0 (α ) . The largest screening length turns out to be on the poles. The reason for that is rather simple. Since the magnetic pressure acting near the poles pushes the magnetic polyions to the membrane it augments near it the concentration of the charges having the same sign as on the membrane and thus diminishes the ecranization of charges on it. An increase of the screening length near the magnetic poles (α=0 and α=π) decreases locally the spontaneous curvature: qualitatively, the oblate shape is favored because it has a lower curvature at the poles than at the equator. Equations ( 4) and ( 8) enable one to derive the spontaneous curvature at any angle α:

( )       + - + ≈ ) ( cos d d d ) ( c o α λ κ α 2 0 1 1 1 3 (9)
We can now calculate precisely the bending energy of the membrane,

( ) ∫       - + + = dS ) ( c R R K K E el b b b 2 0 2 1 1 1 2 1 α
. It involves the two local principal radii of curvature R 1 and R 2 of the bilayer [25]. Due to the dependence of the bending elasticity modulus on the screening constant and its anisotropy described by the relation (8) anisotropic under the action of the magnetic field, the bending elasticity modulus of the membrane is also anisotropic. Nevertheless its account is not necessary at the calculation of the curvature elasticity energy of the vesicle up to the first order terms in 2 e since the direct calculation shows that its contribution in that approximation disappears exactly. That is why for the bending elasticity modulus the is assumed, what does not reflect the dependence of the modulus on the magnetic field strength and its direction. As before, we look at an ellipsoidal deformation of a quasi-spherical vesicle.

Integration leads to:

( ) ( ) ( ) 2 2 0 3 2 1 405 32 e R d d K K E o el b b b κ λ π - + + ≈ (10)
Exact treatment of the electrostatic energy of the vesicle under the action of the magnetic field carried out in [27] shows that besides the anisotropies of the elastic modulus discussed above, there is specific 3D efect leading to the anisotropy of the bending modulus depending on the tangential component of the magnetic field strength.

Howewer its contribution to the electrostatic energy is negligible since it does not contain -in comparison with the expression (10) -a multiplier 2 0 0 ) ( R χ which, as it follows from data in Table 1, is very large. E b enters in the expansion of the total free energy with the same order in e 2 as the magnetic energy E m , but with the opposite sign. Indeed it is possible to explain a prolate deformation or an oblate one depending whether E m or E b is predominant.

Results and discussion

The experimental measurements of R 0 and χ on given liposomes are computed in Eq. (1) to get the slopes in the expansion of E m as a function of e 2 . Then we convert the values of χ into particle concentrations C p , inferring a mean value for the nanoparticles diameters D=12×10 -7 cm, or equivalently a magnetic susceptibility per particle χ p =8.6×10 -19 cm 3 . The salt concentration is deduced by C s =C s 0 C p /C p 0 . From stoichiometry it follows that Z 1 =3

and C 1 =C s for the citrate anions, Z 2 =1 and C 2 =3C s for the sodium cations. Then we calculate the electrostatic decay factor κ o caused by the salt ions only. The two parameters λ and d appearing in the anisotropic term of ) ( i α κ are computed from the values of χ, C p , C s and Z=25, a reasonable estimate of the net charge per particle (including the adsorbed counterions). The surface charge density of the bilayer is inferred to be σ=800 ues cm -2 .

It represents about 1% of negatively charged lipids. Finally we compute E b /e 2 from Eq. (10). The ratio of the two competing terms E b /E m is plotted vs. C s on Fig. 3. For the higher ionic strengths where the screening of charges is the strongest, the magnetic term dominates the electrostatic contribution to the bending energy: the deformation is indeed prolate. The field dependent elongation has been analyzed in a former work in order to determine the bending modulus of the DOPC bilayer, K b =(8.6±0.8)×10 -13 dyn cm [26]. But a dilution of the ions in the bulk by two orders of magnitude is sufficient to reverse the ratio of E b and E m , thus explaining the oblate shape. The critical value for this ratio is 1, and it clearly separates the graph into domains of prolates and oblates. The scatter of Fig. 3 is due to two parameters, other than C s , that vary in a polydisperse suspension of magnetoliposomes: the vesicle radius R 0 and the encapsulated particle concentration C p (or equivalently the λ-parameter). For example the magnetic energy E m varies like 10) because of the polydispersities of R 0 and λ. This gives the residual scatter in the data. However the sensitivity on C s is fully established: C s is varied over two orders of magnitude, and this inverts the weights of E b and E m in the total energy. The set of values for the guessed parameters -Z, χ p (or equivalently D) and σ -is reasonable. The point that needs to be further discussed in our model is the approximation of a flat interface to derive the magnetic field H x (x) as given by Eq. ( 5). Complete treatment of the problem based on the curvature expansion of the electrostatic energy of the membrane of the arbitrary shape is carried out in [27]. Results obtained completely confirms qualitative picture considered above. It should be remarked that results obtained for the anisotropy of the spontaneous curvature in the frame of the complete treatment shows the coincidence with qualitative simple model even up to the numerical coefficient. Therefore current model is sufficient to point out that the electrostatic double layer of membrane is deformed by the magnetic field via its effect on the charged magnetic nanoparticles.

We have reported an original shape transition of giant magnetic liposomes triggered by the ionic strength.

The transition between prolate and oblate ellipsoids has already been observed under an applied ac electric field

[28]: the vesicles are prolate in the 'conductive regime' (low frequency) and oblate in the 'dielectric regime' (high frequency). It has also been predicted for an applied magnetic field: Helfrich showed that the sign of magnetic anisotropy of the lipid molecules themselves leads either to prolate or oblate deformations [25,29].

However the magnetically induced ellipticity has been far too small to be detected experimentally, even in a strong field H 0 =10 4 Oe (1 Tesla). Here we have taken advantage of the magnetic properties of the encapsulated solution to improve the coupling between the magnetic field and the bending energy of the vesicles. The main result of our theoretical description is an anisotropic distribution of particles at the neighborhood of the membrane when a field is applied. Both the Debye screening length at the inner monolayer and the spontaneous curvature of the bilayer maintain this anisotropy. This model corresponds to the regime of electrostatic interaction between a weakly charged phospholipid bilayer and ionic encapsulated species. This might be an important point for the description of biological membranes which are surrounded by charged mesoscopic objects and embedded in salty media. Gravity is along the z axis. The hot spots are lipid-particle aggregates that either settle at the bottom of the cell or make chains attracted by the strongly magnetized giant liposome.

  To get the particle concentration C p , we use the ferrofluid property of superparamagnetism: the susceptibility is the product of C p and the susceptibility per particle



  H . Let a be the semi- axis parallel to the magnetic field  H and b the value of the two other semi-axes perpendicular to  H . The two types of deformations and their amplitudes are characterized by a single parameter of ellipticity: positive for an elongated liposome (prolate ellipsoid) and negative for a compressed one (oblate ellipsoid).

  demagnetizing factor D(e): O. Sandre, C. Ménager, J. Prost, V. Cabuil, J.-C. Bacri, and A. Cebers, Phys. Rev. E 62, 3865 (2000) Author manuscript DOI: 10.1103/PhysRevE.62.

  -Boltzmann theory and the Debye Hückel approximation, a solution is found for a small variation of nanoparticle concentration around the bulk value ( δC x C x ). Then the linear dependence of δC p (x) with the electrostatic potential ψ(x) becomes anisotropic,

3 0 R.

 30 and the anisotropic bending energy E b like 2 0 R . The scatter of the data has been reduced by plotting on Fig. The error bars estimate the uncertainty in our measurements (shape analysis, magnetophoresis and conductimetry). It is impossible to fit all the data with a single simulation of Eqs. (1) and (

Fig. 1 :

 1 Fig. 1: The two competing shapes of magnetoliposomes subjected to a magnetic field (amplitude H 0 =400 Oe, field direction along the arrow); (a) prolate ellipsoid for a high ionic strength ferrofluid (C s 0 =85 mM); (b) oblate

Fig. 2 :

 2 Fig. 2: Cross-sections of an oblate magnetoliposome in the three principal planes xy, xz and yz. A stack of 150 pictures in the sample plane (xy) is taken every 0.5µm along the z axis with a confocal microscope (TCS4D, Leica). The lipid bilayer is labelled with the fluorescent dye Di-6-ASPBS (N-(4-sulfobutyl)-4-(4-(dihexylamino)styryl)pyridinium). Excitation is provided by the 488nm line of an Argon laser, and fluorescence is collected through a fluorescein filter set. The ellipsoid semi-axes are a=19µm and b=32.5µm. The ferrofluid has a low ionic strength (C s 0 =5 mM). The white arrow indicates the magnetic field (intensity H 0 =200 Oe).

Fig. 3 :

 3 Fig. 3: Ratio of the two competing terms E b /E m in the expansion of the free energy as a function of ellipticity. They are calculated for oblate (open markers) and for prolate (filled markers) liposomes using both experimental data (C p , C s , R 0 ) and assumed parameters (see text). Solid line is plotted from Eqs. (1) and (10) for typical values R 0 =15 µm and λ =5×10 -3 .

Table 1 :

 1 Examples of numerical data for oblates (clear cells) and prolates (shadowed cells). R 0 and e 2 are measured by optical microscopy. The parameters d and λ are calculated from experimental measurements of the concentrations C s and C p . The last three columns contain important features of the model: the screening length κ o -

	25.2	-0.397	0.96	3.0×10 -3	45	8.4	280
	13.3	-0.596	0.62	4.8×10 -3	29	2.2	230
	21.4	-1.17	0.62	3.2×10 -3	36	4.0	280
	24.6	-0.845	0.19	3.1×10 -3	20	0.73	380
	22.5	0.496	0.19	5.3×10 -3	15	0.32	290
	50.4	0.343	0.19	2.8×10 -3	21	0.85	400
	21.1	-1.31	0.14	1.5×10 -3	24	1.3	600
	31.6	0.466	0.11	2.3×10 -3	17	0.45	540
	14.1	0.536	0.074	4.5×10 -2	3.2	3.1×10 -3	150
	13.0	0.586	0.056	6.1×10 -3	7.8	4.2×10 -2	440

due to the salt ions only, the ratio of the electrostatic contribution to the bending modulus of membrane, and the local value of the spontaneous curvature radius at the magnetic poles (α=0 and α=π).

Acknowledgments

We thank J. Mertz and E. Karatekin for reading the manuscript and M. Blanchard-Desce for the gift of Di-6-ASPBS. O. Sandre, C. Ménager, J. Prost, V. Cabuil, J.-C. Bacri, and A. Cebers, Phys. Rev. E 62, 3865 (2000) Author manuscript DOI: 10.1103/PhysRevE.62.3865


Baptiste Bouillot
email: bouillot@emse.fr

Son Ho-Van

Jérome Douzet

Ana Cameirao

Daniel Garcia

Jean-Michel Herri

Cameirao Ana

Garcia Daniel

Thermodynamics of cyclopentane hydrates formation from brine: experimental and modelling study

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Hydrates can be used to remove salt from saline water or to clean water as a purification method

Objective

  NaCl, KCl, NaCl-KCl, CaCl 2 , Na 2 SO 4 , MgCl 2 , MgCl 2 -NaCl, MgCl 2 -NaCl-KCl at different concentrations  Develop four models based on Standard Freezing Point Depression (SFPD) equation, Hu-Lee-Sum (HLS) correlation, Kihara and Activity-Based Occupancy Correlation (ABOC) to predict the equilibrium conditions. = 𝑚 × 𝑎 𝑤 2 + 𝑛 × 𝑎 𝑤 + 𝑝 occupancy of cavities correlated to water activity Ho-Van, S. et al., AIChE Journal 64(6) (2018) Ho-Van, S. et al., Industrial & Engineering Chemistry Research 57(43) (2018)


Sébastien Ménigot
email: sebastien.menigot@eseo.fr

Nesrine Houhat
email: n.houhat@crti.dz

Jean-Marc Girault
email: jean-marc.girault@eseo.fr

Optimal Prefiltered Stochastic Transmitted Waves for Fat Inclusion Detection in Milk with Harmonic Ultrasound

Keywords: Nonlinearity, optimal wave, genetic algorithm

In non destructive testing, improvements have been made possible by taking into account the harmonic frequencies, as in agri-food domain. The transmitted signal are often selected empirically as a fixed-frequency Gaussian pulse, by taking into account the transducer bandwidth only. However, waveform should take into account all the features of the ultrasound system and of the medium. To design the waveform, a genetic algorithm looks for the best stochastic wave. However, one of limitations in this optimization process is the high amount of transmitted waves. To reduce this number, instead of transmitting wide band stochastic waves, narrow band stochastic waves limited by the transducer bandwidth are preferred. The optimization was thus applied on the detection of fat cluster in milk by maximizing the signal-to-noise ratio (SNR), while decreasing the amount of transmitted waves. Twelve combinations from different limited bandwidths of transmitted waves were tested. Whereas the low cut-off frequencies did not change the performances, the high cut-off frequencies affected the convergence speed. In this study, it is shown that the best optimization was twelve times faster with the high cut-off frequency of 5.6 MHz and led to a gain of 62% compared to the SNR obtained with a best fixed-frequency sine wave.

I. INTRODUCTION

Nonlinear ultrasound systems have improved the structure/inclusion detection in fluids by enhancing the level of backscattered harmonics. This is particularly the case for detecting fat clusters in the agri-food domain [START_REF] Feng | Ultrasound technologies for food and bioprocessing[END_REF], i.e. in cheese manufacturing. Its simple principle consists in transmitting an ultrasound sinus wave of frequency f 0 into the medium being explored and in receiving harmonic components (2f 0 , 3f 0 ,...) generated during the wave propagation. Several methods of harmonic extraction have been developed by filtering the echoes or by coding successive transmitted waves [START_REF] Burns | Nonlinear imaging[END_REF]. The second kind of methods have improved the signal-to-noise ratio (SNR) while ensuring a good axial resolution, such as pulse inversion imaging [START_REF] Simpson | Pulse inversion doppler: A new method for detecting nonlinear echoes from microbubble contrast agents[END_REF]. Since the most commonly used is the pulse inversion imaging, we only focused our study with this technique.

Even if most of improvements includes post-processing to extract harmonics, the choice of the transmitted wave is decisive. However, the usual settings of the transmitted wave are empirically constrained to a fixed-frequency sine wave, because the optimal settings requires knowledges experimentally inaccessible about the transducer and the fluid. A first solution carries in a suboptimal solution by finding the parameter of the waveform, such as frequency [START_REF] Ménigot | Optimization of contrast to tissue ratio by frequency adaptation in pulse inversion imaging[END_REF]. However, this optimization problem can be solved in an optimal way by sending stochastic waves iteratively, without any assumption on the waveform [START_REF] Ménigot | Optimization of contrast resolution by genetic algorithm in ultrasound tissue harmonic imaging[END_REF].

In this study in the nonlinear imaging context using pulse inversion, a genetic algorithm looks for the best stochastic wave by optimizing the SNR between a reference (milk without any fat cluster) and fat cluster as elements in the milk. One of limitations in this optimization process is the high amount of transmitted waves required to reach the optimal SNR enhancing the fat cluster. To reduce this amount, instead of transmitting wide band stochastic waves by the transducer, stochastic waves were prefiltered by a narrow band filter similar to the transducer bandwidth. Here, we propose to break the principle of not taking into account any a priori information by proposing a prefiltering based on the spectral parameters of the transducer in an acceleration of the optimization process, while reaching the global optimum.

II. OPTIMIZATION WITH PREFILTERED STOCHASTIC WAVES

Finding the best prefiltered stochastic waveform w filtered is a hard optimization problem, which is more complex with the nonlinear imaging context using pulse inversion. To solve it, a pulse inversion system is modified as an iterative close loop system including a metaheuristics in the feedback. Among the metaheuristics, genetic algorithm [START_REF] Haupt | Practical Genetic Algorithms[END_REF] can optimize a cost function by selecting parameters concatenated in a vector. Therefore, it seems to be adjusted to our problem, since prefiltered stochastic waveforms can be written as a vector containing N s samples. The problem can be simply written:

w filtered = arg max w filtered (SNR), (1)
where cost-function allows us to enhance the backscattered power of fat cluster P fat cluster in comparison with the backscattered power of the surrounding milk P milk , such as: SNR = 10 log 10 P fat cluster

P milk . (2)
First, following the procedure of the genetic algorithm [START_REF] Haupt | Practical Genetic Algorithms[END_REF], 12 prefiltered stochastic waveforms w filtered are randomly generated from white noise (magnitude between -1 and 1) by the genetic algorithm. They are prefiltered by a finite impulse response (FIR) filter where their bandwidths at -3 dB are between the low cut-off frequency f low and the high cut-off frequency f high [START_REF] Smith | Spectral audio signal processing[END_REF]. Before the transmission, the power of the prefiltered stochastic waveforms w filtered are normalized by adjusting the amplitude A, so that They are constant to the power P x ref of the impulse response of the transducer driving pressure A 0 :

x q = (-1) q • A • w filtered , (3)
with q = 1, 2 and

A = A 2 0 • P x ref w T filtered • w filtered , (4)
where w T denotes the transpose of the matrix w. Note that the number of samples N s of each waveform is set so that the duration T of the stochastic waveform w corresponded to 100% of the fractional bandwidth of the transducer. Finally, the transmitted and prefiltered stochastic waves x 1 and the wave with their opposite phases x 2 (x 1 = -x 2 ) are sequentially transmitted into the milk including fat cluster. Second, for each echoes of the 12 waveform couples, even harmonic components are respectively extracted by the sum z of their two respective echoes y 1 and y 2 . The 12 respective SNR are assessed on that RF line z.

Finally, a genetic algorithm generates a new prefiltered stochastic waveform w to optimize the SNR. Among the 12 prefiltered stochastic waveforms, only the 6 best individual waveforms which maximized the SNR are selected. To construct them, the crossover operator mixed the best previous waveform with one of the 5 remaining better waveforms. Finally, 40% samples were mutated to obtain robust optimization. Moreover, the 12 new stochastic waveform are prefiltered by the FIR filter to guarantee their bandwidth.

III. SIMULATION MODEL

As described previously, the prefiltered stochastic waves x and their opposite phases were generated digitally with Matlab (Mathworks, Natick, MA, USA). They were transmitted by a transducer centred at 4 MHz with a bandwidth of 80% at -3 dB

Fat Cluster

Fig. 1. Grid of properties of the medium constituted of milk (medium 1) with a fat cluster (medium 2): c is the speed of sound wave velocity, ρ the density and B/A the nonlinearity parameter. The ultrasound transducer was at a depth of 0 mm, here at the top. As an illustration, a harmonic RFline z is depicted in green.

(2.4 -5.6 MHz) and 200% at -10 dB (0.25 -14 MHz). The pulse waves generated were sequentially propagated nonlinearly into the attenuating medium of milk [START_REF] Anderson | A 2d nonlinear wave propagation solver written in open-source matlab code[END_REF]. In the milk, a 4 mm-diameter fast cluster was inserted into the milk at 12 mm below the surface (Fig. 1). Note that the medium properties are shown in table I. Finally, the echoes were filtered by the same transducer than this one in transmission. The genetic algorithm proposes new prefiltered stochastic waves and the procedure is repeated until the end of the optimization.

IV. RESULTS

The optimization process was tested with 12 combinations from different limited bandwidths of transmitted waves (4 low and 3 high cut-off frequencies). The low and the high cutoff frequencies are listed in table II. Note that these values were chosen from the cut-off frequencies of the transducer. Each setting of bandwidth is repeated 10 times, in order to remove the random effect of the algorithm and to check if the optimization converges to the same value. Note that the other settings were always the same:

• the driving pressure A 0 was set to 400 kPa; • the waveform was described by N s = 40 samples, since its duration had to represent 100% of the fractional bandwidth of the transducer. Moreover, as an illustration, the results were compared to a standard fixed-frequency Gaussian pulse for which its frequency was set to optimize the SNR. First, the low cut-off frequency of the prefiltered waveform was fixed at the same value than the low cut-off frequency of the transducer, i.e. 2.4 MHz. The high cut-off frequency of the prefiltered waveform was changed to the value described in table II. Note that these values was chosen to have at least the same bandwidth of transducer at -3 dB and thus to keep the spatial resolution. Fig. 2 show the optimization of the SNR as a function of the genetic algorithm iteration for the 3 high cut-off frequencies of the prefiltered waveform. Whatever the high cut-off frequency, the SNR increased during the genetic algorithm iterations. At the end of the optimization process, the SNR was increased by 3 dB in comparison with the SNR obtained with a best fixed-frequency Gaussian pulse. Moreover, the maximal SNR were higher by using prefiltering waveform than the optimization without prefiltering. However, the narrower bandwidth allowed us to reach better SNR. Finally, with the better prefilter settings (a high cut-off frequency of 5.6 MHz), the convergence was quicker, since a SNR of 8.25 dB can be reach with only 4 iterations, compared to the 12 iterations required without prefiltering.

Second, the high cut-off frequency of the prefiltered waveform was fixed at the same value than the high cut-off frequency of the transducer, i.e. 5.6 MHz. The low cut-off frequency of the prefiltered waveform was changed to the value described in table II. Fig. 3 show the optimization of the SNR as a function of the genetic algorithm iteration for the 4 low cut-off frequencies of the prefiltered waveform. Whatever the low cut-off frequency, the SNR increased during the genetic algorithm iterations. At the end of the optimization process, the SNR was increased by 3.4 dB in comparison with the SNR obtained with the best standard fixed-frequency Gaussian pulse. Moreover, the maximal SNR were always higher by using prefiltering waveform than the optimization without prefiltering and than the waveform prefiltered with the high cut-off frequencies described in Fig. 2. However, whatever the low cut-off frequency, the SNR is equivalent. Nevertheless, a lower cut-off frequency slightly smaller than the cut-off frequency of the transducer makes it possible Fig. 2. Median of optimized SNR from 10 simulations as a function of the genetic algorithm iteration for different high cut-off frequencies of the prefiltered waveform. Note that the SNR value shown here are the median for 10 simulations in the same conditions. The low cut-off frequency of the prefiltered waveform is fixed at the same value than the low cut-off frequency of the transducer, i.e. 2.4 MHz. As an illustration, the optimization without any constraint on the bandwidth is depicted in black. Moreover, the performance using the best fixed-frequency Gaussian pulse is depicted in dash line. to statistically obtain a slightly higher SNR. Note that this improvement was not significant for the 40 simulations (4 settings repeated 10 times).

To summarize, table III shows the iterations to reach a SNR of 8.5 dB and the maximal SNR at the end of the optimization process. As shown in Figs. 2 and2, the optimization process with prefiltering can be about 12 times faster (187/15). After only one iteration, the SNR obtained by prefiltering stochastic waveform was higher than this one obtained by stochastic waveform without prefiltering and by the best fixed-frequency Gaussian pulse. Moreover, after only 6 iterations, the SNR was increased by 3 dB compared to the SNR obtained with the best fixed-frequency Gaussian pulse. Finally, at the end of the optimization process, the maximal SNR was higher with prefiltered stochastic waveform than without prefiltering. 3. Median of optimized SNR from 10 simulations as a function of the genetic algorithm iteration for different low cut-off frequencies of the prefiltered waveform. Note that the SNR value shown here are the median for 10 simulations in the same conditions. The high cut-off frequency of the prefiltered waveform is fixed at the same value than the high cut-off frequency of the transducer, i.e. MHz. As an illustration, the optimization without any constraint on the bandwidth is depicted in black. Moreover, the performance using the best fixed-frequency Gaussian pulse is depicted in dash line.

However, note that the genetic algorithm should be reach the maximal SNR without prefiltering, but it would require a higher number of iterations.

As an illustration, Fig. 4a shows the optimal transmitted stochastic wave x obtained by the genetic algorithm at the transducer output when the prefiltered stochastic waveform w is optimal. Note that this signal was propagated in the milk. The harmonic echoes z obtained by using the optimal transmitted stochastic wave x is shown in Fig. 4b. As expected, it was higher to the harmonic echoes obtained by using the better fixed-frequency Gaussian pulse (shown in red).

V. DISCUSSION AND CONCLUSIONS Prefiltered stochastic waveform were automatically selected in order to optimize the SNR assessed in a pulse inversion imaging system to detect fat cluster in milk. To help the optimization process, the stochastic waveform were prefiltered to constrain the research space. Thus the SNR could be increased by 5% compared to the optimization without prefiltering and 62% compared to a fixed-frequency Gaussian pulse.

However, the settings of the cut-off frequencies had an influence on the better SNR obtained after optimization. First, the high cut-off frequency had to be close to the cut-off frequency of the transducer. If the high cut-off frequency were higher than the high cut-off frequency of the transducer, the harmonic generated during the propagation could be detected, because out of bandwidth. Second, the low cut-off frequency had no effect on the better SNR. Therefore, the best prefilter for stochastic waveform is a low-pass prefilter where its cut-off frequency is the high cut-off frequency. Finally, as the optimization is done in a small area by taking the transducer features into account, the genetic algorithm looks twelve times faster for the optimal waveform. Moreover, the optimization including a prefiltering makes possible to obtain higher SNR than an optimization without prefiltering, in just one iteration.

In conclusion, this method improved the SNR without strong assumptions on the waveform. This settings should be suitable to find more easily the optimal waveform in non destructive testing application.

  Fig.3. Median of optimized SNR from 10 simulations as a function of the genetic algorithm iteration for different low cut-off frequencies of the prefiltered waveform. Note that the SNR value shown here are the median for 10 simulations in the same conditions. The high cut-off frequency of the prefiltered waveform is fixed at the same value than the high cut-off frequency of the transducer, i.e. MHz. As an illustration, the optimization without any constraint on the bandwidth is depicted in black. Moreover, the performance using the best fixed-frequency Gaussian pulse is depicted in dash line.

Fig. 4 .

 4 Fig.4. (a) Optimal transmitted stochastic wave x obtained by the genetic algorithm at the transducer output when the prefiltered stochastic waveform w is optimal. (b) Harmonic backscattered pressure obtained by pulse inversion and optimal prefiltered waveform. As an illustration, the harmonic RF line in red is obtained by using the better fixed-frequency Gaussian pulse.

TABLE I MECHANICAL

 I PROPERTIES OF EACH EXPLORED. N IS THE SYMBOL OF THE GAUSSIAN DISTRIBUTION. N (1531 m/s, 0.008 m 2 /s 2 ) N (1060 m/s,75 m 2 /s 2 ) Density N (1025 kg/m 3 ,0.008 kg 2 /m 6 ) N (928 kg/m 3 , 75 kg 2 /m 6 )

				Milk [9]	Fat Cluster [10]
		Speed of sound c
		Nonlinearity parameter B/A	5.1	10.3
			TABLE II
	LIST OF LOW AND 4 HIGH CUT-OFF FREQUENCIES OF THE PREFILTERED
	STOCHASTIC WAVEFORM. NOTE THAT THE CUT-OFF FREQUENCIES OF THE
	TRANSDUCER ARE WRITTEN IN BOLD
	Low cut-off freq.	0.44 MHz 0.88 MHz 1.17 MHz 2.4 MHz
	High cut-off freq. 5.6 MHz	11 MHz	17 MHz

TABLE III NUMBER

 III ITERATIONS REQUIRED TO REACH 8.5 DB WITH THE GENETIC ALGORITHM AND THE MAXIMAL SNR AT THE END OF THE OPTIMIZATION PROCESS.

	Low cut-off High cut-off	Iterations	Maximal
	frequency	frequency	to reach 8.5 dB	SNR
	0.88 MHz	5.6 MHz	15	8.9 dB
	None	None	187	8.6 dB


Nesrine Houhat
email: n.houhat@crti.dz

Vincent Tournat
email: vincent.tournat@univ-lemans.fr

Sébastien Ménigot
email: sebastien.menigot@eseo.fr

Tarek Boutkedjirt
email: tboutkedjirt@usthb.dz

Redouane Drai
email: r.drai@crti.dz

Jean-Marc Girault
email: jean-marc.girault@eseo.fr

Optimization of Bifrequential Pump Excitation Parameters by Genetic Algorithm in Nonlinear Vibro Acoustic Wave Modulation Technique

Abstract-The nonlinear Vibro-Acoustic wave Modulation (VAM) is a widely used technique in nonlinear nondestructive testing. Many research, in this area, are devoted to find the best combination of the excitation parameters, especially the right frequencies to the optimization of the damage detection sensitivity. In a previous study, we have proposed a novel method using the optimal command principle for a VAM system. Indeed, this method has permitted to find automatically the best pump frequency maximizing the nonlinear modulation (NM) effects in a multiple scattering sample. The cost function to be optimized thanks to the pump wave frequency is the correlation coefficient (CC) between a reference output signal without pump and an output modulated probe signal with the presence of the pump wave. In the present work, we aim to improve the nonlinear damage detection by exciting simultaneously two resonance modes of the medium. We consider a new excitation waveform for the pump wave, composed by the sum of 2 frequencies, instead of a monochromatic sine wave. We propose, then, to find automatically the frequencies fp1 and fp2 giving the optimal sensitivity for damage detection. The tested sample is an aluminium bar where nonlinear scatterers can be controllably added or removed. For both fp1 and fp2, the genetic algorithm (GA) seems to stuck on one of the resonance frequencies until it converges to the optimal combination. The corresponding values of the CC are quite similar to those obtained by a single pump frequency.

I. INTRODUCTION

The nonlinear Vibro-Acoustic wave Modulation (VAM) is a widely used technique in nonlinear nondestructive testing of materials where the conventional methods reach their limits in the detection of small damage such as a micro-crack. In the VAM tests, the sample is excited simultaneously by a low frequency pump and a higher frequency probing waves [START_REF] Donskoy | Vibro-acoustic modulation nondestructive evaluation technique[END_REF], [START_REF] Pieczonka | Nonlinear vibroacoustic wave modulations for structural damage detection: an overview[END_REF]. In the presence of damage, the two frequencies are mixed giving rise to new frequency components often referred to as sidebands. Generally, the pump frequency is chosen to correspond to one of the vibration modes of the studied sample in order to amplify the vibration response [START_REF] Pieczonka | Nonlinear vibroacoustic wave modulations for structural damage detection: an overview[END_REF]- [START_REF] Liu | Influence of low-frequency parameter changes on nonlinear vibro-acoustic wave modulations used for crack detection[END_REF]. A preliminary modal analysis is, then, necessary for identifying and selecting the resonance frequencies of the sample. It is important to underline that such a procedure requires a further experimental setup and time. Furthermore, a possible drawback of a modal excitation is that the crack or the damaged zone can be located at a strain node of the pump or probe wave, which compromises the generation of the nonlinear modulation.

Many efforts in this area are devoted to find the best combination of the excitation parameters, especially the right frequencies, for the optimization of the VAM sensitivity [START_REF] Liu | Influence of low-frequency parameter changes on nonlinear vibro-acoustic wave modulations used for crack detection[END_REF], [START_REF] Dziedziech | Efficient swept sine chirp excitation in the non-linear vibro-acoustic wave modulation technique used for damage detection[END_REF]. In a previous study [START_REF] Houhat | Optimal pump excitation frequency for improvement of damage detection by nonlinear vibro acoustic modulation method in a multiple scattering sample[END_REF], we have presented a combination of the optimal command [START_REF] Ménigot | Optimization of contrast-to-tissue ratio by frequency adaptation in pulse inversion imaging[END_REF]- [START_REF] Houhat | Optimal stochastic excitation for linear flaw detection in a solid material[END_REF] and the VAM technique. A closed loop VAM system able to find automatically the optimal pump frequency that maximizes the nonlinear modulation effects has been presented. The correlation coefficient between a reference output probe signal without the pumping wave and an output modulated probe signal with a pumping wave was considered as the cost function. The main outcome was the existence of optimal pump frequency, corresponding to a one resonance frequency among several resonance frequencies of the sample. This study has permitted to obtain a suboptimal solution by fixing the pump excitation waveform.

The main idea of the present study, is to attempt to improve the previous results by exciting the sample at more of one of his optimal resonance frequencies. Since the medium present several resonance frequencies, we suppose that it should amplify the nonlinear modulation effects between the probe and the pump waves. We propose to test the simplest case; those of exciting simultaneously two resonance frequencies of the sample. A novel suboptimal pump excitation waveform is, then, considered. It is a sinusoidal bifrequential signal composed from 2 frequencies f p1 and f p2 . An experimental modal analysis is performed to identify the resonance modes of the sample. The two continuous sine pump waves are generated simultaneously with the probe wave. Experimental configurations including different positions and amount of localized damage have been tested.

II. VAM SYSTEM FOR BIFREQUENTIAL PUMP

OPTIMIZATION

The proposed method is an optimal command method using a closed loop, in order to optimize the VAM sensitivity of crack detection. It makes it possible to find the pump frequency which optimizes the nonlinear modulation effects. The cost function and its parameters need to be adequately chosen. The correlation coefficient ρ between the received probe coda signal without pumping and with the pump excitation constitutes our cost function. With the presence of a nonlinear scatterer in the sample, the probe signal is expected to be modulated by the pump excitation which induces to a decorrelation between the two signals. A correlation coefficient equal to unity means that there is no influence of the pump on the probe wave, i.e., no nonlinear modulation effect in the sample. On the contrary, for the same pump amplitude value, a deviation of ρ from 1 indicates the presence of nonlinear damage leading to a nonlinear modulation effect. In a theoretical point of view, the problem consists in calculating:

[f * p1 , f * p2 ] = arg min fp1,fp2 (ρ(f p1 , f p2 )). (1)
During the optimization process, for each iteration k, the same probe excitation signal x s (t) is transmitted to the sample without the pumping signal, and the received reference probe coda signal y r,k (t) is recorded for a next use. In a second time, a pump excitation x p,k (t) = A(sin(2πf p1,k t) + sin(2πf p2,k t)) at a frequencies f p1,k and f p2,k is transmitted simultaneously with the probe excitation to the medium. The modulated coda probe signal y m,k (t) is also recorded. The added feedback consists in evaluating and optimizing the CC [START_REF] Houhat | Optimal pump excitation frequency for improvement of damage detection by nonlinear vibro acoustic modulation method in a multiple scattering sample[END_REF]. This coefficient quantifies the resemblance between the reference probe signal and the modulated probe signal by the pump signal. The higher the nonlinear modulation between the probe and the pump signals is, the lower the value of the correlation coefficient is.

The genetic algorithm is used to find a new pumping frequencies f p1,k+1 and f p2,k+1 minimizing the CC [START_REF] Haupt | Practical genetic algorithms[END_REF]. The frequencies are then modified and all the process described above is reiterated until the algorithm converges toward the best solutions. The first step (called generation 1) consists of choosing randomly N pumping frequencies couples from a uniform distribution on a given frequencies interval. In our case, we have chosen 10Hz ≤ f p1,p2 ≤900Hz. This choice is directly related to the experimental setup. The correlation coefficient is evaluated for each pumping frequencies couple and sorted in descending order. To prepare the next step, the N/2 best pumping frequencies couples that minimize the correlation coefficient are kept for the next generation k+1 and become parents. N/2 new pumping frequencies couples named offspring are generated by the crossover operator by mixing the best previous pump frequencies couple with one of the remaining better couples. Finally, 40% of the pumping frequencies couples were mutated to obtain a robust optimization. In our case the population size N=6.

III. MATERIAL

A. Specimen Description

An aluminium bar (600 mm ×15 mm× 3mm), with density ρ=2700 kg/m 3 , Poisson ratio σ=0.33, and Young modulus E=69 GPa, is used as the specimen for our experiments. The probe wave propagation velocity in the bulk of the bar is estimated to 4838 m/s. 10 tapped holes with a 4 mm diameter and localized at distances of 10 mm, 30 mm, 50 mm, 100 mm and 200 mm from each sides of the bar center are drilled. Identical screws can be placed in the tapped holes to mimic nonlinear solid contacts (cracks) such as in [START_REF] Legland | Evaluation of crack status in a meter-size concrete structure using the ultrasonic nonlinear coda wave interferometry[END_REF]. These holes constitute linear scatterers when no screw is present. According to the number of nonlinear scatterers on the bar, different levels of "effective" damage can be obtained. In the present study, two configurations have been more particularly studied: the first one (Config 1) corresponds to the case where all screws are placed (Fig. 1). Configuration 2 (Config 2) corresponds to the case where only 2 screws are positioned at S2 and S9. Note that the nonlinearity level in the sample depends on the screw number and the nuts tightening.

B. Experimental setup

The closed loop pump frequency optimization process requires the experimental setup depicted in Fig. 2. 250 kHz central-frequencies for transmitting and receiving the probe signal are glued to the ends of the sample. For the transmission of the probe signal, a 100 mVpp sinusoidal burst of 3 periods of 250 kHz frequency, repeated every 20 ms, was transmitted by a function generator (AFG3022, Tektronix, Beaverton, Oregon, USA) and amplified to 60 dB (100 Vpp) by a power amplifier (Type 2713, Bruël & Kjaer, Naerum, Danemark). Simultaneously, two lower frequencies f p1 and f p2 continuous sine pump signals are generated by the computer-controlled In order to improve the signal to noise ratio, an average of 300 successive acquisitions is carried out, and a coda averaged signal is recorded. Each measurement lasts about 10 s. Both the function generator and the oscilloscope are controlled by MATLAB (Mathworks, Natick, MA, USA). The pump signal is desynchronized from the probe so that the nonlinear effects are correctly distributed over the successive acquisitions and averaged over the acquired signal.

IV. EXPERIMENTAL RESULTS

A. Empirical optimization: Correlation coefficient dependency on the pump frequency

Figs. 3a and3c show the experimental Frequency Response Functions (FRFs) obtained for Configs 2 and 3, respectively, in a frequency range from 10Hz to 900 Hz with a step of 1 Hz. We can see that the magnitude of the resonance peaks and the corresponding frequencies are different for the two configurations. The results show that there is a shift in the natural frequencies and an overall change in the frequency response due to the effect of the added screws. The first experiment so called "empirical optimization" is to check the existence of global or local minima for the correlation coefficient vs. the pump frequency. The prominent observation is the presence of local minima peaks localized at specific vibration frequencies corresponding well to the resonances of the bar, and ρ 1 elsewhere. As illustrated in Fig. 3, the pump frequencies that give minima peaks of ρ coincide well with the resonance frequencies of the sample for the Configs 1 and 2. This confirms the expected effects that at I shows the accurate pump frequency optimization values. The results of the monochromatic pump optimization are also depicted for a comparison purpose. For the Config 1, it can be noticed that the GA converges to f p1,opt =770 Hz and f p2,opt =145 Hz after the 14 th generation with ρ opt = 0.983. For the monochromatic pump excitation, the GA converges to the global minimum of the empirical cost function at f p,opt =22 Hz (Fig 4 .a). For the Config 2, we can see in Fig. 5.a that the GA converges aften the 20 th generation, to f p1,opt =149 Hz and f p2,opt =150 Hz with a corresponding correlation coefficient ρ opt =0.994 (Fig 5 .b). The results of the monochromatic pump optimization obtained previously [START_REF] Houhat | Optimal pump excitation frequency for improvement of damage detection by nonlinear vibro acoustic modulation method in a multiple scattering sample[END_REF] are also depicted in Fig. 5. It can be noticed that the the optimal pump frequency f p,opt =153 Hz is very close to f p1,opt and f p2,opt and to the optimal frequency obtained by the empirical optimization( f p,EO =150 Hz). V. DISCUSSION AND CONCLUSIONS Optimal bifrequential pump excitation was automatically selected in order to maximize the nonlinear modulation effects in VAM tests. Two experimental configurations including different positions and amount of localized have been tested. The main result is that the GA seems to propose only a solutions close to the resonance frequencies of the sample, until it converges to the optimal couple. The corresponding values of the CC are quite lower to those obtained by a single pump frequency. Nevertheless, we have expected that exciting the sample at two optimal pumping frequencies should give rise to a more important de-correlation value. This observation may be due to the GA precision. We think that the CC value may be refined by an hybridization of the GA with a more accurate algorithm like the Descent Gradient Algorithm. Also, as a perspective, the optimization of the phase between the two sinusoidal pump will be considered. To summarize, an automatic adaptation of a bifrequential pump excitation by a closed loop VAM system giving the optimal damage detection is proposed. The optimal pumping frequencies corresponded to a resonance frequencies of the studied sample.

Fig. 1 .

 1 Fig. 1. Specimen schemes of the two configurations. One Aluminium bar (600 mm ×15 mm× 3mm) containing 10 tapped holes. Config 1: all the screws are placed on the bar. Config 2: 2 screws are placed at S2 and S9.

Fig. 2 .

 2 Fig. 2. Experimental setup for the closed loop Vibro-Acoustic Modulation (VAM) technique.

Fig. 3 .

 3 Fig. 3. Comparison between the resonance frequencies obtained from the experimental Frequency Response Function of the Config 2 (a) and the Config 1 (c). Experimental results of correlation coefficient vs. pump frequency at the sample resonance frequencies areas for the two sample configurations: b) Config 2 and d) Config 1. The FRF peaks coincide with the minimum peaks of the correlation coefficient vs. pumping frequency.

  Figs.4 and 5show the genetic algorithm optimization results of the pump frequencies for the Config 1 and Config 2 respectively. TableIshows the accurate pump frequency optimization values. The results of the monochromatic pump optimization are also depicted for a comparison purpose. For the Config 1, it can be noticed that the GA converges to f p1,opt =770 Hz and f p2,opt =145 Hz after the 14 th generation with ρ opt = 0.983. For the monochromatic pump excitation, the GA converges to the global minimum of the empirical cost function at f p,opt =22 Hz (Fig 4.a). For the Config 2, we can see in Fig.5.a that the GA converges aften the 20 th generation, to f p1,opt =149 Hz and f p2,opt =150 Hz with a corresponding correlation coefficient ρ opt =0.994 (Fig 5.b). The results of the monochromatic pump optimization obtained previously[START_REF] Houhat | Optimal pump excitation frequency for improvement of damage detection by nonlinear vibro acoustic modulation method in a multiple scattering sample[END_REF] are also depicted in Fig.5. It can be noticed that the the optimal pump frequency f p,opt =153 Hz is very close to f p1,opt and f p2,opt and to the optimal frequency obtained by the empirical optimization( f p,EO =150 Hz). Fig 5.b shows that the CC value obtained for the monochromatic pump is quite similar then those obtained for the bifrequential pump excitation. The GA converges to the global minimum of the empirical cost function (Fig 3.b), for both f p1 and f p2 , in a pump frequency range between 10 and 900 Hz, including four (4) local minima of the cost function. Moreover, until the optimization process, the GA stucks on values very close to the resonance frequencies of the sample. We can see clearly, in Fig 4.a, that f p1 = 731 Hz at the first generation. It is equal to 174 Hz from generations 2 to 4, then 365 Hz until generation 12. Finally, it reaches the optimal frequency at =770 Hz.

  (Hz) fp,opt (Hz) fp1,opt (Hz) fp2,opt (Hz) fp,EO (Hz) fp,opt (Hz) fp1,opt(Hz) fp2,opt (Hz)

Fig. 4 .

 4 Fig. 4. Genetic algorithm (GA) optimization results for the Config 1: (a) The pump frequencies vs. generations. The GA converges to fp,opt=22 Hz after the 7 th generation for the monochromatic pump. For the bifrequential pump, f p1,opt =145 Hz and f p2,opt =770 Hz. (b) The corresponding correlation coefficient ρ vs. generations for the bifrequential and the monochromatic pump excitation.

Fig. 5 .

 5 Fig. 5. Genetic algorithm (GA) optimization results for the Config 2: (a) The pump frequencies vs. generations. The GA converges to fp,opt=153 Hz after the 4 th generation for the monochromatic pump. For the bifrequential pump, f p1,opt =149 Hz and f p2,opt =150 Hz. (b) The corresponding correlation coefficient ρ vs. generations for the bifrequential and the monochromatic pump excitation.

TABLE I THE

 I PUMP FREQUENCIES OBTAINED FROM THE MONOCHROMATIC AND THE BIFREQUENTIAL PUMP OPTIMIZATION, AND THE GLOBAL MINIMUM OF THE EMPIRICAL COST FUNCTION f p,EO .


Andrea Mendizabal
email: andrea.mendizabal@inria.fr

Eleonora Tagliabue
email: eleonora.tagliabue@univr.it

Jean-Nicolas Brunet

Diego Dall'alba

Paolo Fiorini

Stéphane Cotin

Physics-based Deep Neural Network for Real-Time Lesion Tracking in Ultrasound-guided Breast Biopsy

Keywords: Ultrasound-guided Breast Biopsy, Deep Neural Networks, Real-time Simulation

come   L'archive ouverte pluridisciplinaire

Introduction

Breast biopsy is the preferred technique to evaluate the malignancy of screeningdetected suspicious lesions. To direct the needle towards the target, biopsy procedures are performed under image guidance, normally done with ultrasound (US) probes due to their ability to provide real-time visualization of both the needle and the internal structures [START_REF] O'flynn | Image-guided breast biopsy: state-ofthe-art[END_REF]. However, proper needle placement with US remains a challenging task. First, malignant lesions cannot always be adequately visualized due to the poor image contrast of US. Furthermore, navigation towards complex 3D lesion geometries is commonly achieved using 2D freehand US (FUS) systems, which provide information in a lower-dimensional space [START_REF] Krücker | Clinical utility of real-time fusion guidance for biopsy and ablation[END_REF]. Since highly sensitive pre-operative images (such as MRI or CT) can provide accurate positions of the lesions, finding a method to update these positions from real-time US images during an intervention would highly benefit current biopsy procedures. Several commercial and research platforms have implemented image fusion techniques that align pre-operative and intra-operative data, exploiting rigid or affine registration methods [START_REF] Guo | Ultrasound imaging technologies for breast cancer detection and management: A review[END_REF]. However, when dealing with breast anatomy, large deformations arise due to compression forces applied by the US probe. To provide accurate probe-tissue coupling and acceptable image quality, an appropriate alignment procedure of the pre-operative and US data is required.

Accurate modelling of soft tissue deformation in real-time is a far-from-beingsolved problem. Biomechanical models relying on the finite element method (FEM) realistically calculate soft tissue deformations by using a mathematical model based on continuum mechanics theory. Although these models have been successfully employed for multimodal breast image registration, they have never been applied to registration between pre-operative data and intra-operative US, due to difficulties in providing a prediction within real-time constraints [START_REF] Hipwell | A review of biomechanically informed breast image registration[END_REF]. This is especially true when considering large, non-linear deformations which involve hyperelastic objects, as it is the case for the breast.

In order to meet real-time compliance, various techniques have been proposed to simplify the computational complexity of FEM. Some of them have focused on optimizing linear solvers (the main bottleneck of FEM) or the formulation itself, such as corotational [START_REF] Felippa | A unified formulation of small-strain corotational finite elements: I. theory[END_REF] and multiplicative jacobian energy decomposition [START_REF] Marchesseau | Multiplicative jacobian energy decomposition method for fast porous viscohyperelastic soft tissue model[END_REF]. Very efficient implementations also exist, like Total Lagrangian explicit dynamics (TLED) [START_REF] Miller | Total lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation[END_REF], which can achieve real-time performances when coupled with explicit time integration and GPU-based solvers [START_REF] Joldes | Real-time nonlinear finite element computations on gpu-application to neurosurgical simulation[END_REF]. Another possible option to lower the simulation time is through dimensionality reduction techniques, like Proper Orthogonal Decomposition (POD), where the solution to a highdimensional problem is encoded as a subset of precomputed modes. The most optimized approach used to model breast biomechanics is the one proposed by Han et al. in [START_REF] Han | A nonlinear biomechanical model based registration method for aligning prone and supine mr breast images[END_REF], which relies on GPU-based TLED formulation. Despite the significant simulation speedup achieved, solving the FE system took around 30s, which is still not compatible with real-time. Modelling methods that do not rely on continuum mechanics laws have also been used to approximate soft tissues behavior. Among these, the position-based dynamics (PBD) approach has been used to predict breast lesions displacement due to US probe pressure in realtime, providing comparable accuracy with FE models [START_REF] Tagliabue | Position-based modeling of lesion displacement in ultrasound-guided breast biopsy[END_REF]. However, not being based on real mechanical properties, such model requires an initial optimization of simulation parameters to obtain a realistic description of the deformation.

An emerging approach which has the potential of being both accurate and fast, exploits neural networks to estimate soft tissue behavior. Machine learningbased methods have proven successful to predict the entire 3D organ deformation starting either by applied surface forces [START_REF] Morooka | Real-time nonlinear fem with neural network for simulating soft organ model deformation[END_REF][START_REF] Tonutti | A machine learning approach for realtime modelling of tissue deformation in image-guided neurosurgery[END_REF] or by acquired surface displacements [START_REF] Pfeiffer | Learning soft tissue behavior of organs for surgical navigation with convolutional neural networks[END_REF][START_REF] Brunet | Physics-based deep neural network for augmented reality during liver surgery[END_REF]. Being networks trained with synthetic data generated from FE simulations, they can reproduce a realistic physics-based description of the organ mechanical behavior. Using FE simulations for model training in the context of MRI-US deformable image registration has already been proposed in [START_REF] Hu | Mr to ultrasound registration for image-guided prostate interventions[END_REF], where the authors build a statistical model of prostate motion which can account for different properties and boundary conditions. In the case of the breast, the potentiality of employing machine learning techniques has been already shown in [START_REF] Martínez-Martínez | A finite element-based machine learning approach for modeling the mechanical behavior of the breast tissues under compression in real-time[END_REF], where several tree-based methods have been employed to estimate breast deformation due to compression between biopsy plates. These methods have been trained on 10 different patient geometries with a very specific FE simulation, where the upper plate is displaced vertically towards the lower one.

Similarly to works in [START_REF] Pfeiffer | Learning soft tissue behavior of organs for surgical navigation with convolutional neural networks[END_REF][START_REF] Brunet | Physics-based deep neural network for augmented reality during liver surgery[END_REF], we propose an approach where a neural network is trained to predict the deformation of internal breast tissues starting from the acquired surface displacements induced by the US probe. Our network can be seen as a patient-specific model. We train it on a single patient geometry before surgery, with a relatively small amount of training data. However, in contrast to the work of [START_REF] Martínez-Martínez | A finite element-based machine learning approach for modeling the mechanical behavior of the breast tissues under compression in real-time[END_REF], FE simulations that compose the training set are generated with several random input displacements, making our approach able to generalize to different probe positions and compression extents.

The proposed method consists in a U-Net architecture, described in Sect. 2.2, and an immersed boundary method used for generating patient-specific simulations, described in Sect. 2.3. Results presented in Sect. 3 show the efficiency of the method when applied to both synthetic and ex vivo scenarios. Our contribution consists of a novel method to generate a real-time capable soft tissue model to improve target visualization during needle-based procedures. The position of lesions identified beforehand on pre-operative images can be updated from intra-operative ultrasound data and visualized by the surgeon in real-time.

Methods

This work presents a data-driven method to estimate in real-time the displacement of the breast internal structures due to probe pressure during US scanning. In our pipeline, we assume to have a patient-specific geometric model of the breast, obtained from pre-operative imaging such as MRI, and to know the position and orientation of the US probe at each time, thanks to a spatial tracking system. If the tracking coordinate system and the coordinate system of pre-operative imaging are registered, knowledge about the 3D pose and the geometry of the US probe directly allows to identify the contact surface between the breast and the probe. Since the US probe is represented as a rigid body, we can reasonably assume that when the anatomy is deformed by the probe during the image acquisition process, points on the breast surface below the US probe will be displaced to the same exact extent as the probe itself. As a consequence, our method can predict the displacements of all the points within the anatomy given as input the displacement of the surface nodes in contact with the US probe. The decision of relying on surface displacement inferred from the spatial tracking of the US probe instead of directly tracking surface deformations (through, for example, an RGBD camera) was taken from the fact that probe-induced deformations are large but local, and the probe itself would occlude most of the deformed surface to the sensor, thus preventing an accurate estimation of the contact surface displacements.

The U-Net architecture

The objective of our work is to find the relation function f between the partial surface deformation under the US probe and the deformation inside the breast. Let u s be the surface deformation and u v the volumetric displacement field. In order to find f a minimization is performed on the expected error over a training set

{(u s n , u v n )} N n=1 of N samples: min θ 1 N N n=1 f (u s n ) -u v n 2 2 (1)
where θ is the set of parameters of the network f . We propose to use the same architecture as in [START_REF] Brunet | Physics-based deep neural network for augmented reality during liver surgery[END_REF], that is a U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] adapted to our application (see Fig. 1).

The network consists of an encoding path that reduces the high dimensional input into a reduced space, and a decoding path that expands it back to the original shape. The skip connections transfer features along matching levels from the encoding path to the decoding path through crop and copy operations. As Fig. 1 shows, the encoding path consists of k sequences (k = 3 in our case)

of two padded 3 × 3 × 3 convolutions and a 2 × 2 × 2 max pooling operation. At each step, each feature map doubles the number of channels and halves the spatial dimensions. In the lower part of the U-Net there are two extra 3 × 3 × 3 convolutional layers leading to a 1024-dimensional array. In a symmetric manner, the decoding path consists of k sequences of an up-sampling 2 × 2 × 2 transposed convolution followed by two padded 3 × 3 × 3 convolutions. At each step of the decoding path, each feature map halves the number of channels and doubles the spatial dimensions. There is a final 1 × 1 × 1 convolutional layer to transform the last feature map to the desired number of channels of the output (3 channels in our case). The design of the U-Net is based on a grid-like structure due to this up-and down-sampling process. Hence we directly mesh our deformable object with regular hexahedral elements as explained in the next section.

Simulation of breast tissue using hexahedral grids

The training data set consists of pairs of (u s , u v ) where u s is the input partial surface displacement and u v is the volumetric displacement field. Even though the data generation process takes place in an offline phase, in order to generate enough training data with FE simulations within clinically acceptable times (the intervention can be performed on the day after pre-operative scan is acquired), it is important to have simulations that are both accurate and computationally efficient. We consider the boundary value problem of computing the deformation on a domain Ω under both Dirichlet and Neumann boundary conditions. Let Γ be the boundary of Ω (in our case, Γ corresponds to breast external surface, while Ω represents the entire breast volume). We assume that Dirichlet boundary conditions are applied to Γ D and are a-priori known, whereas Neumann boundary conditions are applied to Γ N , a subset of Γ that represents probe-tissue contact area and changes depending on current US probe position. In this work, training data for the network are generated by solving the discretized version of the following boundary value problem, exploiting the FE method:

   -∇ • σ = 0 in Ω u = 0 on Γ D σn = t on Γ N (2)
where σ is the Cauchy stress tensor, n is the unit normal to Γ N and t is a traction force applied to the boundary. Note that in [START_REF] Cignoni | MeshLab: an Open-Source Mesh Processing Tool[END_REF] we neglect all timedependent terms and we do not apply any body force like gravity, since our geometric model already accounts for the effect of gravity force. The relation between stress and strain is described through the Saint Venant-Kirchhoff model, which is the simplest and most efficient extension of a linear elastic material to the nonlinear regime. This choice is motivated by the fact that a simple linear elastic model would not be able to appropriately describe the large deformations undergone by the breast. An iterative Newton-Raphson method is used to solve the non-linear system of equations approximating the unknown displacement.

We choose to discretize the domain into 8-node hexahedral elements not only for their good convergence properties and lock-free behavior, but also because it is the required structure for the input to the network. To do that, the 3D breast geometry is embedded in a regular grid of hexahedral elements (see Fig. 2) and we use an immersed-boundary method to correctly approximate the volume of the object in the FE method computations.

Data generation

The input to the network corresponds to the displacement u s of the points belonging to the breast-probe contact area. The punctual displacements are spread to the nodes of the surrounding cuboid cell through a barycentric mapping and the corresponding volume displacement u v is obtained by the previously explained FE approach in response to u s . The data used to train the network must be representative of the application scenario and must allow the network to extract the pertinent features of the tissue behavior. In order to train our model to estimate breast volume deformation in response to pressure imposed with the US probe, we simulate several random probe-induced deformations using the following strategy:

-Select a random node p in the breast surface -Select an oriented bounding box A centered in point p and normal to the breast surface, whose dimensions match those of the US probe lower surface, which represents current probe-tissue contact area -Select all the surface points P falling within the box A -Select as force direction d the normal to the surface at point p plus a random angle α (α ∈ -π 4 , π 4 ) -Apply the same force f of random magnitude (|f | ∈ [0.0, 0.8]) along direction d to the P selected points simultaneously -Store the displacement at the set of points P (input to the network) and the displacement of all the points in the volume (output to the network)

-Repeat the procedure until N + M samples are generated

The choice of applying force f allowing some angle deviation from normal direction enables us to include in our dataset samples where the probe compression is not precisely normal to the surface, as it can be the case in freehand US acquisitions. The maximal force magnitude (e.g. 0.8 N ) is set such that the amount of maximal deformation reproduced in the training dataset never exceeds too much that observed in real clinical settings. The described strategy is used to generate the set {(u s n , u v n )} N n=1 of N samples which is used to train the network, and the set {(u s n , u v n )} M n=1 of M samples which is left for validation. The training dataset is generated with the SOFA framework [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF] on a laptop equipped with an Intel i7-8750H processor and 16GB RAM.

Experiments and Results

The network presented in this work is used to predict US probe-induced deformations of a realistic multi-modality breast phantom (Model 073; CIRS, Norfolk, VA, USA). The 3D geometry model of the phantom surface and 10 inner lesions (diameter of 5-10mm) is obtained by segmenting the corresponding CT image, relying on ITK-SNAP and MeshLab frameworks [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability[END_REF][START_REF] Cignoni | MeshLab: an Open-Source Mesh Processing Tool[END_REF]. A Freehand Ultrasound System (FUS) based on a Telemed MicrUs US device (Telemed, Vilnius, Lithuania) equipped with a linear probe (model L12-5N40) is used to acquire US images of the 10 segmented lesions. The dimension of the probe surface is (5x1cm). For each lesion, we acquire US images in correspondence of four different input deformations. The MicronTracker Hx40 (ClaronNav, Toronto, Canada) optical tracking system is used to track US probe in space (Fig. 3(a)). The overall probe spatial calibration error is below 1mm (±0.7147), estimated through the PLUS toolkit [START_REF] Lasso | Plus: open-source toolkit for ultrasound-guided intervention systems[END_REF]. Landmark-based rigid registration is performed to refer the CT-extracted 3D model, the US probe and the US images to the same common coordinate system, exploiting 3D Slicer functionalities [START_REF] Fedorov | 3d slicer as an image computing platform for the quantitative imaging network[END_REF]. The registration process does not only enable us to extract the breast-probe contact area, as described in Section 2, but also to know in real-time the 3D position of any point belonging to the US image. In this way, it is possible to refer lesions position extracted from US images to the 3D space.

Predict displacement on synthetic data sets

Elastic properties of the physics model used to generate training data are set in accordance with the values estimated in [START_REF] Visentin | Iterative simulations to estimate the elastic properties from a series of mri images followed by mri-us validation[END_REF] for the same breast phantom considered in this study. However, as we are imposing surface displacements, the values of the elasticity parameters do not affect the displacement field inside the simulated volume as long as the ratio of the different stiffness values is maintained [START_REF] Miller | On the prospect of patient-specific biomechanics without patient-specific properties of tissues[END_REF], thus making the method reliable for any patient specificity. Dirichlet boundary conditions are imposed by constraining the motion of all the nodes belonging to the lowest phantom surface. Using the method described in Sect. 2.2 and 2.3, we discretized the breast phantom into 2174 hexahedral elements and we simulated several probe-induced displacements. Overall we generated N = 800 samples for training and M = 200 samples for testing. The U-Net is trained in a GeForce GTX 1080 Ti using a batch size of 4, 100000 iterations and the Adam optimizer. We used a Pytorch implementation of the U-Net. To assess the learning capability of the network, we perform a statistical analysis of the mean norm error e over the testing data set. Let u v m be the ground truth displacement tensor for sample m generated using the finite element method described in Sect. 2.2 and f (u s m ) the U-Net prediction. The mean norm error between u v m and f (u s m ) for sample m reads as:

e(u v m , f (u s m )) = 1 n n i=1 |u v m i -f (u s m ) i |. ( 3
)
where n is the number of degrees of freedom of the mesh. We compute the average e, standard deviation σ(e) and maximal value of such norm over the testing data set. The obtained results are shown in Table 1. The maximal error is of only 0.266 mm and corresponds to the sample shown in Fig. 4(b). The most striking result is the small computation time required to make the predictions: only 3.14±0.56 ms. In contrast, the FE method takes on average 407.7±64 ms to produce the solution. Obviously, the resolution of the FE mesh could be reduced to accelerate the computations but at the cost of an accuracy loss.

Predict displacement on phantom data

In our experiments, we consider one lesion at a time and we reposition the US probe on the surface of the breast such that the lesion considered is visible on the US image. In order to validate our model, we manually extract lesions position from US image acquired at rest (i.e., without applying any deformation, when the probe is only slightly touching the surface) and we consider it as a landmark to track. We then impose four deformations of increasing extent for each lesion, and we compare the U-Net-predicted displacement with real displacements extracted from US images. The comparison is performed computing target registration error (TRE) between the predicted position of the lesion and its ground-truth position. The performance of our method is compared to that of the FE model used for data generation. In Table 2 are shown the target registration errors for each phantom lesion with respect to the applied deformation. The input deformations are classified into five ranges based on the probe displacements. Displacement ranges indicated as D15, D20 and D25 have a fixed length of 5 mm each and are centered at 15, 20 and 25 mm respectively. D10 and D30 contain the extreme cases under 12.5 mm or above 27.5 mm.

Values in Table 2 highlight that the average TRE for all the tumors and for all the deformations is smaller than 6.194 mm which is comparable to the maximum value obtained with the FE method (6.080 mm). The average error increases with the deformation range just like in the FE method. There is no significant difference between the values of the two tables, meaning that in terms of accuracy, our method is comparable to the data generation method used to train it. In order to compute each deformation, the FE method needs about 407.7 ms whereas the U-Net predicts the deformation in only 3 ms.

Conclusion

In this work we have proposed to use a deep neural network to learn the deformable behavior of the breast from numerical simulations based on the finite element method, in order to bypass the high computational cost of the FEM. Our approach represents an interface between precise biomechanical FE modeling (not capable of real time) and clinical applications requiring both high accuracy and very high speed. We have shown that our framework allows for extremely fast predictions of US probe-induced displacements of the breast during US scanning, achieving comparable accuracy to other existing methods. Therefore, it has the potential to be employed to update in real-time the estimated position of breast lesions identified on a pre-operative scan on US images, en-abling continuous visualization of the biopsy target, even when sonography fails to render it.

Although the FE model used to train our network does not perform in realtime, its prediction delay of less than 1 s might be considered already acceptable for our specific application. However, such good computational performance is achieved since in this preliminary evaluation we use a very simplistic model, that does not account for heterogeneity or complex boundary conditions happening in clinical cases. Usage of a more complex FE model will certainly cause an increase of computation load. On the contrary, an important feature of our approach is that the prediction time remains close to 3 ms regardless of the and of the biomechanical model used for the data generation process. This means that increasing the complexity of the model used to generate the data set will not affect the prediction speed. Moreover, our pipeline allows the method to be insensitive to patient specific elastic properties as it imposes surface displacements. It is worth noting that for inhomogeneous objects, the displacement field still depends on the ratio of the different stiffnesses [START_REF] Miller | On the prospect of patient-specific biomechanics without patient-specific properties of tissues[END_REF]. Another advantage of our method is the easy meshing process. Any geometry can be embedded in a sparse grid and through the use of immersed boundary simulations the deformations are correctly estimated.

The main limitation of our method remains the training process, which is burdersome and has to be repeated for every new geometry or application. However, we have shown that a limited amount of training data can be sufficient to train a U-Net such that it obtains accurate prediction within clinically acceptable times. As a future work, we plan to use a more general training strategy leading to a network model able to predict deformations induced by any type and number of compression tools (for example, different probe shapes or the two biopsy compression plates).

Fig. 1 .

 1 Fig. 1. U-Net architecture for a padded input grid of size 32 × 24 × 16.

Fig. 2 .

 2 Fig.2. Breast surface mesh obtained from a pre-operative CT scan immersed in a hexahedral grid for FEM computations.

Fig. 3 .

 3 Fig. 3. (a)Experimental setup. From left to right: monitor showing real-time US images; CIRS breast phantom during FUS acquisition; optical tracking system that allows to map the real positions of the CIRS breast phantom and the US probe to the preoperative geometry model. (b)External surface and inner lesions of the CIRS breast phantom.

eTable 1 .Fig. 4 .

 14 Fig. 4. (a) Sample with maximal deformation (79.09 mm). (b) Sample with maximal mean norm error (0.266 mm). The green mesh is the U-Net prediction and the red mesh is the FEM solution. The initial rest shape is shown in grey. (c) U-Net prediction on phantom data.

Table 2 .

 2 Target registration errors in millimeters for different tumors and different deformation ranges in the breast phantom. The first table is for the proposed method, while the second table reports results obtained with the FE model used for data generation. Not-acquired data is reported as (-).

			U-Net predictions
	TumorID D10 D15 D20 D25 D30	Mean STD
	1	-	1.936 2.002 1.506 3.053 2.124 0.569
	2	3.211 2.905 4.068	-	4.137 3.580 0.534
	3	2.032	-	4.709 7.134 10.90 6.194 3.262
	4	0.505 2.225 5.313 5.903	-	3.486 2.217
	5	0.932 2.768 3.454	-	4.893 3.012 1.425
	6	3.923 6.349 5.625	-	6.724 5.655 1.075
	7	3.454 3.864 4.543 6.710	-	4.643 1.255
	8	2.422 3.261 4.320 5.136	-	3.785 1.030
	9	-	3.928 4.214 4.578 4.858 4.394 0.353
	10	5.529 3.272 3.940 4.846	-	4.397 0.860
	Mean	2.751 3.390 4.219 5.116 5.761
	STD	1.638 1.294 1.007 1.854 2.788
				FE method
	TumorID D10 D15 D20 D25	D30	Mean STD
	1	-	1.326 2.151 2.075 3.759 2.328 0.887
	2	1.956 2.738 3.945	-	4.025 3.166 0.865
	3	1.595	-	4.748 7.044 10.932 6.080 3.404
	4	0.755 1.991 4.544 5.120	-	3.103 1.795
	5	1.029 2.863 3.330	-	4.541 2.941 1.262
	6	2.579 3.409 2.871	-	2.337 2.799 0.400
	7	2.605 3.219 4.095 6.750	-	4.167 1.582
	8	2.695 2.748 4.321 5.411	-	3.794 1.139
	9	-	2.745 2.497 2.510 4.193 2.986 0.704
	10	2.916 2.542 3.015 3.868	-	3.085 0.485
	Mean	2.016 2.620 3.552 4.682 4.964
	STD	0.765 0.593 0.856 1.803 2.757


Volker Schröder

Marie-Catherine Vs

Le Jumel

Baron Jérôme

The birth and beginnings of Madame d'Aulnoy

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

Her math may be off by a month or so, but "13 years and 8 months" in July 1666 points unambiguously to a birth in the fall of 1652. The writer is indeed already "Dame Donoy," following her marriage to François de la Motte, baron d'Aulnoy, on March 8, 1666. A first daughter was born in January 1667; three more children followed in quick succession, before Marie-Catherine even turned 18.

Barring further discoveries, this marginal inscription appears to be the earliest piece of writing that has come down to us from Madame d'Aulnoy. Proud, feisty, and eloquent, it proclaims the author's identity and warns the future reader of "her book" not to underrate what's "inside." D'Aulnoy's first publications did not appear until 1690, but her handwriting was already on the wall in Barneville in 1666.