V Magnier

E Roubin

J B Colliat

P Dufrénoy

Methodology of porosity modeling for friction pad: Consequence on the squeal instability

Keywords: Squeal, Morphological modeling, Finite Element Analysis, Heterogeneous material

Mechanical property is an important parameter of friction material as it influences contact distribution (tribological and thermo-mechanical behavior), dynamic behavior (and then squeal propensity). Due to the process, the percentage of porosity in the volume is near 10%. So, a friction pad has been characterized by a topographical study using X-ray. This step permits to obtain some information concerning the porosity. In a second step, a simplified pin-on-disc is modeled using finite element method considering initially homogenized friction material and secondly heterogeneous friction material including the pores. The pore distributions are generated using a representative random fields technique. The introduction of pores shows clearly some differences on dynamical behavior between homogenized and heterogeneous materials consideration.

Introduction

A brake system can generate noise emissions. This noise can be classified into different categories based on the frequencies and excitation sources that occur like judder, groan, moan, howl and squeal [START_REF] Akay | Acoustics of friction[END_REF]. The latter remains a major concern of the transportation industry in spite of much research and development work performed over the last few decades. Squeal is characterized by high frequencies ranging from 1 to 10 kHz associated with a high sound pressure level up to 120 dB [START_REF] Akay | Acoustics of friction[END_REF]. Several general and global theoretical models have been developed to understand the acoustic emission mechanisms of sliding systems. Some mechanisms lead to vibrations, like flutter instabilities [START_REF] North | Disc brake squeal -a theoretical model[END_REF], while others sustain vibration, e.g., producing a limit cycle, like stick-slip [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF] and sprag-slip [START_REF] Spurr | A theory of brake squeal[END_REF]. Flutter instability, in a case without damping and gyroscopic forces is the coalescence of two eigenfrequencies, called mode lock-in. The instability is associated with a non null real part of the eigenvalues that are complex due to the asymmetry of the matrix stiffness due to contact terms in the modal analysis. Consequently, mode lock-in implies instabilities in systems and can lead to vibrations in the brake system and noise generation [START_REF] Ouyang | Numerical analysis of automotive disc brake squeal: a review[END_REF] [START_REF] Kinkaid | Automotive disc brake squeal[END_REF].

From a numerical point of view, brake models are mostly carried out with homogeneous and mostly isotropic materials. In the literature, minimal models or models based on the real geometry using finite element method are found to gain understanding of the role of macroscopic parameters like coefficient of friction, mechanical properties, etc. [START_REF] Liu | Analysis of disc brake squeal using the complex eigenvalue method[END_REF], [START_REF] Massi | Brake squeal: Linear and nonlinear numerical approaches[END_REF], [START_REF] Hoffmann | A minimal model for studying properties of the mode-coupling type instability in friction induced oscillations[END_REF], [START_REF] Trichês | Analysis of brake squeal noise using the finite element method: A parametric study[END_REF] etc. Recently, some studies have been performed by introducing heterogeneities in friction pad like [START_REF] Alart | Numerical study of a stratified composite coupling homogenization and frictional contact[END_REF] [24] [START_REF] Mbodj | Modeling of Carbon/Carbon Composites under tribological Solicitations[END_REF] [7] [START_REF] Bonnay | Influence of geometry imperfections on squeal noise linked to mode lock-in[END_REF] where the heterogeneities of friction pad are considered as homogeneous material using homogenization technique. In [START_REF] Pernin | A model of heat transfer in brake pads by mathematical homogenisation[END_REF], a heat transfer analysis of the brake pads with porous materials is performed using mathematical homogenization. To date, most modeling practices consider heterogeneity qualitatively and their choice of input parameters is subjective. Recently, some work [START_REF] Adler | Some new random field tools for spatial analysis[END_REF] has found some applications of a mixture of geometry, topology, probability, and statistics data to many problems.

But it is well known that the lining materials of automotive brakes are usually composites formed by compaction of coarse powders including many different components (typically 20 to 30). The development of friction materials is mainly based on an experimental approach with empirical trial-error methods.

But in nowadays context, this way is long and expensive especially with the norms (environmental, safety etc.) becoming more and more restrictive. The mechanical properties of the components span over a wide range, from very soft and easily worn phenolic resins to hard and wear resistant fibers and abrasive particles. Types and relative amounts of ingredients have been determined by empirical observations including binder resin, reinforcing fibers, solid lubricants, abrasives, fillers and friction modifiers. In [START_REF] Chan | Review of automotive brake friction materials[END_REF], typical brake materials are reviewed and their advantages and disadvantages in contemporary brake applications are discussed. In addition to the basic brake components, some porosity (5 to 10% or more) is normally present from the process. Pores seem to be important to control the tribological and mechanical characteristics as well as other ingredients contained in a brake lining. However the actual effect of the porosity is empirical and not really known specifically its influence on vibration induced by friction and squeal. To consider the description of a porous solid, the principal methods found in the literature are stereology, radiation scattering, pycnometry, adsorption, intrusion, suction, maximum bubble pressure, fluid flow, immersion or adsorption calorimetry, thermoporometry, size exclusion chromatography, Xenon NMR and ultrasonic methods as mentioned by [START_REF] Rouquerol | Recommendations for the characterization of porous solids[END_REF]. For brake pads, some studies have been performed to determine the porous network volume quantification using mercury porosity [START_REF] Rouge | Characterization of the Open Porosity of Brake Pads. II. Correlations Between Volume Porosity and Surface Area Porosity[END_REF]. However, this method is insufficient to analyze, to understand and to develop other characteristics like shape ones, number of pores etc. On the other hand, X-ray makes possible to observe and to analyze the characteristics of pores such as shape, distribution and the positional relation to other materials in the lining. [START_REF] Kimura | Visualization of the Pores of the Brake Lining[END_REF] Kimura et al. analyze the brake lining pores and propose an improvement for more accurate data by a higher contrast resolution image using gallium.

The aim of this paper is to consider explicitly heterogeneity of friction material induced by the porosity on dynamic instabilities. A methodology is proposed starting from a characterization of heterogeneities of a friction pad. The microstructure description is provided by X-ray microtomography. For this study, a reduced formulation friction pad with six components has been chosen. In a second step, a morphological modeling framework for the pores generation using information found in the tomography study is described. This heterogeneous description allows including the porosity in a pin-on-disc model which is performed to identify mode coupling with a complex eigenvalues analysis. Due to the random distribution of pores, some draw of pores dispersal are performed to have a tendency on the dynamical behavior results.

Formulation of friction pad

A typical brake pad material contains quite a number of components, usually 20 to 30. The pad lining materials are categorized as metallic, semi-metallic, organic and carbon-based, depending on their composition. The components are classified into four broad groups: binders, reinforcing fibers, fillers, and frictional additives/modifiers, based on the major function they perform apart from controlling friction and wear performances. The binder hold the ingredients together, to maintain structural integrity of the brake lining under varying mechanical and thermal loading. The structural components provide the mechanical reinforcement of the composite matrix; fillers make up the volume of the brake lining, while keeping the cost down; friction modifiers stabilize the coefficient of friction and wear rate. Apart from these basic considerations, these components perform synergistically in controlling friction and wear performances of the brake pad. Most studies on squeal use standard friction material primarily issued from the automotive industry. It is not rare to find in a commercial friction between 20 and 30 components. This kind of formulation makes difficult a fine analysis of dynamic behavior. For this study, a reduced formulation friction pad has been developed. Here, friction pad contains only six components (see Table 1). This material is developed in collaboration with the company FLERTEX.

For the process, all constituents are firstly accurately weighed and added into mixer. Then, the mixture is loaded into a mold. The molded pad samples are hot pressed (160 o C). The pressure of the press is released several times to let go of the gases generated during the cross linking reaction of the phenolic resin.

The brake samples are subsequently post cured. Figure 1 components introduced before the process.

Another phase appears during the process step: the porosity arises mainly during the degassing step.

Characterization of heterogeneities in friction pad

To characterize the friction pad, a cylindrical sample has been extracted from a friction pad along the thickness. This sample has a diameter of 10 mm and a height of 14 mm. Thereafter, this sample has been submitted to a morphology analysis using X-rays. With this technique, the interactions with matter permit to have information about the composition of the sample, and particularly the porosity. Due to the density of components and the high atom number of the binder (baryte), it is necessary to have a powerful energy to pass through the material. This step has been performed at European Synchrotron Radiation The result of the X-ray tomography is illustrated in figure 2(a) with a height viewport fixed to 4mm where the porosity seems to be important at the mesoscale.

To highlight this feeling, a segmentation of each pictures using two colors is performed using the free software ImageJ [START_REF] Schneider | NIH Image to Im-ageJ: 25 years of image analysis[END_REF]. A 3D rendering of pore distribution is described in figure 2(b).

Firstly in the complete sample, in terms of characterization information, the number of pores is close to 17000 with size which varies from some microns to millimeters. The volume fraction obtained for this sample is close to 10%. Secondly, it is interesting to see that the morphological aspect of porosity is close to a flattened spheroid, with the main axis perpendicular to the pad thickness direction. Figure 3(a) (resp. 3(b)) shows the pore shape for clipping plans in ydirection (resp. z-direction). This special shape is essentially due to the process procedure where a uniaxial pressure on the top of the pad is applied. In view of this result, it is necessary to model the porosity which seems to be of the first order on the mechanical behavior. Indeed, the understanding of the impact of this heterogeneity should permit to improve the process without changing the formulation. Here, the target performance will be the dynamic behavior.

Morphological modeling of pores

In this section, a short presentation of a morphological modeling framework for the pores generation is described. The complete approach is described in [START_REF] Roubin | Modélisation numérique multi-échelles des matériaux cimentaires : apport des modèles morphologiques[END_REF]. Indeed, the porosities information described before is characterized for the sample where the shape and dimensions are fixed. So, a tool to generate a pores distribution for all shapes of the friction pad keeping information found in the tomography study must be developed.

The framework deals with correlated Random Fields (RF), which are random spatial continuous functions. By thresholding such a function, one obtains a so-called excursion that is a binary morphology of random shape. This shape clearly relies on the properties of the chosen RF, as well as on the value of the threshold. Several recent results [START_REF] Adler | Some new random field tools for spatial analysis[END_REF] coming from the geometry of RF create a mathematical and theoretical link between this excursion shape and the RF characteristics, giving the possibility of controlling the excursion as morphological model and to match the morphology extracted from X-ray tomography.

From a numerical point of view, several methods have been developed to simulate correlated RF. Here the Turning-band method is chosen, which transforms a 3D simulation into a 1D case. An analytic formula links the random field characteristics with geometrical and topological quantities (volume, surface area, Euler characteristic...) of the underlying excursion set. The last set of parameters is extracted from tomography description performed before.

For a realistic modeling of a concrete like material, both in terms of geometrical and topological quantities, due to the Gaussian case limitation, an application of the random field is made.

To illustrate the approach, let be a realization γ defined as above in figure 4 and u (dashed red line in figure 4 In a N-dimensional space N+1 LKCs can be defined where each can be thought of measures of the "j-dimensional sizes" of A u . In our three-dimensional case, the four LKCs, denoted by L j , j=0..3, provide both geometrical -L 1 , L 2 , L 3and topological -L 0 -descriptions of the morphology A u . They are defined by:

-L 3 (A u ): the three dimensional volume of A u .

-L 2 (A u ): half the surface area of A u .

-L 1 (A u ): twice the caliper diameter of A u .

-L 0 (A u ): the Euler characteristic of A u , which contrary to the other LKCs is a topological measure.

For example, the Euler characteristic (L 0 ) for a ball or a cube is fixed to 1 but differs from a hollow ball (L 0 =2) or a ring torus (L 0 =0). These L j are specified from tomography data.

The chosen shape friction pad is a parallelepiped of 20*20*10mm 3 . An example of pores generation (grey parts) is illustrated in figure 5 with

L 3 (A u ) = 10%, L 2 (A u ) = 1.8, L 1 (A u ) = 1.

Introduction of pores in finite element analysis

In this section, the finite element model is described. The geometry is inspired by an experimental set-up developed at the Laboratoire de Mécanique de Lille and described further in figure 6(a). This experimental set-up is dedicated to squeal noise study and is focused on the role of friction material by 185 simplifying its dynamic behavior with a reduced number of components. The sensitivity to material formulation and geometry on squeal occurrence has been evidenced by [START_REF] Duboc | Influence of contact conditions and pad geometry on disc brake squeal noise[END_REF]. However the impact of porosity is experimentally very difficult to determine and modeling is required.

Homogenization 190

In this part, the homogenized elastic properties for friction pad included porosity is performed. From now, only pores of more than one hundred micrometers are taken into account due to the limit of the meshing process and CPU time consideration. Indeed, here the pores are explicitly meshed requiring very fine elements (∼ 10µm for mesh size near the pores). The volume fraction of porosity considering with this size and obtained with this fixed sample is 9.6%.

The average strains are introduced as degrees of freedom of Reference Points and Homtools [Homtools] automatically generates the linear constraints between these additional degrees of freedom and the displacements on the boundary of the RVE. Thanks to duality properties, the reaction forces at the Reference Points are the average stresses multiplied by the volume of the RVE. With this method, prescribing average strains or stresses become as easy as prescribing displacements or concentrated forces at a node. Here, the kinematic uniform boundary conditions are used. This can be expressed by the displacement u which is imposed at point x belonging to the boundary ∂V such that:

u i = E ij x j ∀x ∈ ∂V (1)
E ij are the Cartesian components of a given symmetrical second-rank tensor that does not depend on x. This implies that the mean strain over V is:

ε ij = 1 V V ε ij dV = E ij (2)

Description of the pin-on-disc model

The test bench is made of a rigid stand, a thin plate, a pad-housing, a friction pad and a disc as shown in figure 6(a). The friction pad is glued to a pad housing and the latter one is tied with a thin plate. The thin plate extremities are embedded with a rigid block which permits to pilot the displacement.

200

This experimental set-up allows a translation of the pad but also tilting with the flexibility of the thin plate. This experimental set-up has been modeled by the finite element method to obtain the eigenfrequencies of the system. They are To illustrate the strategy of study, the analysis has been firstly carried out with a homogenized and an isotropic friction-pad material.

Dynamic analysis of homogenized case and heterogeneous case without contact

In this section, an eigenvalue analysis is performed on the pin-on-disc previously described without contact between pad and disc. First, a homogenized friction-pad is considered with a value of 2680 MPa for the Young's modulus which is the average value of the material with porosities previously obtained.

The results of the eigenfrequencies associated to eigenmodes are summarized in table 3. Here, only the first three pad eigenmodes are presented in view of their 225 contribution in the configuration performed in section 5.5.

A modal assurance criterion on the pad modes is performed between homogenized and heterogeneous cases and shows a quasi-similarity between modes. The relative error (table 3) shows a quasi-similarity of the eigenfrequencies values between heterogeneous and homogenized friction pad. However,

230
the eigenfrequencies values obtained with the heterogeneous cases are always smaller than the ones obtained with the homogenized case, this tendency can be explained by the difference of the mesh size localized at the friction material.

Indeed, the mesh is the same between the different cases except on the friction pad material where the element size is equal to ∼ 200µm (resp. ∼ 10µm) for homogeneous (resp. heterogeneous) case.

Results with friction pad homogenized material with contact

For this simulation, a homogenized friction-pad is considered with a value of 2680 MPa for the Young's modulus which is the average value of the material with porosities previously obtained. The complex eigenvalues analysis of the complete system has been performed with a variation of the coefficient of friction from 0.1 to 0.8. The complex eigenvalue analysis provides the tool to trace the regions of the parameter space that leads to an instability of the system. Figure 7 shows the complex eigenvalues of the system function of the coefficient of friction. The solid lines represent the complex part i.e. the frequencies while the dashed line represents the maximum real part of the system. According to the results, the eigenfrequencies associated to disc mode are slightly modified. This is due to the high rigidity of the disc relative to stiffness variations of the system influenced by the contact conditions. The pad eigenfrequencies are more affected by the evolution of coefficient of friction. This is mainly attributed to the increase in rigidity induced by the increase in contact length [START_REF] Duboc | Influence of contact conditions and pad geometry on disc brake squeal noise[END_REF]. The tilting pad mode (blue curve) and the twisting pad mode (orange curve) coalesce for a friction coefficient between 0.1 and 0.55, leading to a mode lock-in mechanism and an unstable mode characterized by positive real part (red dashed curve).

Results with heterogeneous friction pad material with contact

In this part, the same problem is solved with an explicit introduction of It can be noted that the eigenfrequencies associated to disc mode are quasi unchanged. But, some eigenfrequencies associated to the pad friction material are clearly modified. Indeed, the tilting pad mode (blue curve in figure 8)

decreases with a drop of ∼ 600Hz while the one associated to twisting pad mode decreases with a drop of ∼ 250Hz in comparison to the results obtained with a homogenized friction pad. So, with the introduction of porosity in the system, any mode lock-in appears, illustrated by a real part equal to 0 for all coefficients of friction (red dashed curve).

Extension with differents draws of pores distribution

In this part, the same problem is solved considering a different draw of pores distribution. In this paper, it is not a question of stochastic analysis but this study permits to know if the tendency is global. Here, one hundred random configurations for porosity distribution are treated using the mathematical tool previously described. In figure 9, the range of variation according to one hundred tested configurations for tilting (resp. twisting) mode is plotted in blue (resp.

orange).

For all tested configurations, the eigenfrequencies associated to disc mode are quasi-unchanged in contrast of the frequencies associated to pad mode. It can be noted that the eigenfrequencies variation associated to tilting and twisting pad mode are around 250 Hz according to different configurations. For all configurations, no mode lock-in appears for all tested cases. The probability distribution of eigenfrequencies are plotted for each coefficient of friction in figure 9 using histograms. Here, the sample mean is equal of 3017Hz for the tilting pad mode (3349Hz for the twisting pad mode) with a standard deviation equal to 66Hz (resp. 76Hz). For the worst case (twisting pad mode), the 95% confidence interval is equal to [3319Hz; 3379Hz]. So, in view of the values of this confidence interval endpoints, the number of random distributions can be sufficient.

This simulation shows that the dynamical behavior of the system is completely affected by the introduction of porosity, even if, as shown, the mechanical properties are quasi-unchanged between homogenized and heterogeneous cases.

Therefore the difference of dynamic behavior cannot be explained by this way. Another possible way to explain this decrease of eigenfrequencies is the modification of the contact pressure distribution.

Interpretation of the results

In figure 10(a), a contact pressure distribution for the homogenized case with a coefficient of friction equal to 0.3 is illustrated. A classic distribution 300 of the contact pressure is due to the friction sliding in the tangential direction and to the static bending deformation of the disc, because of the normal load.

Therefore, a more important contact pressure is located at the leading edge, and to the inner contact radius. In figure 10(b), two contact pressure distributions for heterogeneous cases with the same value of coefficient of friction are plotted.

305

It can be noted that the contact pressure distribution is strongly affected and can lead to change the dynamic behavior of the system. The contact pressure is logically less continuous in the heterogeneous cases, but also less concentrated at the leading edge. Indeed, the presence of pores, especially the ones near the contact surface, modifies the distribution of the contact pressure. In fact, the local stiffness is equivalent or smaller that found with the homogeneous case which permits regarding to the global system to reduce the rotation of the friction pad in the heterogeneous cases. Moreover, the alteration of the pressure distribution can explain the dispersion of the frequencies in function of the geometric pores raw where systematically the frequencies (associated the tilting and twisting modes) obtained in the heterogeneous cases are smaller than the ones found in the homogeneous case. In function of the geometric pores raw, the rotation of the pad can be more or less important.

Finally, the last results of contact pressure could explain the change of frequencies level. This kind of work can be linked to other such as [START_REF] Wang | Experimental and numerical studies of friction-induced vibration and noise and the effects of groove-textured surfaces[END_REF] who shows the effect of groove on squeal by both experimental and numerical methods where the pressure distribution is strongly impacted that can change the frequencies of the problem.

Conclusions and prospects

In this paper, a complex eigenvalue analysis of a pin-on-disc system including heterogeneous friction material with porosity has been developed. The porosity has been firstly identified at the ESRF synchrotron using X-ray tomography. A framework study of morphological generation of pores has been secondly detailed and used to include porosity in a finite element model. Results show clearly that porosity consideration has an impact on the dynamic behavior. Comparing to a homogenized case which is mode lock-in, the introduction of heterogeneities leads to mode lock-out even if the homogenized compressive modulus is the same range. The way of explanation is focused on the difference of contact pressure distribution on the pin contact surface which is more evenly distributed for the heterogeneous cases. This can lead to a frequency decrease of the tilting mode involved in the mode lock-in. In prospects, it would be interesting to describe the other components of the formulation to extend the heterogeneous description.

For that, it is necessary to use an advanced finite element technique to describe all components. This approach could also be extended to a full brake pad of complete brake disc configuration.

Facility. Figure 2 Figure 1 :Figure 2 :

 212 Figure 1: SEM of different components: (a) Steel fibers, (b) Mineral fibers, (c) Twaron fibers, (d) Resin, (e) Natural graphite, (f) Barite

  115

Figure 3 :

 3 Figure 3: Tomography images illustrating pore shape: (a) Clipping plan y=0, (b) Clipping plan z=10mm

Figure 4 :

 4 Figure 4: Illustration of the thresholding technique in 1D and 2D cases.: (a) Thresholding technique in 1D case, (b) Thresholding technique in 2D case

  5 and L 0 (A u ) = 17000.

Figure 5 :

 5 Figure 5: Generation of porosity in a parallelepiped friction-pad using tomography data extracted from an elementary cylinder.

Figure 6 :

 6 Figure 6: Experimental set-up and the associated finite element model + Boundaries conditions: (a) Experimental set-up, (b) Finite element model

Figure 7 :

 7 Figure 7: Eigenfrequencies of the pin-on-disc system according to coefficient of friction

  porosity in the friction pad as previously determined. The mechanical property of matrix is fixed to 3000 MPa for Young modulus and 0.3 for Poisson coefficient although the mechanical properties of the pores are considered to zero for Young modulus, Poisson ratio and density. The results of complex eigenfrequencies with the introduction of pores are illustrated in figure 8 for one set of pores distribution.

Figure 8 :

 8 Figure 8: Eigenfrequencies of the pin-on-disc according to coefficient of friction

Figure 9 :

 9 Figure 9: Eigenfrequencies of the pin-on-disc according to coefficient of friction for one hundred random porosity configurations

Figure 10 :

 10 Figure 10: Contact pressure distribution with a coefficient friction equal to 0.3: (a) Pressure distribution at the pin-contact surface for homogenized material, (b) Pressure distribution at the pin-contact surface for heterogeneous material (2 cases)

Table 1 :

 1 shows the different Formulation of the simplified friction pad

	Component	Category
	Resin	Matrix
	Baryte	Binder
	Twaron fibers	Structural materials
	Steel fibers	Abrasive materials
	Mineral fibers Structural and abrasive materials
	Graphite	Lubricant

Table 2 :

 2 Homogenized compressive modulus of the material considering porosity

Table 3 :

 3 Eigenfrequencies for homogenized and heterogeneous cases without contact between pad and disc

Preprint submitted to Tribology InternationalDecember 12, 2016

Acknowledgments

The present research work has been supported by the international Campus on Safety and Intermodality in Transportation, the Nord-Pas-de-Calais Region, the European Community, the Regional Delegation for Research and Technology, the Ministry of Higher Education and Research, and the National Center for Scientific Research. The authors gratefully acknowledge the support of these institutions. The investigations and results reported herein are supported by the National Research Agency (France) under the COMATCO research program.


B Denis

Senneville A Novell

C Arthuis

V Mendes

P A Dujardin

F Patat

‡ A Bouakaz

J M Escoffre

F Perrotin

ANALYSIS OF CONTRAST-ENHANCED ULTRASOUND USING FLUID DYNAMIC

72 pt 1 in 25.4 mm 54 pt 0.75 in 19.1 mm 54 pt 0.75 in 19.1 mm 54 pt 0.

INTRODUCTION

Contrast-enhanced ultrasound (CEUS) is a non-invasive imaging technique extensively used for blood perfusion imaging of various organs [START_REF] Jakobsen | Ultrasound contrast agents: clinical applications[END_REF]. This modality is based on the acoustic detection of gasfilled microbubble contrast agents used as intravascular flow tracers. Recent efforts aim at quantifying parameters related to the enhancement in the vascular compartment using time-intensity curve (TIC), and at using these latter as indicators for several pathological conditions [START_REF] Tranquart | Perfusion quantification in contrastenhanced ultrasound (CEUS)-ready for research projects and routine clinical use[END_REF]. However, this quantification is mainly hampered by two reasons: first, the quantification intrinsically solely relies on temporal intensity variation, the explicit spatial transport of the contrast agent being left out. Second, the exact relationship between the acquired US-signal and the local microbubble concentration is hardly accessible. In this study, we introduce the use of a fluid dynamic model for the analysis of dynamic CEUS (DCEUS) [START_REF] Denis De Senneville | Development of a fluid dynamic model for quantitative contrast-enhanced ultrasound imaging[END_REF], in order to circumvent the two above-mentioned limitations. A new kinetic analysis is proposed in order to quantify the velocity amplitude of the bolus arrival.

RESULTS

Fig. 1 reports typical microbubble transport fields obtained from two dynamic contrast-enhanced experiments (one clip was selected in the control and one in the ligature population, respectively). The estimated microbubble transports are visually larger for the control population (Fig. 1a-c) as compared to the ligature population (Figs. 1d-f). The spatio-temporally averaged microbubble transport amplitude was significantly higher (p-value=1.1 × 10 -3 , AUROC=0.93) in the control group (9 rats) than in the ligature group (11 rats) .

This work was supported by grants from Inserm (FP, AB), the University F. Rabelais of Tours (FP, AB), the AERO PERINAT Association (FP), the French Society of Perinatal Medicine (FP, CA) and the Laboratory of Excellence TRAIL ANR-10-LABX-57 (BDDS). This study has also been carried out with the financial support of the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx Bordeaux-CPU (ANR-10-IDEX-03-02) (BDDS). Fig. 1. Example of microbubble transport estimates from the DCEUS of two rats using the proposed approach. The first and the second rows display results associated to a control and a ligature rats, respectively. The 1 st , 2 nd and 3 rd columns display DCEUS images acquired 10 s, 15 s and 20 s after microbubble arrival, respectively.

CONCLUSION

This study introduces the use of a fluid dynamic model for the analysis of DCEUS. The proposed classification criterion, estimated from DCEUS echography, was demonstrated to be a good binary classification criterion for ligature/non-ligature rat placentas. Our methodology opens great perspectives for the evaluation of the proposed technique for the clinical diagnostic of obstetrical disorders.


Baudouin Denis De Senneville

Anthony Novell

Chloé Arthuis

Vanda Mendes

Paul-Armand Dujardin

Frédéric Patat

Ayache Bouakaz

Jean-Michel Escoffre

Franck Perrotin
email: franck.perrotin@med.univ-tours.fr

Development of a fluid dynamic model for quantitative contrast-enhanced ultrasound imaging

Keywords: Ultrasound, Perfusion imaging, Computer-aided detection and diagnosis, Quantification and estimation. I

Contrast-enhanced ultrasound (CEUS) is a noninvasive imaging technique extensively used for blood perfusion imaging of various organs. This modality is based on the acoustic detection of gas-filled microbubble contrast agents used as intravascular flow tracers. Recent efforts aim at quantifying parameters related to the enhancement in the vascular compartment using time-intensity curve (TIC), and at using these latter as indicators for several pathological conditions. However, this quantification is mainly hampered by two reasons: first, the quantification intrinsically solely relies on temporal intensity variation, the explicit spatial transport of the contrast agent being left out. Second, the exact relationship between the acquired US-signal and the local microbubble concentration is hardly accessible.

The current study introduces the use of a fluid dynamic model for the analysis of dynamic CEUS (DCEUS), in order to circumvent the two above mentioned limitations. A new kinetic analysis is proposed in order to quantify the velocity amplitude of the bolus arrival. The efficiency of proposed methodology is evaluated both in-vitro, for the quantitative estimation of microbubble flow rates, and in-vivo, for the classification of placental insufficiency (control vs. ligature) of pregnant rats from DCEUS. Besides, for the in-vivo experimental setup, we demonstrated that the proposed approach outperforms the performance of existing TICbased methods.

ultrasound scattering of blood from the microcirculation (i.e. arterioles, venules, capillaries) as well as general circulation (i.e. arteries, veins). The rheology of microbubbles in the blood circulation is similar to that of red blood cells, thus demonstrating that the microbubbles do not interfere with hemodynamics [START_REF] Lindner | Microvascular rheology of definity microbubbles after intra-arterial and intravenous administration[END_REF] and have a good safety profile in both cardiac, vascular and abdominal ultrasound applications [START_REF] Wei | The safety of definity and optison for ultrasound image enhancement: a retrospective analysis of 78,383 administered contrast doses[END_REF].

To achieve this objective, computer tools are required for the spatio-temporal and quantitative assessment of organ perfusion. Practically, once the bolus injection of the ultrasound contrast agent is performed, the enhancement in the vascular compartment is screened using contrast-specific US-imaging sequences. Recent efforts are focused on the quantification of parameters related to this enhancement and on the use of these latter as indicators for several pathological conditions [START_REF] Tranquart | Perfusion quantification in contrast-enhanced ultrasound (CEUS)-ready for research projects and routine clinical use[END_REF]. To this end, the average intensity within a region of interest (ROI) is analyzed as a function of time, in the form of a time-intensity curve (TIC). The wash-in and wash-out of the contrast agent within the ROI are subsequently quantified by calculating amplitude and temporal parameters, such as the peak enhancement (PE), the wash-in rate (WiR), the time to peak (TTP), the area under the curve during wash-in (WiAUC) [START_REF] Tranquart | Perfusion quantification in contrast-enhanced ultrasound (CEUS)-ready for research projects and routine clinical use[END_REF], [START_REF] Rognin | Parametric imaging of dynamic vascular patterns of focal liver lesions in contrastenhanced ultra-sound[END_REF]- [START_REF] Strouthos | Indicator dilution models for the quantification of microvascular blood flow with bolus administration of ultrasound contrast agents[END_REF]. More recently, pixel-wise TIC-based approaches have been proposed. That way, a spatial kinetic analysis can be subsequently performed within the organ of interest [START_REF] Mischi | Angiogenesis imaging by spatiotemporal analysis of ultrasound contrast agent dispersion kinetics[END_REF], [START_REF] Schalk | 4-d spatiotemporal analysis of ultrasound contrast agent dispersion for prostate cancer localization: a feasibility study[END_REF].

One of the most recent applications in this field is quantitative assessment of placental perfusion [START_REF] Burton | What is the placenta?[END_REF], [START_REF] Hua | Effects of diagnostic contrast-enhanced ultrasound on permeability of placental barrier: a primary study[END_REF]. Insufficient development or remodeling of placental maternal vascular supply is a major cause of placental insufficiency, a unifying pathogenetic concept that connects various obstetric syndromes, particularly intrauterine growth retardation (IUGR) and preeclampsia (PE), which are the leading cause of iatrogenic preterm birth and maternal morbidity [START_REF] Brosens | The "great obstetrical syndromes" are associated with disorders of deep placentation[END_REF], [START_REF] Burton | Rheological and physiological consequences of conversion of the maternal spiral arteries for uteroplacental blood flow during human pregnancy[END_REF]. In contrast to other circulatory beds of the systemic circulation where blood traverses from arteries via capillary beds to veins, placental circulation acts as an open system with the spiral arteries opening in a large lake of blood, namely the intervillous space, with little or no impedance to blood flow. That way, using dedicated small animal ultrasound scanner (Vevo R 2100) and TIC-based post-processing software (Vevo R CQ), a recent study reported on real-time monitoring of microbubble perfusion in the uteroplacental unit, across the gestation of rats [START_REF] Arthuis | New insights into uteroplacental perfusion: quantitative analysis using doppler and contrast-enhanced ultrasound imaging[END_REF]. The enhancement of ultrasound contrast in the mesometrial triangle and in the maternal spiral arteries was visualized 1 s after intravenous injection of contrast agents. These agents were then detected through the maternal vessel in the center of the placenta to the chorionic plate. The perfusion of these agents slowly expanded in the placenta to the basal plate. During the gestation of rats, the blood flow increased from days 14 to 20 in the mesometrial triangle and the placenta. In other recent studies, the CEUS of uteroplacental perfusion in nonhuman primates were conducted using clinical scanners (e.g., Sequoia Systems) and custom-designed dynamic CEUS (DCEUS) analysis program or image post-processing software were provided by the manufacturers (VueBox R , Bracco, Milan, Italy). Keator et al., reported that CEUS was able to detect the primary placental disc and underlying vessels two days earlier than DUS identified endometrial thickening [START_REF] Keator | Contrast-enhanced ultrasound reveals real-time spatial changes in vascular perfusion during early implantation in the macaque uterus[END_REF]. Using CEUS, spatial differences in vascular perfusion between the endometrium, myometrium and endometrial-myometrial zone were observed and quantified. Thus, myometrium displayed a higher blood flow than that detected endometrium and junctional zone.

However, using the TIC-based approach, the quantification is typically hampered by the following two main reasons:

1) First, the quantification intrinsically solely relies on temporal intensity variations, the explicit spatial transport of the contrast agent being left out. 2) Second, the exact relationship between the acquired USsignal and the local microbubbles concentration is hardly achievable. In the present study, we introduce the use of a fluid dynamic model for DCEUS analysis. In order to estimate the microbubble transport (noted V ), we employed the following transport equation in homogeneous environnement:

I t + V • ∇I = 0 (1)
where I denotes the grey level intensity on DCEUS images and I t the partial temporal derivative of I. The left part of this equation is composed by a transient term (I t ) and a transport term ( V • ∇I), which stand for any temporal and spatial grey intensity variations, respectively. The estimated transport field V thus accounts simultaneously for spatio-temporal grey level intensity variations occurred during the dynamic imaging sequence. This feature may in turn be advantageous to address the first above-mentioned limitation of the TICbased approach. Now, let c be a continuous and differentiable function that associates each US-grey level intensity I to its corresponding microbubble concentration c (I). The transport equation of fluid dynamic, applied on concentration-converted US-images, can be expressed as follows:

∂c(I) ∂t + V • ∇c(I) = 0 (2)
Hence, one obtains:

c ′ (I) I t + V • ∇I = 0 (3)
By noticing that c ′ (I) is independent of the unknown V , one can verify that Eq. ( 3), after division by c ′ (I), meet the transport Eq. (1). Interestingly, the latter does not take into account the function c at all, demonstrating thereby the intrinsic independence of this equation to microbubble concentration. This feature may thus be advantageous to address the second above-mentioned limitation of the TIC-based approach.

Our contribution is four-fold:

• A new kinetic analysis is proposed in order to quantify the velocity amplitude of the bolus arrival in DCEUS. • A complete pipeline has been established, including motion compensation of the imaged tissue, image prefiltering and estimation of the apparent microbubble transport using Eq. (1). • A numerical resolution of the fluid dynamic model, based on a variational approach, has been set. It requires a reduced number of input parameters that need to be calibrated and the employed cost function renders itself compatible with fast numerical schemes, while providing a dense pixel-by-pixel transport field. • The efficiency of the proposed methodology was assessed both in-vitro, for the quantitative estimation of microbubble flow rates, and in-vivo, for the classification of ligature/non-ligature rats' placentas, from DCEUS. In particular, the sensitivity of the algorithm to empirically determined choices of calibration parameters was analyzed. Besides, the in-vivo classification performance was compared to results obtained using the classical TICbased approach.

II. MATERIALS AND METHODS

The proposed processing pipeline takes benefit of simultaneous acquisition of DCEUS and B-mode images as follows. The instantaneous apparent microbubble transport is estimated using the transport equation of fluid dynamic applied on each acquired DCEUS image. B-mode images are dedicated to the manual delineation of the imaged tissue and to the estimation of possible periodic, spontaneous or drift displacements of the tissue. The latter may be induced either by physiological activity (breathing or peristaltic) or by motion of the ultrasound probe.

Since our study has been conducted using a 2D ultrasound imaging system, the equations provided in the current paper refer to the 2D implementation of the algorithms.

A. Proposed kinetic analysis of DCEUS

The proposed processing pipeline, detailed in Fig. 1, includes the following steps: 1) A B-mode image is first used to manually segment the region where the vascularization needs to be analyzed (for example the placenta for the analysis of placenta abnormalities), as described in section II-A1. This region may also exclude surrounding tissues hampering the DCEUS analysis (such as the fetus location in case of DCEUS in the placenta); 2) Physiological activity and motion of the ultrasound probe are first compensated on the DCEUS images (Task #2), since the latters may hamper the estimation of the apparent microbubble transport (see section II-A2) [START_REF] Klein | Adaptive stochastic gradient descent optimisation for image registration[END_REF] [18]; 3) A spatial low-pass filter is then applied on each DCEUS image (Task #3) in order to mitigate the impact of USspeckles on the employed transport equation of fluid dynamics (see section II-A3) [START_REF] Shapiro | Computer Vision[END_REF]; 4) The instantaneous apparent pixel wise microbubble transport is subsequently estimated (Task #4). The step of the temporal derivative is referred to as δt. The transport between DCEUS images acquired at instants t and t+δt is estimated. Section II-A4 details the proposed numerical implementation [START_REF] Horn | Determining optical flow[END_REF]; 5) The instant time of the microbubble arrival, which for practical reasons is likely to vary from one experiment to the other, is estimated (Task #5), as described in section II-A5; 6) Finally, the quantitative microbubble transport amplitude during the bolus (Task #6) may be assessed, as described in section II-A6.

1) Manual delineation of the imaged tissue (Task #1): A region of interest encompassing the imaged tissue is first manually drawn on a B-mode image, and a binary mask (noted Γ) is constructed. Pixels of the image inside the mask have a value of one, and outside a value of zero.

2) Compensation of the displacement of the imaged tissue (Task #2): Using the proposed microbubble transport estimation method, perturbations occur predominantly when the true transport field violates the physical model used for its estimation. Indeed, any temporal change in image intensity is likely to be attributed to microbubble transport. Therefore, an additional challenge arises from the fact that DCEUS image intensity changes are not only due to contrast change during bolus passage but also to the physiological activity. In the scope of this study, the motion estimation process between t and t+δt was a gradient driven descent algorithm maximizing an inter-correlation coefficient and applied on B-mode images, assuming a translational displacement restricted to the binary mask Γ encompassing the imaged tissue [START_REF] Maintz | A survey of medical image registration[END_REF]. The estimated spatial transformation was subsequently applied to compensate motions of the imaged tissue on the DCEUS image acquired t + δt.

3) Spatial low-pass filtering of DCEUS images (Task #3): A spatial low-pass Butterworth filter (order 1) was then applied on DCEUS images acquired at instants t and t + δt in order to mitigate the impact of local image structures on the estimated microbubble transport, due to the differential nature of terms involved in the transport equation used afterwards (i.e. Eq. ( 1)). Tested cutoff-frequencies for the low-pass filter are provided later in the manuscript.

4) Estimation of instantaneous apparent microbubble transport (Task #4):

a) Implemented variational approach:: The transport model of Eq. ( 1) is intrinsically under-determined, and thus leads to an ill-conditioned numerical scheme. The seek transport field V was thus estimated on a pixel-by-pixel basis through the following minimization process:

arg min V Ω I t + V • ∇I + α ∇u 2 2 + ∇v 2 2 d r (4)
where Ω ⊆ ℜ 2 is the image coordinates domain, (u, v) the estimated pixel wise transport vector components, and r ∈ Ω the spatial location. Note that I t is the partial temporal derivative of I after motion compensation (see section II-A2) and spatial low-pass filtering (see section II-A3). The minimized functional accounts for the two following additive contributions:

• A data fidelity term (left part of the integral in Eq. ( 4)) that optimizes, through a L 1 norm, the transport model of Eq. [START_REF] Correas | Ultrasound contrast agents: properties, principles of action, tolerance, and artifacts[END_REF]. Given that a L 1 penalizer is employed, transient variations act identically regardless the amount of grey level intensity [START_REF] Zachiu | An improved optical flow tracking technique for real-time mr-guided beam therapies in moving organs[END_REF]. • A regularization term (right part of the integral in Eq. ( 4)) designed to provide a sufficient conditioning to the numerical scheme. The regularization term is given by

∇u 2 2 = u 2 x + u 2 y and ∇v 2 2 = v 2
x + v 2 y with u x , u y , v x and v y being the partial spatial derivatives of u and v, respectively. Physically, this regularization term assumes that the transport between neighboring pixels is moderate. α is a pre-defined weighting factor designed to link these two contributions. The interested reader is referred to [START_REF] Corpetti | Dense estimation of fluid flows[END_REF] for additional information concerning the implementation of dense estimation of fluid flows.

b) Numerical resolution:: In order to render Eq. ( 4) differentiable, we replace |s| by ψ(s) = √ s 2 + ǫ 2 , with ǫ = 10 -3 and s = I t + V • ∇I.

Then, by applying the Euler-Lagrange equations on Eq. ( 4) on a pixel-by-pixel basis, one can derive the following system of equations for each r ∈ Ω:

s ψ(s) I x -2α∆u = 0 s ψ(s) I y -2α∆v = 0 (5)
where I x,y denote the partial spatial derivatives of I, and ∆(.) the Laplacian operator. Neumann boudaries conditions were employed.

From here, we have a set of 2 × |Ω| non-linear equations with common unknowns u and v. As suggested in [START_REF] Horn | Determining optical flow[END_REF], we approximated ∆(.) in the discrete domain with ∆u = u -u, u being the 3 × 3 local average of u. That way, two additional implicit linear contributions (along u and v, respectively) are obtained. This approximation allows linearizing the system as follows:

Au + Bv = C Du + Ev = F (6)
with:

               A = I 2 x + 2αψ (s) B = I x I y C = 2αψ (s) u -I x I t D = I x I y E = I 2 y + 2αψ (s) F = 2αψ (s) v -I y I t
Solutions u and v can be found through a fixed-point scheme for which A, B, C, D, E and F are explicitly calculated, as follows:

   u n+1 = u n -I x u n Ix+v n Iy+It I 2 x +I 2 y +2α √ (It+Ixu n +Iyv n ) 2 +ǫ 2 v n+1 = v n -I y u n Ix+v n Iy+It I 2 x +I 2 y +2α √ (It+Ixu n +Iyv n ) 2 +ǫ 2 (7)
where n + 1 denotes the new iteration.

Note that the fixed-point scheme of Eq. ( 7) is designed to get the maximum benefit from implicit terms that can be expressed linearly along u and v, while leaving an explicit expression for u, v and

(I t + I x u n + I y v n ) 2 + ǫ 2 . The residual V n+1 -V n 2
was compared to a maximum allowed tolerance of 10 -3 in order to ensure the convergence. c) Coarse-to-fine scheme: The differential nature of terms involved in Eq. ( 1) hampers the estimation of transport of high amplitude. Displacements that are larger than the size of one pixel cannot be estimated. In order to overcome this limitation, we adopted a coarse-to-fine strategy [START_REF] Pratikakis | Robust multiscale deformable registration of 3d ultrasound images[END_REF], which iterated the registration algorithm from a 16-fold down sampled image step by step to the original image resolution.

In addition, an iterative refinement of the transport estimates was performed within each resolution. This implies running the algorithm several times at the same resolution, initializing the displacement field at the current run of the algorithm with the displacement field that resulted during the previous run. In this manner, the stability of the numerical scheme is improved and at the same time a better quality of the estimates is obtained. It was considered that the numerical scheme in Eq. ( 7) converged when the average variation of the transport magnitude from one iteration to the next was smaller than 10 -3 pixels.

5) Determination of the microbubble arrival time (Task #5):

The instant time of the microbubble arrival was then estimated. For this purpose, each dynamically acquired image was iteratively enumerated and the average pixel intensity over the binary mask Γ encompassing the imaged tissue was calculated. Once this value exceeded a typical pre-defined threshold of 1 % of the maximal intensity saturation value, the associated time point is considered as the arrival time (referred to as t 0 throughout the rest of the manuscript). The sensitivity of the above-mentioned threshold is discussed later in the manuscript.

6) Quantitative evaluation of the overall microbubble transport amplitude (Task #6): At this stage, a set of pixel wise transport fields is obtained, each one being associated to each ultrasound image. The spatio-temporal averaged microbubble transport (noted γ), over the imaged tissue and during a time window (covering a duration ∆T , starting from the bolus arrival time t 0 ) can be then calculated as follows:

γ = t0+∆T t=t0 r∈Γ I( r, t) V ( r, t) 2 t0+∆T t=t0 r∈Γ I( r, t) (8)
Eq. ( 8) simplifies all the pixel-wise transport data V ( r, t) down to a single averaged value γ. Note that the amplitude of each transport vector V ( r, t) 2 is weighted by the amount of grey level intensity in the DCEUS image I( r, t). That way, similar γ values are obtained for scenarios with identical microbubble transports, but various pixel intensities and/or various amounts of pixels exhibiting microbubbles. V ( r, t) is converted in mm/s beforehand so as to express γ in a common metric unit.

7) Calibration of the algorithm: At this point, it is important to underline that four user-defined parameters may impact the performance of the kinetic analysis process:

1) The duration of the observation window ∆T , which defines the time interval after bolus arrival in which the kinetic analysis is performed; 2) The step of the temporal derivative δt, which defines the time interval for the estimation of the instantaneous microbubble transport; 3) The cutoff-frequency f c of the low-pass filter applied on the input DCEUS images. Throughout the rest of the manuscript, f c is expressed as a fraction of the original image sampling frequency of the pixels (noted f 0 ); 4) The regularization term of the transport equation α, which infers the stability of the numerical scheme. In the rest of the manuscript, a special attention will be drawn to the impact of these four parameters on the overall results.

8) Hardware and implementation: Our hardware platform was an Intel 2.5 GHz i7 workstation (8 cores) with 32 GB of RAM. The implementation was performed in C++ and parallelized through multi-threading.

B. Experimental evaluation

1) Definity R ultrasound contrast agents: Definity R microbubbles (Lantheus Medical Imaging, Billerica, MA) are a second-generation clinically approved contrast agent, com-posed of octafluoropropane gas encapsulated in a thin and flexible monolayer of PEGylated phospholipids [START_REF] Raymond | Broadband attenuation measurements of phospholipid-shelled ultrasound contrast agents[END_REF], [START_REF] Hyvelin | Characteristics and echogenicity of clinical ultrasound contrast agents: An in vitro and in vivo comparison study[END_REF]. The mean diameter ranges from 1.1 to 3.3 µm. The physicochemical properties of Definity R microbubbles are described in [START_REF] Kitzman | Efficacy and safety of the novel ultrasound contrast agent perflutren (definity) in patients with suboptimal baseline left ventricular echocardiographic images[END_REF]. Definity R microbubbles were prepared according to the manufacturer's instructions. Briefly, a single vial of Definity R was set to room temperature and was then activated using a Capmix R device (3M-ESPE, Cergy-Pontoise, France) for the full 45-second activation cycle.

2) High frequency ultrasound scanner: A MS-250 probe (21 MHz center frequency, 13-24 MHz bandwidth, 75 µm axial and 165 µm lateral resolutions) connected to the Vevo R 2100 ultrasound scanner (VisualSonics Inc., Toronto, Canada) was used to acquire in-vitro and in-vivo images.

3) In-vitro flow experiments:

a) Evaluation platform: A laboratory-made flow system was an open circuit consisting of 3.59 mm internal diameter flexible silicone tubing, through which a suspension of Definity R microbubbles circulated. The tubing flow system was submerged in a degassed water tank. One end was connected to reservoir filled with degassed water, in which Definity R microbubbles were diluted in physiological serum. The diluted microbubble solution was subsequently delivered by a peristaltic pump (MCP Process IP65, Cole-Parmer GmbH, Wertheim, Germany). The center of the flow tube was positioned at a distance of 21 mm from the MS-250 probe. For each experimental condition, a video clip of 30 s was recorded at 20 frames per second (pixel size=0.027 × 0.027 mm 2 ) using the Vevo R 2100 scanner.

b) Assessment of the estimated DCEUS-based microbubble transport:

After image acquisition, microbubbles perfusion was quantitatively analyzed as follows. The DCEUSbased microbubble transport (γ) was computed for four regularly sampled flow rates delivered by the peristaltic pump (i.e. 1.2, 2.4, 3.6 and 4.8 mL/min). All experiments were independently replicated (N = 3). The DCEUS analysis was performed on a mask Γ (see Task #1 in Fig. 1) encompassing the microbubble circulation area (i.e. the inner region of the tube, no wall). In order to assess quantitatively the relationship between the DCEUS-based microbubble transport and the delivered flow rate, a correlation coefficient (R 2 ), the slope (which is expected to be close to 1) and the y-intercept (which is expected to be close to 0) of a linear regression were calculated. The complete evaluation process was repeated for two different microbubble concentrations: 1:1000 and 1:2000 (diluted in physiological serum). These dilution values were selected on the basis of high scattering of the microbubbles and low attenuation. Note that these values are only slightly different from clinical situations where a volume of 1 or 2 mL of microbubbles is diluted in about 5 liters of blood. The delivered flow rates were converted into flow velocities (expressed in mm/s) as follows: considering that the tube diameter is 3.59 mm, the cross sectional area of the tube is equal to 10.12 mm 2 . For a typical volume flow rate of 1.2 mL/min, the plug flow velocity is equal to (1.2 × 1000 mm 3 /min) / (10.12 mm 2 ) = 118.6 mm/min = 1.98 mm/s. The same calculation was done for all other delivered flow rates. 4) In-vivo procedure: All procedures were performed in accordance with the ethical guidelines and approved by the French Committee (n • 19) for Animal Care and Ethics in Animal Experiments (APAFIS#3879-2016020117195710). Twenty pregnant Sprague-Dawley rats were purchased from Janvier Labs CERJ (Le Genest Saint-Isle, France). They were maintained at constant room temperature with 12h light cycle. The rats were 10-12 weeks old at the beginning of the experiments, weighing in average 377 ± 17 g. a) Surgery: As previously described by Wigglesworth et al., a ligature of the uterine horn pedicle was performed on eleven rats to induce an IUGR on the seventeenth days of gestation [START_REF] Wigglesworth | Experimental growth retardation in the foetal rat[END_REF]. Under gaseous anesthesia (Aerrane R r, Baxter, Deerfield, IL), the pregnant rats were placed on a thermostatically controlled pad in order to maintain body temperature at about 37 • C. Before the surgery, analgesia was performed by intraperitoneal (i.p.) injection of buprenorphin (0.05 mg/kg; Buprecare R , Axience SAS, Pantin, France). After shaving of the abdomen, a surgeon performed a midline laparotomy incision with sterile precautions. The number of implantation sites was checked in each uterine horn. A 5-0 Ethilon nylon suture (Ethicon, Somerville, NJ) was placed around the uterine vessels near the caudal end of one horn. The non-ligated horn served as a control. The abdominal incision was repaired in layers using standard surgical procedures. Five hours later, a single 5-mg/kg buprenorphine was intraperitoneally administered to manage postoperative pain.

b) CEUS examination: Under gaseous anesthesia (Aerrane R ), the pregnant rats were placed on a thermostatically controlled pad and their abdomens were shaved before CEUS examination on the nineteenth day of gestation [START_REF] Arthuis | New insights into uteroplacental perfusion: quantitative analysis using doppler and contrast-enhanced ultrasound imaging[END_REF]. A 24-gauge catheter was placed in a tail vein to inject Definity R contrast microbubbles. Ultrasound B-scans were used to image fetal-placental unit in crosssection. A bolus of 200 µL of contrast agent (0.5 mL/kg) was intravenously injected. Subsequently, a video clip of 150 s was recorded at 10 frames per second (pixel size = 0.035 × 0.035 mm 2 ) to investigate the utero-placental perfusion. The 2D US-imaging plan was oriented to cross the central part of the placenta around the cordonal insertion, as preconized in [START_REF] Brosens | The "great obstetrical syndromes" are associated with disorders of deep placentation[END_REF].

c) In-vivo DCEUS analysis: After image acquisition, utero-placental perfusion was quantitatively analyzed from the DICOM video data using the proposed methodology. The perfusion was also analyzed using the existing TIC-based approach. For this purpose, we used the CEUS analysis software Vevo-CQ TM , which is directly integrated in the ultrasound scanner. The four following parameters were extracted from the TIC: PE, WiR, TTP and WiAUC. The interested reader is referred to [START_REF] Tranquart | Perfusion quantification in contrast-enhanced ultrasound (CEUS)-ready for research projects and routine clinical use[END_REF] for a complete description of the above mentioned TIC-based parameters. Note that, for both of our approach and tested TIC-based methods, the DCEUS analysis was performed on an identical ROI (Γ) encompassing the imaged tissue (see Task #1 in Fig. 1) in order to clarify the benefits on the final results.

For the two rat populations, all indicators (i.e. the overall microbubble transport amplitude γ as well as the 4 parameters extracted from the TIC) were expressed as medians and interquartile ranges, and were compared by using the unpaired Mann-Whitney U test. The results were considered significant when the p-value was lower than 0.025.

In addition, considering each indicator as a classification criterion of the two rat populations (i.e. ligature vs nonligature), the performance of a binary classifier system was also assessed using a "receiver operating characteristic" (ROC) curve. The area under the ROC curve (AUROC) was subsequently computed: while a binary classifier acts like a completely random guess for AUROC = 0.5, the best possible prediction method would yield to a value of 1.

d) Assessment of image-based measurement errors:

Potential errors on γ may be ascribed to the manual delineation process (Task #1) and to the organ displacement estimation (Task #2). Errors were individually applied on the manually delineated mask Γ (the binary mask was iteratively eroded until its surface reached the seek error) and on the organ displacement estimates (the error was added on each component of the 2D translation estimated in section II-A2). For each scenario, the mean, the standard deviation, as well as the first, second and third quartiles of γ were calculated.

III. RESULTS

A. In-vitro flow experiments

Fig. 2 reports typical findings obtained using the employed in-vitro experimental platform. A temporal window size ∆T = 15 s was employed here. In order to supply as much as possible the high inherent velocity of the microbubble in the silicone tubing, the transport was estimated between consecutive images. In other terms, since the imaging framerate was equal to 20 Hz, the temporal derivative step δt was set to a value of 1/20 s. High instabilities of the numerical schemes were observed using the proposed implementation for α<0.3. A value of α=0.3 was therefore selected to generate the presented results. The cutoff-frequencies f c of the spatial low-pass filtering of DCEUS images was set to the original image sampling frequency f 0 (i.e. no image filtering was applied) in order to mitigate the deterioration of the apparent local shapes of regions exhibiting microbubbles. Fig. 2 displays typical examples of DCEUS-based microbubble transport fields (displayed in form of vector fields) estimated using the proposed methodology for microbubble dilutions of 1:1000 (2a) and 1:2000 (2b). Vectors of homogeneous amplitude and oriented along a common horizontal direction (accurately matching that of the delivered microbubble flow) are observable for each tested bubble dilution. Note that Fig. 2a and 2b display identical transport amplitude (the size of yellow arrows is similar in 2a and 2b) but different pixel intensities (image intensity is higher in 2a than in 2b) and different number of pixels exhibiting microbubbles (arrow density is higher in 2a than in 2b). The averaged transport amplitude γ remains however visually similar (see Fig. 2c and2d). As shown in Fig. 2c and2d, the proposed estimated DCEUSbased transport amplitude (γ) is a good quantitative indicator of the flow rate within the tube, for the two tested bubble concentrations (R 2 >0.96, y-intercepts <1 mm/s). Note that a slight overestimation of the flow velocity (slopes > 1.4) stemmed from fluid friction considerations close to the tube walls, the latter being further hampered by the practical design of the pump (our pump contained a funnel located in the center of the imaged tube, outside the FOV, in order to regulate the flow). It can be noticed that, for a bubble dilution 1:2000, the precision of γ diminished, especially for low delivered flow rates (i.e. 1.2 and 2.4 mL/min). The DCEUS-based microbubble estimation process was in this case impaired by few local image contrasts (see the reduced US contrast within the tube in Fig. 2b, as compared to Fig. 2a), which in turn hampered occasionally the conditioning of the employed 8)) estimated over each rat populations using the proposed approach, (b-e): distributions of PE, WiR, TTP and WiAUC estimated over the placenta for the two rat populations using the existing TIC-based approach. The median is shown by the central mark, the first and the third quartiles are reported by the edges of the box, the whiskers extend to the most extreme time points that are not considered as outliers, and outliers are individually plotted in red. Note that the term "outliers" refers here to values which benefitted from a specific display (i.e. red crosses), and not to rats which would have been removed in the statistical analysis. It corresponds to values that are more than 1.5 times the interquartile range away from the top or bottom of the plotted box. (f): ROC curve obtained using γ, (g): ROC curves obtained using PE and WiR, (h): ROC curves obtained using TTP and WiAUC. I: Analysis of the impact of unwanted errors on γ occurred in the targeted organ delineation process (Task #1) and in the organ displacement estimation (Task #2), during the in-vivo experiments. The first raw of each cell reports the mean ± standard deviation over the 20 rats, and the second raw provides the first, second and third quartiles ([0.25, 0.5, 0.75]).

numerical scheme of Eq. ( 7), especially for low-resolution levels in the coarse-to-fine scheme (see section II-A4). The proposed algorithm also failed to find numerically the global optimum of the functional of Eq. ( 4) for high transport data (i.e. for flow rates > 8 mL/min).

B. In-vivo flow experiments

Fig. 3 reports typical microbubble transport fields obtained from two dynamic contrast-enhanced experiments (one clip was selected in the control and one in the ligature population, respectively). For each rat, the mask encompassing the placenta, manually drawn in Task #1, is displayed in Figs. 3a and3e. The estimated microbubble transports are visually larger for the control population (Figs. 3b-d) as compared to the ligature population (Figs. 3f-h) (note that the estimated vector fields follow visually the dynamic contrast-enhancement in associated DCEUS images). This observation is confirmed in Fig. 4: the instantaneous transport of the apparent microbubble was larger for non-ligature rats, especially for images acquired during the first 10 seconds after the bolus arrival in the placenta. The spatio-temporally averaged microbubble transport (a) (b) Fig. 6: Analysis of the impact of the temporal window size (∆T ) and temporal derivative step (δt) on the overall results of the proposed approach. The bolus velocity was estimated for various values of ∆T and δt (α and f c being fixed to typical values of 0.1 and f 0 /16, respectively). The significance of the difference between bolus velocities obtained for the two rat populations (i.e. ligature and control) is reported using the pvalue of an unpaired Mann-Whitney U test (a) and the AUROC (b). amplitude (i.e. coefficient γ in Eq. ( 8)) was significantly higher (p-value = 1.1×10 -3 ) in the control group than in the ligature group (see Fig. 5a). The associated AUROC, displayed in Fig. 5f, was equal to 0.9343. On the other hand, none of tested TIC-based parameters led to a significant difference between control and ligature populations (p-values obtained using PE, WIR, TTP and WiAUC were equal to 0.098, 0.348, 0.972 and 0.285, respectively), as shown in Fig. 5b-e. The best achievable AUROC was obtained using PE (AUROC=0.72), as shown in Fig. 5g and5h.

Fig. 6 and Fig. 7 analyze the sensitivity of the proposed methodology against the four calibration parameters listed in section II-A7. It can be observed in Fig. 6 that γ (i.e. the overall microbubble transport amplitude) is a good classification criterion of the two rat populations for all tested temporal window size (∆T ) and temporal derivative step (δt) (p-value constantly < 10 -2 , AUROC > 0.89), especially for low values of ∆T (i.e. < 20 s). Fig. 7a shows that the performance of γ as . The significativity of the difference between bolus velocities obtained for the two rat populations (i.e. ligature and control) is reported using the AUROC. a classification criterion for the two rat populations improved when reduced cutoff-frequencies f c were employed in the spatial low-pass filtering of DCEUS images. This justifies our choice to set f c = f 0 /16 in the results previously shown. In addition, it is important to report that high instabilities of the numerical schemes were observed for α < 0.1. This justifies our choice to set a value of 0.1 as a default parameter for α for the minimized functional of Eq. ( 4) in the scope of this experiment. Fig. 7b shows that comparable classification performances were assessed for an α value 15 times higher.

Table I analyzes the bias on γ arising from errors applied on Γ (Task #1) as well as on placenta displacement estimates (Task #2). This bias has to be examined in the light of Fig. 5a, which shows that a threshold on γ around 0.05 allowed classifying the two rat populations. A bias of 0.0003 ± 0.0055 mm/s occurred when Γ was eroded by 20 % of its original surface. It can also be observed that 5 % of error in the targeted organ motion compensation process (which still corresponds to a systematic error around 7 pixels with the used pixel size) led to a bias on γ of 0.0034 mm/s.

IV. DISCUSSION

In this study, we introduce the use of a fluid dynamic model for the analysis of DCEUS. A complete kinetic analysis was proposed in order to quantify the velocity amplitude of the bolus arrival. The estimated transport was shown to be highly correlated with the microbubble velocity within our in-vitro system, for the delivered flow rates and for the two tested microbubble concentrations (see Fig. 2). Moreover, the proposed classification criterion γ, estimated from DCEUS, was demonstrated to be a good binary classification criterion for ligature/non-ligature rat placentas. In particular, it outperformed TIC-based methods that provided best results using the parameter PE (AUROC=0.72, as shown in Fig. 5g and5h), but remained however well below our fluid dynamic approach (AUROC=0.9343, see Fig. 5f).

Conceptually, the proposed fluid dynamic approach analyses pixel wise microbubble velocity fields (both spatial and temporal derivative of the DCEUS image intensity are simultaneously involved in Eq. ( 1)). While existing compartment models or TIC-based approaches [START_REF] Tranquart | Perfusion quantification in contrast-enhanced ultrasound (CEUS)-ready for research projects and routine clinical use[END_REF] intrinsically solely rely on temporal intensity variations, our method performs a simultaneous spatio-temporal analysis. Consequently, our method provides an additional information for the kinetic analysis of DCEUS (see the quantitative estimates on the delivered flow rate in Fig. 2, which are not accessible when using the TIC-based approach together with a continuous infusion of microbubbles). That being mentioned, while a difficult conversion from grey level intensity to microbubble concentration is a common necessary prerequisite when using the existing TICbased approach, this task is circumvented here. This advantage of the proposed methodology is illustrated on the presented invitro experiment (see Fig. 2), for which quantitative indicators of the flow rate could be accurately estimated, for the two tested microbubble concentration. Practically, microbubbles generate a hyper-intense signal, which saturates grey level intensities. The histogram of the DCEUS images within the imaged organ consistently depicts two main peaks: while high values correspond to pixels where microbubbles are present, low values are associated to the others. "Binary" DCEUS images are thus obtained which allows the estimation of a wave front displacement using the proposed methodology, regardless microbubble concentration considerations. It can be noticed that the parameter γ in Eq. ( 8) is weighted by image intensity. A density conversion task will thus impact the value of γ in case of a non-linear function c (we recall that c associates each US-grey level intensity I to its corresponding microbubble concentration c(I)). However, as mentioned previously, microbubbles generate in practice a hyper-intense signal, which saturates grey level intensities. c is thus likely to provide near values for saturated pixels. Here again, we therefore anticipate a moderate impact of the density conversion task on γ.

The proposed methodology is sensitive to the four calibration parameters listed in section II-A7. Fig. 6 and Fig. 7 demonstrate that a compromise has to be found for each parameter, as follows. First, hile the duration of the observation window (∆T ) must be high enough to screen a significant microbubble bolus (around 10 seconds were found to be sufficient in the scope of the in-vivo experiment), this value must be also reduced in order diminish potential long-term motion bias (e.g., motion of the ultrasound probe, fetus motion...), as shown in the AUROC map of Fig. 6b. Second, while the step of the temporal derivative (δt) must be high enough to observe a significant microbubble transport between instants t and t + δt (several seconds provided good results in the employed experimental setup), this value must also be reduced in order to diminish the potential occurrence of physiological motions between the same instants. Third, while the cutoff-frequency of a low-pass filter applied on the input images (f c ) has to be low enough to filter regionalized ultrasound artifacts which may be attributed to "false transport" (f c = f 0 /16 provided best results on the presented in-vivo data), this value must also be high enough to prevent discarding accurate microbubble information. Finally, the regularization parameter (α) in the optimized functional (Eq. ( 4)) must be high enough to avoid ill conditioning of the numerical scheme (≥ 0.1 with the proposed implementation for the proposed in-vivo protocol), it must also be low enough to allow estimating local transports. Fig. 6 and Fig. 7 demonstrate that the four calibration parameters ∆T , δt, α and f c may be fixed permanently for a specific DCEUS imaging protocol. Consequently, a great advantage of the proposed methodology is that an intervention of the user is only required for the manual determination of the ROI encompassing the imaged tissue.

Using the proposed approach, it is shown that a moderate imaging session duration is mandatory (< 20 seconds after the bolus arrival in the placenta, as shown in Fig. 6). This is a great advantage since it diminishes the potential occurrence of long-term bias, including motion of ultrasound probe, fetus motion, etc... It must be reported that the use of an image intensity weighting in Eq. ( 8) mitigates the impact of potential errors in the delineation of the placenta (Task #1): errors in regions exhibiting permanently low US-intensities have little impact on γ and can be thus be evenly added/removed (see the moderate variation of γ for different sizes of Γ in Table I). However, it must be underlined that our method failed analyzing the vascularity of the organ of interest if vascular areas are not taken into account in Γ. Another source of uncertainty may come from potential errors in the targeted organ motion estimation process, as shown in Table I. For the general case, a moderate impact on the classification performance may thus be anticipated for a sub-pixel registration accuracy.

A small positive threshold on US-intensities was empirically determined to detect the bolus arrival time t 0 (see section II-A5). The threshold value has a direct impact on t 0 and, in turn, on the computed γ. However in practice, the stiffness of the time intensity slope moderated the impact of the threshold value: over the presented in-vivo experiments, intensity thresholds of 5%, 10% 20% induced averaged t 0 delays of 0 s, 0.15 s and 1.52 s, respectively. These delays had no impact on the corresponding AUROC, which was constantly equal to 0.9343.

It is interesting to note that, in the in-vitro experiments, the chosen pixel size was very small (0.035 mm) as compared the speckle spot size (several millimeters). That way, a spatial low-pass filter could be applied (i.e. Task #3) which mitigated in turn the impact of the pixel size on the overall results. The use of a common cutoff-frequency for the low-pass filter may facilitate the generalization of our method to US-systems with different pixel sizes, array geometries and frequencies.

It must also be underlined that the proposed implementation opens good perspectives toward a real-time diagnostic: the employed variational cost function of Eq. ( 4) renders itself compatible with a fast linear numerical scheme (see Eq. ( 7)), while providing a dense pixel wise transport field. Calculations are consequently low time consuming: on our test platform, less than 100 ms were mandatory for the complete processing of one single frame (i.e. for Tasks #2, #3 and #4, including the coarse to fine scheme described in section II-A4). This computation time was measured when using the entire image FOV in order to provide an upper bound of what is expected for calculations restrained to the mask Γ. Taking into account that the employed imaging frame-rate was equal to 10 Hz for the in-vivo experiments, the processing of one image could be achieved within the interval of time available between two successive acquisitions. An embedded implementation on an ultrasound scanner system for an immediate diagnostic is thus conceivable using available hardware.

Several limitations must be emphasized in this retrospective proof-of-concept study, including the fact that a 2D ultrasound imaging system was employed. This may induce potential out-of-plan motions and may hamper the screening of the bolus arrival in the image field-of-view [START_REF] Ta | 2-tier in-plane motion correction and out-of-plane motion filtering for contrast-enhanced ultrasound[END_REF]. Moreover, the perfusion and bubble circulation is a 3D process and a 2D model thus provides only an incomplete view of the reality. It must be underlined that our choice to perform 2D imaging was purely dictated by the available US-hardware, since the proposed processing pipeline itself fully generalizes in 3D. In our future work, we plan to use a 2D probe to generate 3D data (as shown in [START_REF] Schalk | 4-d spatiotemporal analysis of ultrasound contrast agent dispersion for prostate cancer localization: a feasibility study[END_REF]) and to extend the implemented model to fit 3D measurements. That being said, the absolute need of a 3D process can be moderated for specific clinical scenarios. For example, if one considers the analysis of the perfusion in the placenta, the imaging plan can be oriented to cross the central part of the placenta around the cordonal insertion (as done in our in-vivo experiments). That way, an accurate screening of the perfusion within the placenta can be assessed using a 2D US-probe, as shown in [START_REF] Brosens | The "great obstetrical syndromes" are associated with disorders of deep placentation[END_REF]. It must be also underlined that CEUS quantification is typically hampered by several factors causing variations of the image intensity including propagation/attenuation through tissue and non-linear propagation effects from the microbubble (a good review of these sources of variability and attempts for correction can be found in [START_REF] Tang | Quantitative contrast-enhanced ultrasound imaging: a review of sources of variability[END_REF]).

Moreover several remarks can be made about the stability of the numerical scheme of Eq. ( 7) and its robustness/accuracy with respect to the input data. This directly relies on the US-data content and thus on the specific experimental conditions. In our in-vitro and in-vivo studies, microbubbles are immersed in highly different environments (water for the in-vitro case, placenta tissue for the in-vivo one). While microbubble dilutions were optimized in-vitro to screen clear moving microbubble aggregates in the FOV, no dilution was employed in-vivo in order to generate a distinct hyper-intensive microbubble wave-front. Image contents in turn logically differ between the two scenarios, and the anticipated impact on the system conditioning of Eq. ( 7) is not straightforward. In particular, as reported in section II-A4, the proposed coarseto-fine scheme begins by solving Eq. ( 7) using 16-fold downsampled data. While this is of great advantage on our in-vivo data (which exhibits a local compact region with microbubbles surrounding the artery), this becomes a limitation on our invitro data due to the inherent filtering of the small local regions with microbubbles. This filter deteriorates the information contained in the 16-fold down sampled image, degenerating in turn the system conditioning of Eq. ( 7). The algorithm would benefit from an individual calibration for specific application scenarios. It can also be anticipated that the proposed algorithm may need several improvements/adaptations in term of stability/robustness of the numerical scheme for other possible scenarios [START_REF] Ackermann | Detection and tracking of multiple microbubbles in ultrasound b-mode images[END_REF], including an analysis of the conditioning of the system with an automatic detection of numerical fails [START_REF] Belsley | Regression Diagnostics: Identifying Influential Data and Sources of Collinearity[END_REF].

In addition, the compensated displacement of the imaged tissue is here restrained to translational displacement and may therefore under-approximate real deformations [START_REF] Woo | Non-rigid ultrasound image registration based on intensity and local phase information[END_REF], [START_REF] Cifor | Hybrid feature-based diffeomorphic registration for tumor tracking in 2-d liver ultrasound images[END_REF]. This may in turns results bias in transport field estimated derived from motion compensated DCEUS images. Moreover, although encouraging results were obtained, a decision threshold for the classification criterion γ is needed in order to perform a diagnostic and its sensitivity must be analyzed. Our ultimate goal is to evaluate the efficiency of the proposed technique for the clinical diagnostic of the placental disorders [START_REF] Hua | Effects of diagnostic contrast-enhanced ultrasound on permeability of placental barrier: a primary study[END_REF], [START_REF] Guttmacher | The human placenta project: it's time for real time[END_REF], [START_REF] Thornburg | Uteroplacental circulation and fetal vascular function and development[END_REF]. To this end, the standardization of the classifier system for various ultrasound scanner systems should also be analyzed in order to allow multi-centric studies.

V. CONCLUSION

This study introduces the use of a fluid dynamic model for the analysis of DCEUS. A complete kinetic analysis was proposed in order to quantify the velocity amplitude of the bolus arrival. The estimated microbubble transport is shown to be highly correlated with in-vitro velocity of a flow system, for the delivered flow rates and for the two tested microbubble concentrations. Besides, the proposed classification criterion γ, estimated from DCEUS, outperformed existing TIC-based methods for the binary classification of ligature/non-ligature rat placentas. Our methodology opens great perspectives for the evaluation of the proposed technique for the clinical diagnostic of obstetrical disorders.

Fig. 1 :

 1 Fig. 1: Processing sequence for the quantitative evaluation of the microbubble transport amplitude occurred in the imaged tissue during a microbubble bolus.

Fig. 2 :

 2 Fig. 2: Assessment of DCEUS-based microbubble transport using the in-vitro experimental setup. (a,b): typical example of DCEUS images of the silicone tubing undergoing a flow rate of 4.8 mL/min with bubble dilutions of 1:1000 (a) and 1:2000 (b). The associated transport field, estimated using the proposed methodology, is superimposed using yellow arrows (for an easier visualization, only arrows associated to pixels with DCEUS intensity higher than 10 % of the maximal saturation value are displayed). (c,d): DCEUS-based transports flow rate estimates, as a function of the pump flow rate, independently replicated (N =3) for bubble dilutions of 1:1000 (c) and 1:2000 (d). The black dashed line represents the linear regression line; its equation and the correlation coefficient are indicated in the text box in the bottom-right of each graph.

Fig. 3 :

 3 Fig.3: Example of microbubble transport estimates from the DCEUS of two rats using the proposed approach. The first and the second rows display results associated to a control and a ligature rats, respectively (referred to as rats #11 and #7 in the scope of this study). The first column displays a B-mode image with the manually drawn mask superimposed (yellow dashed line) underlying the placenta. The 2 nd , 3 rd and 4 th columns display DCEUS images acquired 10 s, 15 s and 20 s after the determined microbubble arrival in the placenta, respectively. The corresponding transport field estimates are superimposed using yellow arrows. The maximum amplitude of the instantaneous transport is reported on the upper-right of each panel. In these results, ∆T , δt, α and f c were fixed to common typical values of 60 s, 3 s, 0.1 and f 0 /16, respectively.

Fig. 4 :

 4 Fig.4: Temporal evolution of the bolus velocity (averaged over the ROI encompassing the placenta, and weighted by the amount of grey level intensity in the DCEUS image) estimated for each control (a) and ligature (b) rat. Note that, for an easier visualization, a temporal average filter (kernel size=30) has been applied on each individual curve before display.

Fig. 5 :

 5 Fig.5: Assessment of proposed and TIC-based approaches for the classification of the two rat populations (i.e. control and ligature). (a): boxplots of overall microbubbles velocities (i.e. γ in Eq. (8)) estimated over each rat populations using the proposed approach, (b-e): distributions of PE, WiR, TTP and WiAUC estimated over the placenta for the two rat populations using the existing TIC-based approach. The median is shown by the central mark, the first and the third quartiles are reported by the edges of the box, the whiskers extend to the most extreme time points that are not considered as outliers, and outliers are individually plotted in red. Note that the term "outliers" refers here to values which benefitted from a specific display (i.e. red crosses), and not to rats which would have been removed in the statistical analysis. It corresponds to values that are more than 1.5 times the interquartile range away from the top or bottom of the plotted box. (f): ROC curve obtained using γ, (g): ROC curves obtained using PE and WiR, (h): ROC curves obtained using TTP and WiAUC.

Fig. 7 :

 7 Fig. 7: Analysis of the impact of the cutoff-frequency f c of the low-pass spatial image filtering and of the regularization parameter α of the numerical scheme on the overall results of the proposed approach. (a): bolus velocities estimated for various values of f c (∆T and δt and α and being fixed to typical values of 60 s, 3 s and 0.1, respectively), (b): bolus velocities estimated for various values of α (∆T , δt and f c and being fixed to typical values of 60 s, 3 s and f 0 /16, respectively). The significativity of the difference between bolus velocities obtained for the two rat populations (i.e. ligature and control) is reported using the AUROC.

  AUTHORS' CONTRIBUTIONS B.D.D.S. implemented the fluid dynamic model, performed image processing, analyzed the data and wrote the manuscript. A.N. designed and performed in-vitro experiments. C.J.A., A.N. and V.M. designed and performed in-vivo experiments. J.M.E. wrote the manuscript. P.A.D. and F.P. reviewed the manuscript. F.P. and A.B. reviewed the manuscript and supervised this work. All authors gave final approval for publication.

TABLE

ACKNOWLEDGEMENTS

This work was supported by grants from Inserm (AB, FP), the University F. Rabelais of Tours (AB, FP), the AERO PERINAT Association (FP), the French Society of Perinatal Medicine (FP) and the Laboratory of Excellence TRAIL ANR-10-LABX-57 (BDDS). This study has also been carried out with the financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx Bordeaux-CPU (ANR-10-IDEX-03-02).