S Nadler

O Bonnefoy
email: bonnefoy@emse.fr

J.-M Chaix

G Thomas

J.-L Gelet

Parametric study of horizontally vibrated grain packings Comparison between Discrete Element Method and experimental results

Keywords:

Numerical and experimental studies have been undertaken to analyze three parameters controlling the compaction of granular media submitted to sinusoidal horizontal vibrations. We have characterized the influence of the dimensionless acceleration , the geometry of the container and the friction coefficients on the grain velocities and on the packing densities. Above a critical acceleration , the velocities increases with . For low values of , the surface layers are compacted, whereas the bottom layers remain at their initial density. For high values of , the bottom layers get compacted, the surface layers are fluidized so that the bulk dynamic and relaxed densities decreased. In the same way, the effect of the dimensions of the container and of the friction coefficients on the packing properties has been studied for given heights of sand, acceleration and frequency. It has been shown that the influence of the two last parameters is similar to that of acceleration. The numerical results given by the Discrete Element Method appear to be in good agreement with experimental results.

Introduction

Industrial context

Vibrations-based processes are frequently used in industry to compact granular packings. For instance, in building and civil engineering works, vibrating rolls or vibratory plate compactors are used to compact and stabilize soils. Vibrations are also used in food industry; by lowering the apparent volume of granular packings, the packaging and transport costs are reduced. In electrical protection, a field that has inspired this study, vibrations or taps are used to densify and homogeneously arrange sand grains around a silver blade so that the electrical arc is extinguished properly in case of overload.

The industrial know-how, often empirically obtained in a narrow operating window, usually suffers from poor scientific understanding and thus cannot be optimized or surely extrapolated to new designs. It is hence far from easy to choose the vibrations parameters that optimize the densification process in a given situation (grains population, container geometry and dimensions . . . ).

Literature survey

From a scientific point of view, numerous works have explored granular systems submitted to periodic mechanical solicitations, either with experiments or computer simulations.

Their phenomenology appears highly governed by the intensity of the mechanical excitation, which can be quantified by , the maximum acceleration adimensionnalized by the earth gravity . Vertical and horizontal vibrations exhibit the same sensitivity to : For gently driven systems, i.e. below or close to 1, a densification process is usually observed. During a seminal experiment lead by the Chicago laboratory on spherical beads submitted to vertical taps, a logarithmic decrease of the global packing density has been reported (Nowak et al. [START_REF] Nowak | Reversibility and irreversibility in the packing of vibrated granular material[END_REF], Knight et al. [START_REF] Knight | Density relaxation in a vibrated granular material[END_REF]). In the commonly accepted explanation, the mechanical energy input induces a local and temporary stress release that allows the granular column to unjam. In this state, grains can roll, slide or move freely, thus exploring a new configuration that will be selected if the gravitational potential energy of the associated packing is lowered. The kinetic behaviour of such a process has been compared to glass relaxation [START_REF] Nowak | Reversibility and irreversibility in the packing of vibrated granular material[END_REF]. This experiment has been further analyzed in details by many authors. For instance, Philippe et al. [START_REF] Philippe | Numerical model for granular compaction under vertical tapping[END_REF][START_REF] Philippe | Granular medium under vertical tapping: Change of compaction and convection dynamics around the liftoff threshold[END_REF] have shown by Monte-Carlo simulations that a negative vertical gradient of the packing fraction was induced by vibrations.

For vigorously driven systems, i.e. above 1 and below a fuzzy threshold around 3, densification is not anymore the predominant phenomenon. In some or all regions of the granular packing, because of the high energy input level, the system is not moving continuously in the configuration space towards a denser state but is "re-initialized" at each period, loosing memory of any previous minimization of its potential energy. In this state, convection can arise.

For higher accelerations, i.e. above 3, convection becomes the predominant effect, the grain packing is fluidized and the global density decreases [START_REF] Raihane | Densification of a 3D granular bed by horizontal vibrations[END_REF][START_REF] Raihane | Experimental study of a 3D dry granular medium submitted to horizontal shaking[END_REF][START_REF] Raihane | Etude du comportement des milieux granulaires vibres horizontalement -Application au remplissage des fusibles[END_REF].

Numerous driving mechanisms have been identified to explain the convection: a sufficiently high negative vertical gradient of the packing as described above [START_REF] Nowak | Reversibility and irreversibility in the packing of vibrated granular material[END_REF][START_REF] Knight | Density relaxation in a vibrated granular material[END_REF], hysteretic phenomena like those induced by sawtooth shaped vibrating walls [START_REF] Farkas | Transitions in the horizontal transport of vertically vibrated granular layers[END_REF], collision inelasticity evidenced for instance in ratchet convection [START_REF] Van Der Meer | Spontaneous ratchet effect in a granular gas[END_REF]. In the present study, convection appears because of another mechanism, namely avalanching. This driving mechanism has been reported by Tennakoon et al. [START_REF] Tennakoon | Onset of flow in a horizontally vibrated granular bed: Convection by horizontal shearing[END_REF] in the case of horizontal vibrations and precisely detailed by Raihane et al. [START_REF] Raihane | Densification of a 3D granular bed by horizontal vibrations[END_REF][START_REF] Raihane | Experimental study of a 3D dry granular medium submitted to horizontal shaking[END_REF][START_REF] Raihane | Etude du comportement des milieux granulaires vibres horizontalement -Application au remplissage des fusibles[END_REF].

The convective patterns of vibrating granular media are amazingly various. Here, for sake of conciseness, we only report literature results concerning horizontally vibrated systems. In a numerical study on a 2D granular medium combining a hard sphere model and particle dynamics, Liffman et al. [START_REF] Liffman | Granular convection and transport due to horizontal shaking[END_REF] showed that, under horizontal vibrations, there is one single convection roll for < 1 and two counter-rotating convection rolls for > 1.2. In this second case, grains are moving upwards along the lateral walls. The experimental work conducted by M. Medved et al. [START_REF] Medved | Convection in horizontally vibrated granular material[END_REF] showed two counter-rotating convection rolls in the upper part of a horizontally vibrated bed, and even four counter-rotating convection rolls for particular situations (highly frictional walls and large values of ). They also noticed that the number of convection rolls depends on the filling height to container length ratio. In our laboratory, vertical and horizontal vibrations have been experimentally investigated by Rouèche et al. [START_REF] Roueche | des parametres de vibrations sur la rheologie d'un milieu granulaire : Application au remplissage des fusibles[END_REF] and Raihane et al. [START_REF] Raihane | Densification of a 3D granular bed by horizontal vibrations[END_REF][START_REF] Raihane | Experimental study of a 3D dry granular medium submitted to horizontal shaking[END_REF][START_REF] Raihane | Etude du comportement des milieux granulaires vibres horizontalement -Application au remplissage des fusibles[END_REF][START_REF] Raihane | Analysis of the densification of a vibrated sand packing[END_REF] respectively. Whatever the vibrations direction, it has been observed that two counter-rotating rolls appear, in the opposite way of [START_REF] Liffman | Granular convection and transport due to horizontal shaking[END_REF]: grains are moving downwards along the lateral walls and upwards in the central region. Concerning the competition between densification and convection processes, experiments carried out on an initially loose sand grains packing vibrated in a rectangular container, [START_REF] Raihane | Densification of a 3D granular bed by horizontal vibrations[END_REF][START_REF] Raihane | Experimental study of a 3D dry granular medium submitted to horizontal shaking[END_REF][START_REF] Raihane | Etude du comportement des milieux granulaires vibres horizontalement -Application au remplissage des fusibles[END_REF] have shown that the global packing density is maximized for a given value of the acceleration. The authors hypothesized the existence of a maximum as follows. With increasing acceleration, a densification front moves downwards and the thickness of the granular packing able to be densified increases. In parallel, convection gradually develops in the upper part and pre-vents the densification process in this area. These conclusions are close to those established by Richard et al. [START_REF] Richard | Slow relaxation and compaction of granular systems[END_REF] in the case of a vertically vibrated granular media. The main difference between vertical and horizontal vibration appears to be the lift up of the granular medium: in vertical solicitations, the whole sand packing unjam easily at low acceleration ( = 3 for f=50 Hz, = 2 for f=20 Hz [START_REF] Roueche | des parametres de vibrations sur la rheologie d'un milieu granulaire : Application au remplissage des fusibles[END_REF]).

After this phenomenological literature survey, let us focus on the numerical method used in this study. The Discrete Element Method (DEM) elaborated by Cundal et al. [START_REF] Cundall | Discrete numerical-model for granular assemblies[END_REF] is now a classical tool in granular media problems and gives good agreements with experiments even for complex processes. For instance, Delaney et al. [START_REF] Delaney | Validation of DEM predictions of granular flow and separation efficiency for a horizontal laboratory scale wire mesh screen[END_REF] studied industrial processes combining a granular flow, segregation and vibration effects and Cleary [START_REF] Cleary | Large scale industrial DEM modelling[END_REF] reported large scale DEM simulations. To our knowledge, few laboratories are using this method to study dense vibrated granular media. We report that Yao et al. [START_REF] Yao | DEM simulation of densification in concrete block vibration forming[END_REF] have used DEM in order to model the influence of high frequency vertical vibrations (600-1200 Hz) on the densification of granular media. Matchet et al. [START_REF] Matchett | Vibrating powder beds: a comparison of experimental and distinct element method simulated data[END_REF] have compared DEM and experimental results concerning the energy loss by vertically vibrated granular media, and have shown that this method is partially satisfying since the calculation of the applied force was not successful.

Close to DEM, another method called "Molecular Dynamics", is also used for studying various phenomena. In this method, Coulomb friction is differently implemented and the rotation of particles is neglected [START_REF] Gallas | Convection cells in vibrating granular media[END_REF]. Gallas et al. have employed this method to study the influence of vertical vibrations on a granular medium. It has been reported that the friction, the geometry, the presence of an obstacle [START_REF] Gallas | Granular media on a vibrating plate -a molecular dynamics simulation[END_REF] and the dissipation coefficient [START_REF] Gallas | Convection cells in vibrating granular media[END_REF] have a large influence on the behaviour of the medium. The segregation phenomenon can also be reproduced by this method [START_REF] Gallas | Molecular dynamics simulation of size segregation in three dimensions[END_REF]: for a big sphere immersed in a medium, no segregation occurs for the lower frequencies (f≤20 Hz), segregation has been evidenced for f≈40 Hz, and a more complex behaviour for f≥50Hz (up and down displacements of the big sphere). Saluena et al. have shown that the dissipation coefficient also has a major influence on the shape and intensity of the convection rolls in horizontally vibrated granular media [START_REF] Saluena | Convection in horizontally shaken granular material[END_REF]: the intensity of the convection rolls, their number and direction can change when varying these parameters. We must notice that numerical results are not systematically compared with experimental results.

As a consequence, we think that the study of a vibrated dense granular medium combining experimental and DEM modeling, is an interesting contribution to the literature.

Issues

In the previous section, we have mentioned that presents a great influence on the vibrated granular system. This influence will be studied in this paper both experimentally and by simulation. Other parameters as the frequency f, the amplitude of the vibration A, the recipient geometry, The influence of these three parameters is presented in this same paper in order to show the similarities of the response when varying one of them. The studies concerning the influence of and of the container size combine laboratory experiments and numerical simulations. The influence of the friction coefficients is mainly studied with DEM.

Materials and methods

Experimental set-up and materials

Our experimental conditions are inspired by the industrial manufacturing process of electrical fuses, which are made out of a ceramic recipient containing a silver blade with sand grains in-between. Here, we use a rectangular recipient driven horizontally by an electromagnetic vibrator. The delivered force is continuously adjusted so that the motion is sinusoidal, with a frequency f=50 Hz. The forcing intensity is described by the peak acceleration. More conveniently, we use the dimensionless acceleration = (2 ) 2 / , where is the amplitude of oscillations and the earth gravity, the values of ranging between 0 and 5. The container is made out of 10 mm-thick transparent PMMA walls. The dimensions of the container are denoted , and respectively along the vibration axis x, the horizontal transverse axis y and the vertical axis z. The vertical walls are labeled north/south, east/west(fig. 1). The granular medium is made of rounded sand grains, composed by more than 99.2% of quartz (controlled by XRD). The number size distribution of the grains, as measured by Laser granulometry, is very well approximated by a skew normal distribution with an average diameter < >= 370 m, a standard deviation = 216 m and a skewness of 3.6. The wall-grain friction coefficient has been measured by shear tests performed with a powder rheometer (FT4, Freeman Technology); its value is / = 0.3. The heap angle of repose of the granular medium has been experimentally determined and is equal to = 30 ∘ . As it will be established in section 3.4, it corresponds to a numerical grain-grain friction coefficient / = 1. Table 1 summarizes these data, as well as other parameters useful for numerical computations like shear ratio = 2(1+ ) (E, : Young's modulus and Poisson ratio, with standard values for silica sand grains [START_REF] Sawicki | Elastic moduli of noncohesive particulate materials[END_REF]), numerical and viscous damping D and D [START_REF] Ollagnier | Contact particles calibration in order to use discrete element code[END_REF].

The sample preparation protocol consists in pouring the granular medium in the container with the help of a funnel of 2 cm diameter and then flatten the free surface by suction (bed height of

). Such a protocol generates initial loose packings with a density of 59.5%±-0.5%.

Experimental data acquisition tools

Two main data are used to characterize the state of the vibrated granular medium: the global packing density and the velocity field of the grains moving along the transparent container walls. To accurately measure the global density, we use a confocal optical pen, whose 200 m diameter light beam scans the free surface of the packing. By integration, we then compute the apparent volume of the granular medium and its global density (averaged density in the granular media). The accuracy of the global density is estimated at 0.5 %, both in dynamic (with vibration) and static conditions. The density measurement device is detailed in [START_REF] Raihane | Analysis of the densification of a vibrated sand packing[END_REF], where it has been extensively used to analyse the densification of a vibrated sand packing. The 3D grains trajectories can be experimentally measured with sophisticated tools such as Positron Emission Particle Tracking (see Hoomans et al. [START_REF] Hoomans | Experimental validation of granular dynamics simulations of gas-fluidised beds with homogenous in-flow conditions using positron emission particle tracking[END_REF]). Here, we have chosen a simpler method, namely the Particle Image Velocimetry. A fast camera takes pictures of the grains at the immediate vicinity of a given face of the recipient at a rate of 50 frames per second, that is to say in phase with the forced oscillations. Then, with a specially designed method, we analyse a set of 50 pairs of images to produce the velocity field. On the north/south walls, this 2D velocity field is "summarized" in a 1D velocity profile, where, for each altitude , we compute the vertical component of the mean grain velocity (the horizontal component is comparatively, negligible).

Modeling method and data analysis approach

The Discrete Element Method is used in this work with the commercial software PFC3D from ITASCA. A sinusoidal horizontal motion is imposed to a rectangular recipient containing a collection of grains. Each grain is modeled by a soft sphere and its position is computed at each time step on the basis of the gravitation force and contact forces that apply on it. For the contact force, we use the standard Hertz model with a friction coefficient and a stiffness coefficient . Results of modeling works found in literature are usually difficult to reproduce because the set of parameters is not extensively indicated. Table 1 presents the values chosen in our study. We point out that, numerical simulations have been carried out using a truncated gaussian Particle Size Distribution almost proportional to the real one. The polydispersity offers the advantage of preventing crystallization [START_REF] Pouliquen | Crystallization of non-Brownian spheres under horizontal shaking[END_REF]. When not specified, the mean diameter is twice that of the experimental one to keep the particles number below 80,000 and hence the calculation time under 150 hours. We assume that the difference of shape does not affect the qualitative behaviour of the system [START_REF] Ribiere | Experimental compaction of anisotropic granular media[END_REF]. The initial state is prepared by sedimentation of spheres randomly created in space. After the force balance is reached, the density is measured. It is rather loose with a value of 59.5%±0.5%.

In order to represent a wide system with the minimum number of particles (to get a reasonable calculation time), we have used periodic boundary conditions along the horizontal y-axis that is orthogonal to the vibration x-axis. The influence of the size of the periodic side on velocities is shown on fig. 2. The minimum length of the unit cell required to avoid finite-size effects appears to be 10 grain diameters. For larger values, the velocity profile along the north face becomes independent of the unit cell size: the container behaves as if the transverse dimension were infinite, i.e. without the influence of the side walls. Let us remark that the maximal velocity on north face can be considered as one of the most sensitive parameter to the periodic size and we assume that if this scalar remains the same, the global behaviour of the system remains unchanged. The reader can notice that the minimum unit cell size depends on the numerical experiment. For instance, [START_REF] Anand | Predicting discharge dynamics from a rectangular hopper using the discrete element (DEM)[END_REF] found that a unit cell of 2 grains size is enough for describing the discharge of a rectangular hopper.

Moreover, we have transformed the natively used Lagrangian description (grains are followed individually) into an Eulerian description (continuum described by density-, velocity-, stress-, ... fields) in order to better visualize the modeling results and easily compare them with the experimental ones, especially with the velocity fields. The averaging procedures are standard ones. It should be notice that, in order to smooth the fluctuations due to the chaotic nature of the medium behaviour, the data have been averaged over a series of 10 oscillation periods. The velocity of a grain for one oscillation period has been calculated by considering its displacement during one cycle of vibration and dividing it by the period of the oscillation. Afterwards, the velocities have been averaged in space (eulerian procedure). Observation of the counter-rotating convection rolls. =3 and f=50 Hz. In these conditions, the rise of the grains is not located at the wall surface for this acceleration, but in the bulk of the granular packing. =500 m

Results

We have chosen to study the influence of the dimensionless acceleration , the length of the vibrated box, and the friction coefficients upon the global packing fraction and the velocities of the grains along the transparent walls of the container. In section 3.1, we present results where the particles diameters in the simulation are the same as in the experiments. In the following sections, the diameter in simulation is twice the experimental one.

Convective pattern and simulation calibration

For f=50 Hz and an acceleration above the critical value = 0.8 (under which no bulk movement is noticed), the rheological behaviour of the granular medium can be mainly described by two counter-rotating convection rolls in the upper region as described in fig. 1, fig. 3 and fig. 4. This convective motion is driven by the avalanching mechanism: the vibrations imposed along the north-south axis create an intermittent gap between the granular medium and the north or south wall alternatively. Medved et al. [START_REF] Medved | Connections between response modes in a horizontally driven granular material[END_REF] have studied the existence and evolution of this gap with . By gravity, the grains can fall down into this gap. Then, they penetrate into the bulk, are pushed horizontally by the walls and go upwards in the central part of the container. When they arrive at the top of the pile, the grains may fall down again and the loop is closed [START_REF] Raihane | Etude du comportement des milieux granulaires vibres horizontalement -Application au remplissage des fusibles[END_REF].

The validity of the parameters chosen in the model has been checked by two experiments and simulations where most geometrical parameters (size of the grains, width of the box and height of sand, . . . ) were the same. The only difference is the periodic side created in the simulation. Figure 3 and fig. 4 show respectively the experimental and numerical velocity field averaged over 10 cycles of vibration at = 3, f=50 Hz, =20 mm, =20 mm, and =40 mm in experiment and =10 mm periodic in modeling. The reader will notice that the rise of the grains is evidenced in the simulation (fig. 3) but not in the experiment (fig. 4). Effectively, the grains seem motionless at the surface wall and they go up in the core of the medium. In addition, the extent of the region at rest is larger in DEM than in experiments and is attributed to the difference of grain shape (perfect spheres vs. rounded sand grains) and boundary conditions (physical walls vs. periodic boundary conditions).Apart from these two differences, with identical diameters, the convective patterns and the order of magnitude of the velocity fields show a good agreement. Figure 5 shows that, for the same diameter, the numerical and experimental velocity profiles along the north face are qualitatively and quantitatively similar (error≈15 % on the maximal velocity on north side). We can point out that for low z-values, i.e. at greater depth, the grain velocity is larger in the numerical simulations than in the experiments. We attribute this mobility increase to the spherical shape of the modeled grain that allows them to move more easily by rotation or sliding than it would be with more polyedric grains. It also shows that in the given conditions, multiplying the diameter by two implies a division of the maximal velocity by 7.

In the following sections, to avoid excessively long calculation time, simulations are performed with grains twice as large as the sand grains used in experiments. Therefore, we draw the attention of the reader on the fact that the numerical results exposed in the next sections will give a more qualitative than quantitative description of the granular medium behaviour.

Influence of acceleration on packing properties

Density

Experiments

On exactly the same system, Raihane et al. [START_REF] Raihane | Analysis of the densification of a vibrated sand packing[END_REF] have measured the local relaxed density of the packing (fig. 6). For low accelerations, it has been shown that the bottom layers remained at the initial density, whereas the upper lay- ers undergo a densification. For higher accelerations, the bottom layers get compacted, but the density decreases in the upper layers. We can explain this latter drop by the fact that the upper grains acquire a significant velocity. In this case, the free volume is enhanced and when vibrations are stopped, the density keeps low.

Simulations

An assembly of 32000 spheres of density = 59.5% is created. This relaxed state can be characterized by the bulk density map (fig. 7). The behavior of the granular packing during the vibrations depends on the acceleration.

For low accelerations ( < 1.5), the density at the bottom of the container remains equal to the initial density. This occurs because the Coulomb criterion has not been reached in this region, which indicates that the compressive stress is high compared to the shear stress.

For high accelerations ( > 3.0), the shear stress is high enough and the Coulomb criterion is fulfilled, even in the deeper layers. Therefore, the grains can move and the packing can densify even at the bottom of the container. In the upper layers, the granular medium becomes fluid and the density decreases.

The local relaxed density increase ranges between 0 % in the jammed bottom layers and 3 % in the fluidized top layers of the packing. The mapping of the density in the dynamic state is given in fig. 7. Experimentally and numerically, we can observe a further small gain in densification ( ≈ 0.5%), when vibrations are stopped. Simulation results on the density evolution vs. the relative acceleration are in agreement with Raihane [START_REF] Raihane | Analysis of the densification of a vibrated sand packing[END_REF] experimental data.

Velocity of particles at north and south walls

Experiments and simulations have been undertaken in order to obtain the evolution of the maximal velocity along the north face with (fig. 8). Both experimentally and numerically, for 0< <5, the velocity of the grains increases linearly with the acceleration. We notice a good qualitative agreement between modeling and experimental results (fig. 9 and fig. 10). The quantitative differences between the experimental and computed velocities can be attributed mainly to the size of the grain (double in the model), and secondly to the shape of the grain and the choice of a periodic side in simulations.

Influence of the container size on the packing properties

Density

The granular medium behaviour obviously depends on many parameters. Hereafter, we focus on the influence of the size of the container. Indeed, we have to be cautious when considering the influence of on the bulk dynamic density of the system (fig. 6, fig. 7): these results are valid only for a given container dimension. Numerical simulations and experiments have been performed to evaluate the influence of the length (measured along the vibration axis) upon the density (fig. 11, fig. 12). Numerical simulations (fig. 12) show that, for small values of (20 mm), only the upper layers are compacted, whereas for greater values (starting from =60 mm) all the packing is compacted. This behaviour can be explained by a larger constraint on the displacement of the grains in small containers. Comparing fig. 7 and fig. 12, we can see that increasing has got a similar effect to imposing a larger acceleration : the mobility of grains increases and it enables to compact the bottom layers. Experimentally, for =60 mm, we verify that the dynamic or relaxed density increases with the size of the container (fig. 11). The modeling results suggest that the compaction of the bottom layers occurs only if is long enough. This behaviour explains the experimental evolution of the global density with .

Velocity of particles at north and south walls

For a given bed height ( =60 mm), experiments show (fig. 13) that the velocity along the north face and the fluidized zone thickness increases with the size of the con- tainer in the range =20-60 mm (Flemmer et al. have anticipated this type of behaviour in [START_REF] Flemmer | Coherence of a packed bed under lateral oscillation[END_REF]). For small containers ( =20 mm), no motion is observed at the bottom of the container.

In fig. 14, we have plotted the velocity profiles along the north face obtained in simulation. The grain assembly exhibits qualitatively the same behaviour in experimental and modeled systems: the velocity increases with until =60 mm.

Experimentally and numerically, the velocities can proportionally increase with (for small values of ). To explain this phenomenon, we suggest that the steric constraint is weaker when the size of the box increases and the grains can move faster and in the deeper layers. Furthermore, we can state that the geometry and especially the length of the box have an similar influence to that of the acceleration (fig. 9 and fig. 10).

Influence of the friction coefficients on the packing properties

In this section we use simulations to analyze the influence of the friction coefficients / and / on the packing behaviour.

Heap angle

For a given powder, the heap angle depends on many parameters such as the attractive or repulsive forces acting at the friction contact forces due to surface roughness of each individual grains and the contact morphology and associated entanglement. When a specific "real" powder is modeled by perfectly spherical particles, the coefficient of friction / must be chosen carefully. One robust way to do so is to adjust this parameter so that the heap angle of the model packing is the same as the experimental one. Experimentally, the heap angle of

Density

Figure 16 shows that the global dynamic density gently decreases with the grain-grain friction coefficient / . We can also see that the density profile is similar for / = 0.1 and / = 1. Let us just notice a global increase of density of approximately 1 % for lower frictional grains (fig. 17). This observations correlates well with Vandewalle et al. [START_REF] Vandewalle | The influence of grain shape, friction and cohesion on granular compaction dynamics[END_REF], who experimentally found a difference of 1 % on the relaxed density after compaction between smooth and rough particles. The normalized volume fraction is defined as ′ = ( -0 )/( -0 ) where is the density, 0 the minimal density, and the maximal density. Van-

Velocity of particle at north and south walls

We have observed that numerical velocities strongly depend on the values of / and / . In fig. 18, one can see that velocities increase sharply with both friction coefficients. An unexpected result is obtained when grain-grain friction coefficient is very low: the velocity becomes negative for / <0.1 (fig. 18) and the convection rolls are reversed (fig. 19). In this case, the grains climb along the north and south walls and fall in the middle of the granular medium. Liffman et al. [START_REF] Liffman | Granular convection and transport due to horizontal shaking[END_REF] have also shown this type of behaviour, which can be explained by the lack of grain/grain friction coefficient in their model. We can remind that this result is essentially numerical and cannot be reproduced experimentally: at our knowledge, there is no real granular medium with null friction coefficient.

Conclusion and prospect

In this work, we have shown the strong influence of , the length of the container, and the friction coefficients on the density and on the velocity fields of a granular medium submitted to horizontal sinusoidal vibrations. The velocity always increases, when any of these parameters is increased. For high values of and , a densification of the deepest layers occurs, whereas the surface is fluidized. For low value of and , densification of surface layers may occur, whereas the bottom layers remain unjammed. The correlation between experimental and modeling results is qualitatively good, but we have to carefully interpret the numerical results. Particularly, it has been shown that the particle size has a large impact on the behaviour of the system. We have also underlined that the friction coefficient may have a major influence on the rheology of the granular medium. Increasing numerically the friction coefficient leads to convection rolls flowing in opposite direction. It is probably not possible to establish such experimental results. However, modeling results suggest that the roughness influences the velocities of the grains. As the influence of the three parameters appears to be qualitatively the same in simulations and experiments, an experimental study of the influence of the roughness on the velocity of the grains would be interesting to verify if the velocities of particles implied in the rolls are affected by friction coefficient. The influence of the grain size has also been suggested by the simulation: a specific study on this size effect would also complete our knowledge. Further studies on the influence of the particle shape are also planned.

Fig. 1 .

 1 Fig. 1. Counter-rotating convection rolls observed experimentally when f=50 Hz and > 1.5

Fig. 2 .

 2 Fig. 2. Influence of the size of the periodic cell on the velocity of spheres along the north face and focus on the maximum velocities, f=50 Hz, = 3, =40 mm

Fig. 3 .

 3 Fig. 3. East side view. Velocity field measured experimentally by PIV technique and averaged over 10 cycles of vibrations.Observation of the counter-rotating convection rolls. =3 and f=50 Hz. In these conditions, the rise of the grains is not located at the wall surface for this acceleration, but in the bulk of the granular packing. =500 m

Fig. 4 .

 4 Fig. 4. East side view. Numerical velocity field averaged over 10 cycles of vibrations. Observation of the counter rotating convection rolls. =3 and f=50 Hz. = = 500 . Due to periodic boundary conditions, the surface velocity fields coincide with the bulk one.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Experimental and numerical velocity profiles for identical sand height and box length at =3, f=50 Hz

Fig. 7 .

 7 Fig. 7. Computed relaxed density [%] as a function of and for =50 Hz

Fig. 8 .

 8 Fig. 8. Maximum velocity on north face vs. . Experimental velocities averaged on 50 vibration cycles (standard deviation ≈ 15%). Modeled maximum velocity averaged on vibration cycles 90-100 (standard deviation ≈ 30%). f=50 Hz, =40 mm, H=40 mm

Fig. 9 .

 9 Fig. 9. Experimental velocity profile dependence upon for f=50 Hz, =40 mm.

Fig. 10 .

 10 Fig. 10. Simulated velocity profile dependence upon for f=50 Hz, =40 mm.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Density increase due to vibrations as a function of , for =60 mm : experimental and numerical results. Typical error ± 0.5%, = 3, f=50 Hz

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Experimental velocity profile dependence upon for = 60 mm, = 3, f=50 Hz

Fig. 15 .

 15 Fig. 15. Modeling result : heap angle as a function of the grain/grain friction coefficient imposed in the DEM.

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Computed evolution of the dynamic global density vs. grain-grain friction/

Fig. 18 ./.

 18 Fig. 18. Numerical study of the dependence of the velocities on/

Fig. 19 .

 19 Fig. 19. Calculated velocity field indicating a reverse convection roll for/

Table 1 .

 1 Default main parameters

	Parameters			Experimental Model
	Mean grain diameter		[ m]	450	900 or 450 (1 exp)
	Density		[kg m -3 ]	2700	2700
	Shear ratio		[MPa]	40	40⋅10 3
	Poisson's ratio		[-]	0.18-0.27	0.25
	Friction coeff.	/	[-]	0.3	0.3
	Friction coeff.		[-]	?	1 (see sect. 3.4)
	Container size		[mm]	40	40
	Container size		[mm]	40	periodic
	Restitution coeff.		[-]	0.7	0.7
	Viscous damping coeff.		[N.s.m -1 ] -	0.1
	Numerical damping coeff.		[-]	0	0.1
	Number of grains		[-]	≈ 10 7	≈ 5.10 4

/


Ahmad Almhdie
email: ahmad.almhdie@univ-orleans.fr

Patricia Lopes-Pereira

Sandra Même

Caroline Colombier

Véronique Brault

Frédéric Szeremeta

Bich-Thuy Doan

Roger Lédée

Rachid Harba

Yann Hérault

Jean-Claude Beloeil

Christophe Léger

Chan-Vese based method to segment mouse brain MRI images: application to cerebral malformation analysis in trisomy 21

In this paper, a semi automatic active contour method based on Chan-Vese model is proposed for the segmentation of mouse brain MR images. First, a 2 ½ D strategy is applied on the axial images to segment the 3D volume of interest. The method takes into account the special shape of the object to segment. Moreover, the user defines the limits where to search these contours and also provides an initial contour. This semi automatic method makes that human intervention is limited and the tedious manual handling is greatly reduced. Results have shown that the brain volumes estimated by the method are identical to expert manually estimated volumes. Last but not least, the new method was used in the analysis of the cerebral malformations linked to Trisomy 21: no significant difference of the brain volumes between Trisomy 21 mice and the control ones were found.

1.

INTRODUCTION

This work is a part of a project granted by the French Region Centre which associates 3 laboratories of the city of Orléans, France. The project aims to a better understanding of the morphological malformations of the brain associated to Trisomy 21 or Down syndrome. Characterization of morphological malformations associated to trisomic mouse models for particular chromosome segments will help to link genes from these segments to specific dysmorphologies. It would also help to look after the consequences of this disease, and to improve the quality of life of the patients. For this purpose, models of trisomic animals are usually of interest. Among various animals that can fit the study, mouse are often chosen since 80% of the genome is common with that of the Human. In addition, mice are easy to breed. For such studies, 3D images of the mouse brain must be analysed. One of the main strengths of MRI compared to other imaging techniques (optical, PET, CT, ultrasound,…) is the ability to obtain not only anatomical images in 3 dimensions but also functional information with the implementation of new methods of imaging (fMRI, diffusion, perfusion, imaging of chemical shift, imaging with contrast agents...). For this reason, MRI images will be used in our study. 3D mouse brain MRI images must be segmented in order to extract the brain itself. For the moment, most of biologists working in this domain perform a manual segmentation of these 3D data. For that reason, an automatic or semi automat-ic method is needed. Consequently, human intervention will be limited and the tedious manual handling will be greatly reduced. Even if a priori information is taken into account, and even if the brain is surrounded by the cerebrospinal fluid, the automatic segmentation of 3D mouse brain MRI images is a difficult issue since the contrast between different objects in the brain is low. A semi automatic technique that could consider information given by the user could be more efficient and less complicated to develop. A good review of current medical image segmentation methods can be found in [START_REF] Pham | Brain MRI Segmentation with Multiphase Minimal Partitioning: A Comparative Study[END_REF]. For the segmentation process itself, a Chan-Vese based approach of active contour [START_REF] Chan | Active contours without edges[END_REF], which is efficient in the case when image edges are not well defined, seems to be well adapted for our application. In this work, we propose a semi automatic Chan-Vese active contour based method to segment mouse brain MRI images. The method includes some a priori information and takes into account user knowledge. This paper is organized as follow. First, the material will be presented, i.e., the production of trisomic mouse models, and the acquisition of related 3D MRI brain images. Section 3 will describe the method. Next sections present the results and some perspectives. Concluding remarks are drawn in the final section.

MATERIAL

Trisomic mouse models: Different mouse models have been created to elucidate genetic pathways implicated in specific clinical aspect of the trisomy 21 pathology (Figure 1). In addition to the Tc1 model which corresponds to a mouse line carrying an additional quasi-complete copy of HSA21 [START_REF] O'doherty | An Aneuploid Mouse Strain Carrying Human Chromosome 21 with Down Syndrome Phenotypes[END_REF], four other models are trisomies of counterparts regions of chromosome 21 carried by the mouse chromosome 16 (MMU16; MMU for Mus musculus): Ts1Yu, Ts65Dn, Ts1Cje, Ts1Rhr. Two other models, named Ts1Yah and Ts2Yah, have recently been created by the team of Yann Hérault at the IEM laboratory in Orleans. They are trisomic respectively for the part of the chromosome 21 counterpart region on MMU16 absent from the Ts65Dn model and for the chromosome 21 counterpart region in MMU17. Analyses of these different models have revealed alterations of the cerebral and cranial development. Skull and brain defects were mainly analyzed in the Tc1, Ts65Dn, Ts1Cje and Ts1Rhr models and are discussed in [START_REF] Aldridge | Differential effects of trisomy on brain shape and volume in related aneuploid mouse models[END_REF][START_REF] Olson | A chromosome 21 region does not caus specific Down syn-drome phenotypes[END_REF][START_REF] Olson | Down syndrome mouse models Ts65Dn, Ts1Cje, and Ms1Cje/Ts65Dn exhibit variable severity of cerebellar phenotypes[END_REF]. It was shown that Ts65Dn mice have a reduced number of granule cells in the cerebellum and the hippocampus. When considering cerebellum, this reduction is accompanied by a decrease in the cerebellum size which is visible during postnatal development. No study has been made on the recently published Ts1Yu model [START_REF] Li | Duplication of the entire 22.9 Mb human chromosome 21 syntenic region on mouse chromosome 16 causes cardiovascular and gastro-intestinal abnormalities[END_REF]. The main conclusion from these studies was that the brain and skull defects observed in these models cannot be directly associated to one gene or to a chromosomal segment, but suggests that interactions between several genes throughout the chromosome 21 are responsible for the dysmorphologies seen in Down syndrome patients. Further studies are yet needed with new trisomic models for chromosomal segments that have not yet been studied such as in the Ts1Yah and Ts2Yah models in order to determine the importance of each trisomic region to the complete phenotype. In order to further explore the role of each trisomic segment on brain patterning, we will analyse 3D reconstructions of MRI images of brains from Ts1Yah, Ts2Yah and control mice. 5 individuals of each group will be produced by taking adult animals of similar age, sex and genetic background. A contrast agent will be injected to the mice before cerebral MRI imaging, following the protocol already approved by the ethics committee (Decision CL2006-027 granted to Yann Hérault) for preliminary experiments which were conducted in order to define the best experimental conditions for the acquisition of images. Cerebral MRI imagery of the small animal: One of the main strengths of MRI compared to other imaging techniques (optical, PET, CT, ultrasound, etc.…) is the ability to obtain not only anatomical images in 3 dimensions but also functional information with the implementation of new methods of imaging (fMRI, diffusion, perfusion, imaging of chemical shift, imaging with contrast agents etc. ...). However, one of the major limitations of MRI is its low sensitivity.

HSA21

To overcome this defect, it is possible today to use increasingly intense static magnetic fields. Nevertheless, the increase in the magnetic field increases the sensitivity to the effects of magnetic susceptibility and reduces the natural contrast between different tissues. Thus, one of the research objectives in MRI today is to propose new methodological, technological and instrumental developments to improve the contrast, sharpness, speed and spatial resolution at high field. So, the search for new contrast agents of contrast more adapted to the intense magnetic fields will allow the IRM to have a better specificity. For example, the Manganese ion Mn 2+ has paramagnetic properties leading to an enhancement of the MRI signal [START_REF] Aoki | In vivo detection of neuroarchitecture in the rodent brain using manganese-enhanced MRI[END_REF][START_REF] Lee | Manganese-enhanced magnetic resonance imaging of mouse brain after systemic administration of MnCl2: dosedependent and temporal evolution of T1 contrast[END_REF][START_REF] Pinter | Neuroanatomy of Down's syndrome: a highresolution MRI study[END_REF]. In this work, the MR experiments were performed using Manganese Mn2+ as contrast agent within the framework of a cerebral application in order to carry out the quantitative and qualitative follow-up of the impact of the trisomy of certain regions of the chromosome in small animal models of the disease. Several experiments were conducted to obtain images with high spatial resolution (118.288.5118.2 m 3 /voxel) with many details of the neuro-architecture, thanks to the power of the magnetic field and imaging gradient available on the used imaging system, first of its kind in France. MRI experiments were carried out on a spectrometer 9.4T/21 USR Biospec (Bruker, Wissembourg, France) using a stereotaxic setup. The mice were anaesthetized with isoflurane (1.5-2.5%) in an O2/N2O 1:1 mixture. Breath and temperature physiological parameters were monitored throughout the experiment. 3D FLASH experiments were recorded using, TR= 150 ms, a = 60°, duration 2h.

METHOD

Remember that we aim at developing a semi automatic segmentation for 3D mouse brain MR images. We will also take into account the a priori information relative to the specificity of the treated images, and will use some information given by the user. This strategy is explained in the following paragraphs.

Due to the elongated shape of the mouse brain, we have chosen a 2 ½ D contour detection based on axial slices. On axial slices, the brain mostly appears as a round shape that will be easier to detect (Figure 3c), compared to the other two sections (Figure 3b and3d). First, the user determines the starting and ending point of the brain. The search of contours in the axial slices is limited in the region where the brain is present.

After that step, the contours of the brain in the initial axial slice are estimated using a Chan-Vese based approach of active contour segmentation [START_REF] Pham | Brain MRI Segmentation with Multiphase Minimal Partitioning: A Comparative Study[END_REF]. Exhaustive tests, not reported here, have confirmed the superiority of the Chan-Vese approach compared to other methods for our application. Since Chan-Vese method requires an initial contour, an initial ellipse that best fits the desired organ is positioned by the user. These two tasks (determining the starting and the ending slices, and positioning the initial ellipse) are the only manual interventions of the user. By this way, human intervention is limited, which reduces the tedious manual handling.

Even if the brain is surrounded by the cerebrospinal fluid that makes this segmentation not too complex, some difficulties arise when using the Chan-Vese method. The contours given by the Chan-Vese method are not only the correct ones, but others are also detected. To overcome this difficulty we propose the following method. Since we know that the correct contour is a long one, a Mathematical Morphology method operation (a closing) is used on binarized level set images. As a result, small contours disappear. We also know that the right contour is close to the initial contour. Therefore, we only retain the contour which is the closest to the initial one. At this point, the processing of the initial slice is over.

Next slice is then considered. Using as initialization the contour found in the previous image, the contour of the next image is rapidly detected. This process is iterated until the limit given by the user is reached. The precise strategy is detailed and illustrated in the following.

Initialisation: we shall follow the Chan-Vese model. Let u: be the original image of an axial slice to be segmented, and φ be its corresponding level set function, where: φ > 0 inside the contour, φ < 0 outside the contour, φ = 0 on the contour.

Figure 4 shows a typical example of the level set function expressed in terms of a signed distance function to the original contour initiated by an ellipse that best fits the desired organ. This ellipse is given by the user. Contour detection: the Chan-Vese method is used to detect the contour in the axial slice. This step will also produce a level set image. This level set image is binarized with a threshold equal to zero as illustrated in Figure 5.

Contour selection: since the segmentation is based on level set formulation, multiple contours could result. Hence, contour selection is preceded by firstly applying a morphological closing to the binarised level set image to separate attached contours. The structuring element is 5×5 pixels. As show in Figure 6, the closest contour to the initial ellipse is then chosen.

This process is iterated until the limit fixed by the user is reached. The user only defines the original ellipse and the range of search of the various contours.

An example is presented in Figure 7 for the 3 slices of Figure 3, whereas Figure 8 presents the volume of the detected brain structure.

RESULTS

When a user manually determines the brain volume, the processing time is of about 4 hours. Our semi automatic process reduces the time to about 15 minutes. This is a quite interesting time reduction. In the following, two additional results are given. We first compare the estimated brain volume by our semi automatic method to manually estimated volumes measured by an expert using the OsiriX ® software. In this study, 15 mice were considered: 5 controls, 5 TS2Yah and 5 TS2Yah trisomic models (see section 2 for details). Results are shown in Table 1, Table 2 A Student t test shows that there is no statistical significant difference between the measured volumes by an expert and our semi-automatic measurement. This demonstrates that the proposed method is of high quality. Moreover, a Student t test shows that there is no statistical significant difference in brain volume between control mice and trisomic ones. Although this result might seem surprising in a first step, it is already known in the literature of this domain [START_REF] Olson | A chromosome 21 region does not caus specific Down syn-drome phenotypes[END_REF].

FUTURE WORK

It is known that differences between trisomic and control mice brain mainly occurs in the cerebellums and the hippocampus. The cerebellum is responsible of the equilibrium and of the walking process. The hippocampus rules the learning and memory procedure. It is known that there is a significant difference between a trisomic subject and a control one relatively to these various tasks. It will of great interest to have a good knowledge of the relations between the volumes of the cerebellums and the hippocampus and these functions.

However, the cerebellums and the hippocampus are difficult to segment due to their complicated shape. In addition, limits of these objects included in the brain are unclear. The high quality results of segmentation of the brain that we obtained could be used to perform such a task. A possible strategy is described in the following for the segmentation of the cerebellum (it would be identical for the hippocampus).

We start from a reference brain and a reference cerebellum of a given mouse which are manually segmented by an expert. On the mouse under consideration, we first could use a rigid registration process already developed in our institution (only translations, rotations and homothetic deformations are allowed) [START_REF] Almhdie | 3D Registration Using a New Implementation of the ICP Algorithm Based on A Comprehensive Lookup Matrix: Application to Medical Imaging[END_REF]. The reference brain is first registered with the brain segmented using the semi automatic process presented in this paper (first line of Figure 9). The measured rigid transformation between the current brain and reference brain is then applied on the reference cerebellum (second line of Figure 9). The resulting registered cerebellum is used for initializing the segmentation of the true cerebellum. Such a strategy is under active development.

Figure 9 -Reference data are presented in blue color, whereas segmented brain is in purple. Top of the figure, the reference brain is registered with the brain under consideration. Bottom of figure, the rigid registration is applied to the reference cerebellum. This will be the initialisation to the segmentation process of this object.

CONCLUSION

In this work, we have proposed a semi automatic Chan-Vese active contour based method to segment mouse brain MRI images. This method associates the Chan-Vese approach to a priori information and also uses knowledge given by the user. Results first show that the brain volumes estimated by our method were identical to expert manually estimated volumes. Second, no significant difference of the brain volumes between Trisomy 21 mice and control ones were found. We are currently investigating the use of a previously proposed registration method for the segmentation of the cerebellums and the hippocampus.

Figure 1 -

 1 Figure 1 -Mouse models of the Trisomy 21.

Figure 1

 1 Figure 1 gives a schematic representation of the mouse models available. Human chromosome 21 (HSA21 for Homo Sapiens) is shown on top with a bold line with the graph of the Tc1 mouse line. The corresponding homologous regions in the mouse region are indicated labelled as MMU16, MMU17 and MMU10 with genes located at the borders. Respective position of the mouse models, Ts1Yu, Ts65Dn, Ts1Cje and Ts1Rhr are indicated with a black line. New models, Ts1Yah and Ts2Yah are shown with a dashed line.In order to further explore the role of each trisomic segment on brain patterning, we will analyse 3D reconstructions of MRI images of brains from Ts1Yah, Ts2Yah and control mice. 5 individuals of each group will be produced by taking adult animals of similar age, sex and genetic background. A contrast agent will be injected to the mice before cerebral MRI imaging, following the protocol already approved by the ethics committee (Decision CL2006-027 granted to Yann Hérault) for preliminary experiments which were conducted in order to define the best experimental conditions for the acquisition of images.

Figure 2 -

 2 Figure 2 -Considered mouse brain models

Figure 3 Figure 3 -

 33 Figure 3 -An example of a resulting 3D image and its related slices

Figure 4 -Figure 5 -Figure 6 -

 456 Figure 4 -an example of initial contour (a) and the corresponding level set image (b).

Figure 7 -

 7 Figure 7 -three orthogonal sections, after final segmentation.

Figure 8 -

 8 Figure 8 -3D visualisation of the segmented brain.

Table 2 -

 2 and Table3, respectively. volume estimation of brains of 5 trisomic TS2Yah mice

	Data	Vo lume (mm 3 )
	Reference	Automatic	Manuel
	13137	467.96	428.4
	13577	463.17	452.1
	13654	457.41	440.3
	13580	439.20	460.1
	13653	415.55	454.1
	Table 1 -volume estimation of brains of 5 control mice
	Data	Volume (mm 3 )
	Reference	Automatic	Manuel
	13136	446.49	457.6
	13579	467.87	451.6
	13578	452.60	447.9
	13695	380.12	413.8
	13582	415.63	447.8
	Data	Volume (mm 3 )
	Reference	Automatic	Manuel
	13130	437.74	417.7
	13544	483.10	478.6
	13549	475.71	456.6
	13548	414.41	433.4
	13644	454.52	460.0

Table 3 -

 3 volume estimation of brains of 5 trisomic TS1Yah mice

ACKNOWLEDGMENT

The authors thank the French Region Centre for their financial support that made this work possible.


A Marchadier

C Vidal

J.-P Tafani

S Ordureau

R Lédée

C Léger

Quantitative CT imaging for adipose tissue analysis in mouse model of obesity

Keywords: small animal imaging, medical CT, adipose tissue, gaussian mixture model, expectation maximization, level set segmentation

In obese humans CT imaging is a validated method for follow up studies of adipose tissue distribution and quantification of visceral and subcutaneous fat. Equivalent methods in murine models of obesity are still lacking. Current small animal micro-CT involves long-term X-ray exposure precluding longitudinal studies. We have overcome this limitation by using a human medical CT which allows very fast 3D imaging (2 sec) and minimal radiation exposure. This work presents novel methods fitted to in vivo investigations of mice model of obesity, allowing (i) automated detection of adipose tissue in abdominal regions of interest, (ii) quantification of visceral and subcutaneous fat.

For each mouse, 1000 slices (100µm thickness, 160 µm resolution) were acquired in 2 sec using a Toshiba medical CT (135 kV, 400mAs). A Gaussian mixture model of the Hounsfield curve of 2D slices was computed with the Expectation Maximization algorithm. Identification of each Gaussian part allowed the automatic classification of adipose tissue voxels. The abdominal region of interest (umbilical) was automatically detected as the slice showing the highest ratio of the Gaussian proportion between adipose and lean tissues. Segmentation of visceral and subcutaneous fat compartments was achieved with 2D ½ level set methods. Our results show that the application of human clinical CT to mice is a promising approach for the study of obesity, allowing valuable comparison between species using the same imaging materials and software analysis

INTRODUCTION

In human, the accurate determination of adipose tissue deposits throughout the body is an important issue in medical analysis of obesity. Different types of adipose tissue at different locations (subcutaneous versus visceral) play different physiological roles. Visceral adipose tissue is a significant risk factor for coronary artery diseases or diabetes [START_REF] Seidell | Regional distribution of muscle and fat mass in men -new insight into the risk of abdominal obesity using computed tomography[END_REF] and more predictive of obesity-induced pathologies than total or subcutaneous adipose tissues. MRI and CT imaging techniques allow accurate and non invasive analysis of adipose tissue distribution. These imaging modalities differ in term of cost, reproducibility, accuracy and acquisition time. Both MRI and CT, have been used in 2D [START_REF] Summe | Sex differences in visceral adipose tissue volume among African Americans[END_REF][START_REF] Shen | Visceral adipose tissue : relations between single-slice areas and total volume[END_REF][START_REF] Sussman | Automated fat measurement and segmentation with intensity inhomogeneity correction[END_REF][START_REF] Positano | Accurate segmentation of subcutaneous and intermuscular adipose tissue from MR images of the thigh[END_REF][START_REF] Ellis | Visceral and Subcutaneous Adiposity Measurements in Adults: Influence of Measurement Site[END_REF] and in 3D [START_REF] Al-Attar | Semiautomated segmentation and quantification of adipose tissue in calf and thigh by MRI: a preliminary study in patients with monogenic metabolic syndrome[END_REF][START_REF] Brennan | Rapid automated measurement of body fat distribution from whole-body MRI[END_REF][START_REF] Jin | Segmentation and Evaluation of Adipose Tissue from Whole Body MRI Scans[END_REF][START_REF] Jørgensen | Unsupervised Assessment of Subcutaneous and Visceral Fat by MRI[END_REF][START_REF] Makrogiannis | Automated Quantification of Muscle and Fat in the Thigh from Water-, Fat-and Non-suppressed MR Images[END_REF][START_REF] Furukawaa | Evaluation of whole-abdominal fat volume by 700-slice CT scanning and comparison with the umbilical fat area anthropometric indices[END_REF] in human. Better assessment of adipose tissue distribution is provided using a 3D region of interest than using a single slice [START_REF] Ellis | Visceral and Subcutaneous Adiposity Measurements in Adults: Influence of Measurement Site[END_REF].

Equivalent methods in murine models of obesity are still lacking. Micro-MRI technique was used for adipose tissue analysis in rodent but involves prohibitive costs and long acquisition time [START_REF] Tu | Automatic tracing and segmentation of rat mammary fat pads in MRI image sequences based on cartoon-texture model[END_REF][START_REF] Johnson | Reproducible MRI measurement of adipose tissue volumes in genetic and dietary rodent obesity models[END_REF]. Current small animal micro-CT involves long-term X-ray exposure precluding longitudinal studies [START_REF] Hildebrandt | Validation of a high-resoltuion X-Ray computed tomography system to measure murine adipose tissue depot mass in situ and longitudinally[END_REF][START_REF] Judex | Quantification of adiposity in small rodents using micro-CT[END_REF][START_REF] Luu | In vivo quantification of subcutaneous and visceral adiposity by micro-computed tomography in a small animal model[END_REF]. We have overcome this limitation by using a human medical CT which allows accurate and very fast 3D imaging and minimal radiation exposure. The resolution, reaches in the last generation of medical CT, is similar to the resolution needed to image adipose tissue in small rodent [START_REF] Judex | Quantification of adiposity in small rodents using micro-CT[END_REF].

For both human and small animal, manual adipose tissue quantification is often performed using two fixed level of thresholds [START_REF] Hildebrandt | Validation of a high-resoltuion X-Ray computed tomography system to measure murine adipose tissue depot mass in situ and longitudinally[END_REF][START_REF] Judex | Quantification of adiposity in small rodents using micro-CT[END_REF][START_REF] Luu | In vivo quantification of subcutaneous and visceral adiposity by micro-computed tomography in a small animal model[END_REF]. Several automatic methods have been developed to classify adipose tissue voxel using connectedness map [START_REF] Al-Attar | Semiautomated segmentation and quantification of adipose tissue in calf and thigh by MRI: a preliminary study in patients with monogenic metabolic syndrome[END_REF] or FCM [START_REF] Sussman | Automated fat measurement and segmentation with intensity inhomogeneity correction[END_REF][START_REF] Positano | Accurate segmentation of subcutaneous and intermuscular adipose tissue from MR images of the thigh[END_REF] and to segment visceral and subcutaneous compartments [START_REF] Sussman | Automated fat measurement and segmentation with intensity inhomogeneity correction[END_REF][START_REF] Positano | Accurate segmentation of subcutaneous and intermuscular adipose tissue from MR images of the thigh[END_REF][START_REF] Jørgensen | Unsupervised Assessment of Subcutaneous and Visceral Fat by MRI[END_REF][START_REF] Tu | Automatic tracing and segmentation of rat mammary fat pads in MRI image sequences based on cartoon-texture model[END_REF] We present an automatic method that classify adipose tissue id 3D, segment visceral and subcutaneous adipose tissue and detect the abdominal region of interest. Our method uses Gaussian mixture model, expectation maximization algorithm like Positano's method, and Chan-Vese based level sets.

MATERIAL AND METHODS

Method overview

A human medical CT was used for imaging adipose tissue in mice. The global data processing of the data included three studies : 1) a Gaussian mixture model (GMM) designed to identify voxels of adipose tissue and lean tissue, 2) analysis of GMM proportions along the mouse body in order to automatically define the region of interest and 3) segmentation of visceral versus subcutaneous adipose tissue using level sets method (Figure 1).

Mouse model of obesity

Male and female Swiss OF1, 10 weeks old, were obtained from the Charles River laboratory (France). The control group (n=8) received normal diet. Experimental obesity was produced in mice fed with highly energetic food pellets (A4010 rodent food pellets, Safe, Augy, France). Group1 mice (n=4) received a controlled amount of pellets. Group2 mice (n=4) was fed ad libitum (group2). For the CT imaging, mice were briefly anesthetized with isofluorane.

Scanning procedures

We used a medical multi-detector row CT (Aquilon 64, Toshiba Medical Systems Corp., Tochigi, Japan). The X-ray source was set with a voltage peak of 135 kV and a tube current of 400mAs. The imaging time of the mouse whole body was two seconds allowing the acquisition of 1000 slices. The resolution of the reconstructed image was a non-isotropic voxel of 160x160x100 µm. Figure 2 shows different views obtained from the same scan illustrating various organs and tissues visualization.

Voxel classification of fat and lean tissues

We focused the analysis on a subpart of the Hounsfield curve [-400 HU ; 300 HU] to avoid contribution of other tissues such as lung or bone. We used a mixture model of three Gaussians, two for adipose and lean tissues and one for others tissues. Each Gaussian can be expressed as:

(1

)
where µ is the mean and the variance of the Gaussian function . The mixture model is the linear combination of three Gaussian functions as described below: [START_REF] Summe | Sex differences in visceral adipose tissue volume among African Americans[END_REF] where G represents the Gaussian mixture model, is the proportion corresponding to the Gaussian function . The three Gaussian proportions are driven by the following equation:

(3) For each Gaussian, three parameters were estimated by means of the Expectation Maximization algorithm (EM). EM is an iterative method which alternates between performing an expectation (E) step and a maximization (M) step [START_REF] Dempster | Likelihood from incomplete data via the EM algorithm[END_REF]. The parameter values used at the initial step are shown in The Figure 3 illustrates the Hounsfield curve obtained from a slice in the abdominal area and the corresponding fitted curve obtained with the Gaussian mixture model (GMM). This model was applied on each 2D transverse slices along the animal as shown in Figure 4. Analysis of adipose tissue volumes was performed with two methods : 1) the GMM method and 2) the method with a Fixed Thresholding (FT method) range of [-250 HU ; -50 HU].

Automatic detection of the region of interest

We defined a reference abdominal slice of interest (corresponding to umbilical in human) by the image showing the highest amount of adipose tissue compared to amount of lean tissue. This slice was automatically detected using the proportions of the GMM previously computed as described in the following equations :

         i lean i fat i R , ,   (4) ) ( max i i R I  (5)
where R i represents the ratio of slice i, α adipose,i the proportion parameter of the adipose tissue Gaussian model of slice of index i, and α lean,i the proportion parameter of the lean tissue Gaussian model of slice of index i and I the index of the detected slice. The distribution of ratio R along the body is shown in figure 5.

Segmentation of visceral and subcutaneous adipose tissue

Segmentation of visceral adipose tissue (SAT) and subcutaneous adipose tissue (VAT) was performed using active contours. Our method was based on level set without edges [START_REF] Chan | An active contour model without edges[END_REF] and used the dense implementation of the Chan and Vese method [START_REF] Mosaliganti | Level Set Segmentation: Active Contours without edge[END_REF] available in the Insight toolkit (ITK) in 2D to segment the fat compartments limited by the abdominal muscular wall and the skin. The segmentation protocol was initialized on the slice of interest and then extended to the adjacent slices as a 2 1/2 method.

RESULTS

Analysis of whole body adipose tissue volumes

Figure 6 shows the repartition of adipose tissues in obese and control mice using the GMM method and the FT method.

Figure 6 : Fat tissue visualization of a mouse of each group with the GMM method and with the FT method

The quantitative analysis confirms the imaging data showing higher adipose tissue in obese mice than in control mice (Table2). 2: Mean volumes of adipose tissue (in mm3) computed for each group of mice, using FT and GMM The comparative analysis shows minor differences between GMM and FT measurements. When adipose tissue volumes are correlated with the mouse weight, we obtain slightly better correlation with the GMM method than with the FT method (Table 3). 3: Coefficient of correlation between the weight and the adipose volumes computed with FT and GMM As the whole, these results demonstrate the validity of the GMM method as compared to the FT method.

Group

Group

Automatic detection of the region of interest

In our protocol, the automatic localization of the slice of interest was based on the relative amounts of adipose and lean tissues. Interestingly, this criteria corresponded to skeletal landmark : the lumbosacral junction of L6 and S1 vertebras. The mean distance between this landmark and the detected slice was of 1.01 +/-0.45 mm. Figure 7 illustrates the ROI compared to others regions. 4: Coefficient of correlation between the weight and the adipose volume for the three regions of interest. The adipose tissue volumes used for the correlation were computed with the automatic thresholds.

The volumes of adipose tissue were computed on the three regions showing a higher amount of fat tissue in group1 and group2 as compared to control (respectively : 1107.1 mm 3 , 1396.3 mm 3 , 816.9 mm 3 ). Adipose tissue volumes in the ROI were better correlated with the body weight than the two others regions (Table 4).

This result reinforces the validity of the choice of the selected region.

Segmentation of visceral and subcutaneous compartments

Figure 8 illustrates the visceral and subcutaneous adipose tissue compartments. The segmentation protocol did not show any obvious correlation between weight and VAT, SAT and ratio VAT/SAT (Table 5).

Group

CONCLUSIONS

This work presents novel methods fitted to in vivo investigations of mice model of obesity, allowing (i) automated detection of adipose tissue in abdominal regions of interest, (ii) automated detection of the region of interest, (iii) quantification of visceral and subcutaneous fat.

Our results show that medical CT imaging combined with automatic image analysis provide precise and reproducible quantification of adipose tissue in mice in vivo allowing repetitive examinations for longitudinal studies. This new approach in experimental animals is of particularly interest for the study of obesity, allowing valuable comparison between species using the same imaging materials and software analysis

Fig 1 :

 1 Fig 1: Fat tissue classification process

Figure 2 :

 2 Figure 2 : Different views of the same CT scan of a mouse : skin (a), lung (b), adipose tissue (c), mineralization (e) and bone (e) visualisations

Fig 3 :

 3 Fig 3: On the left, example of 2D slice (a) and its localization on whole body mouse. On the right, Hounsfield curve of the slice and corresponding Gaussian mixture model of the distribution of fat and lean tissues (c).

Fig 4 :

 4 Fig 4: Serial slice associated mixture model of the whole body.

Figure 5 :

 5 Figure 5 : The bottom curve shows the ratio R for each slice. The slice of interest has been detected around the slice numbered 300. The up picture illustrates the anatomical localization of the slice. The Region Of Interest (ROI) included 200 slices (20 mm) centered on the detected slice. This ROI was compared to previously described regions of interest [] [], extracted manually in vertebral levels : 1) L1-L5 and 2) L5-C1.

Figure 7 :

 7 Figure 7 : Localization of the detected slice and comparison of the region with others area used in literature (left) ; detected slice of interest (right)

Fig 8 :

 8 Fig 8 : Segmentation of fat compartments limited by the abdominal wall (visceral fat, yellow) and the skin (subcutaneous fat, blue) in the reference slice.

Table 1 :

 1

		Mean	Variance	Proportion
	Gaussian 1 : Fat	-135 HU	2500 HU	0.3
	Gaussian 2 : Lean	80 HU	2500 HU	0.3
	Gaussian 3 : Others tissues	300 HU	2500 HU	0.3

Table 1 :

 1 Values of the initialization of the parameters of the Gaussian model

Mice wheight Total adipose tissue volume

			in mm3
		in g	FT	GMM
	control	33.9 +/-3.5	2212.4 +/-583.6	2063.2 +/-365.5
	group1	40.1 +/-3.8	2795.2 +/-491.5	3000.0 +/-377.5
	group2	42.8 +/-3.3	3564.0 +/-530.2	3308.8 +/-120.1
	Total	38.5 +/-5.1	2768.8 +/-692.8	2725.9 +/-629.6
	Table

Table 5 :

 5 Correlation coefficients between weight and fat volumes in subcutaneous and visceral adipose tissue area and for the ration VAT/SAT This result may be due to our animal model of obesity which exhibits relatively limited visceral fat accumulation as compared to genetically modified mice model of obesity.

			Coefficient of correlation
			weight/fat volume
		VAT	SAT	Ratio VAT/SAT
	control	0.410	0.580	0.061
	group1	0.539	0.682	-0.857
	group2	0.842	0.142	0.876
	Total	0.694	-0.349	0.008


Arnaud Marchadier

Catherine Vidal

Sylvain Ordureau

Roger Lédée

Christophe Léger

M Young

M E Goldberg

3D visualization and quantification of bone and teeth mineralization for the study of osteo/dentinogenesis in mice models

Keywords: Times Roman, image area, acronyms, references

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

INTRODUCTION

Research on bone and teeth mineralization in animal models is critical for understanding human pathologies. Genetically modified mice are highly valuable models for research on bone and tooth mineralization, and for better understanding of human bone and teeth pathologies [START_REF] Ameye | Mice deficient in small leucine-rich proteoglycans: novel in vivo models for osteoporosis, osteoarthritis, Ehlers-Danlos syndrome, muscular dystrophy, and corneal diseases[END_REF]. Current investigations on mice dental and skeletal phenotype use destructive and time consuming methods such as histology and scanning microscopy. Micro-CT imaging is quicker and provides high resolution qualitative phenotypic description.

However reliable quantification of mineralization processes in mouse bone and teeth are still lacking. We have established novel CT imaging-based software for (i) acurate qualitative and quantitative analysis of mouse mandibular bone, (ii) 3D reconstruction of molar pulp chamber, providing a unique tool for the quantitative evaluation of dentinogenesis.

MATERIALS AND METHODS

Method overview

Using a micro-CT normally devoted to industrial applications, we obtained high resolution 3D images of mandibles. Before analysis of these data, we performed a histogram calibration in order to have the best correlation between grey levels and tissues mineralization. Qualitative and quantitative analysis were then computed allowing complex visualization and 3D measurements as described in Figure 1.

Mandible samples

All experiments were performed under an institutionally approved protocol for the use of animals in research (#NIDCR-IRP-09-536). Data obtained from mandibles of 3 week-and 10 week-old Fmod-KO mice were compared with age matched wild type mice (WT) (16 mice/group).

Data were obtained from mandibles of mice lacking the Fibromodulin gene (KOFM mice), which is knowed to be involved in mineralization processes [START_REF] Goldberg | Fibromodulin-deficient mice display impaired collagen fibrillogenesis in predentin as well as altered dentin mineralization and enamel formation[END_REF], and compared to wild type mice (16 mice/group).

Micro-CT imaging

For high-resolution 3D imaging of mandibles, we use a micro-CT originally devoted to industrial applications (Viscom, X8060 NDT). The X-Ray source, designed with a molybdenum target, was set up with a voltage of 60kV and a current of 200µA. The exposure time was 1sec per frame. Mandible 3D images had an isotropic resolution of 13µm and a dynamic range of 16 bits.

Calibration of grey levels

Contrasts in grey levels of CT images rely on differential attenuations of X-Rays through the material []. Theoretically, only monochromatic X-Ray beam, as reach in a synchrotron, can lead to a direct linear relation between grey level and density []. As a conventional X-Ray source, the micro-CT X-Ray beam is polychromatic and may give rise to artifactual variations of the grey scale. The principal artifact : beam hardening effect, leads to variations of grey level with the size of sample and need to be addressed carefully []. To eliminate this artifact, different method are available [] and in particular for scanning bone sample. Moreover, assuming a beam hardening correction, it has been previously reported that the relation between grey level and densities could be calibrated by a linear approximation [].

During the tomographic reconstruction process, we used an option of the dedicated software (Vol64, Viscom) that allowed beam hardening correction. In order to obtain an accurate comparison between the mandibles of the two groups of mice, we performed a systematic calibration of the CT images. Two references were used (i) the enamel at the tip of the incisor, the density of which is constant in all mice, (ii) the background grey level. Histograms of 3D images were then corrected if necessary according to the defined references before proceeding to the 3D visualization. This correction protocol consisted on a shift and scale operation. Implementation of volume rendering with a ray casting method, available in the VoxBox software (UsefulProgress), was used for advanced visualization of mandibles. Comparison between control and defective mouse mandibles was made by applying the same transfer function for each 3D data, thus allowing to detect shape, colour and density discrepancies.

We defined four transfer functions : one for opacity and three for color in a HSV space as illustrated in Figure 3 :

Pulp chamber volume

The pulp chamber of the molar refers to the inner most portion of the tooth that contains the pulp. It constitutes a 3D closed region suitable for region growing method segmentation. We used the connected threshold image filter available in the Insight Toolkit [START_REF] Ameye | Mice deficient in small leucine-rich proteoglycans: novel in vivo models for osteoporosis, osteoarthritis, Ehlers-Danlos syndrome, muscular dystrophy, and corneal diseases[END_REF] to extract the pulp chamber volume.

Grey level : good contrast air/dentin Three parameters were set as input: two thresholds of grey level and a 3D seed point.

RESULTS

Qualitative results

Based on the volume rendering visualization, our first result was to identify three distinct parts forming the mandible as described in Figure 6. The distal or rostral part differed from the central part and represents the most mineralized module. The central part was subdivided into an upper zone containing the three molars implanted in an alveolar bone, whereas in the lower part the basal bone encircled the incisor. The proximal part included endochondral-derived domains: the angular, condylar and coronoid processes. The dark brown color indicates that this area is thinner and consequently less mineralized than the rest of the mandible. With this qualitative analysis of the mineralization distributions, we globally conclude that a severe hypomineralization was induced in KOFM mice. An example of this qualitative difference is shown in Figure 7 : In details, the main target for genetic deficiency was the proximal part of the mandible including condylar, coronoïd and angular processes, which was less dense KOFM compared with the WT. This was also detectable in the alveolar bone (upper middle) but became undetectable in the basal and rostral parts of the mandible.

Volume rendering of the whole mandible in a sagittal orientation bring principally two results : i) that the mandible is a mosaic formed by different domains with an unequally distributed thickness and mineralization pattern and ii) the proximal and alveolar parts of the mandible constitute the specific targets in KOFM mice.

Moreover, transverse virtual slices across the mandible allowed to analyze the inner structure of the alveolar bone surrounding the molar roots. The virtual slices, through the junction of the first and second molars, were similar and even more accurate than those obtained with scanning electron microscopy, e.g. revealing the detailed range of mineralization defects of the alveolar bone owing to the color code as illustrated in Figure 8. The marrow spaces of the medulla (trabecular zone) were larger in the KOFM mice compared with the WT. No difference was noted for the cortices. Therefore, bone hypomineralization was associated with the enlargement of endosteous spaces. There was no spontaneous repair or rescue of this phenotype during the aging process, at least for the period of time studied.

Quantitative results on molars

The segmentation method allowed the 3D reconstruction and quantification of the entire pulp chamber volume of molars. As described in Figure 9, KOFM mice showed a statistically significant larger volume in as compared to control mice thus revealing for the first time a defect in dentinogenesis. Moreover, several morphological changes, as for example the length of the roots, were reported but were not quantified (Figure 9 on right).

CONCLUSIONS

This new method is highly powerful for the study of oro-facial mineralization defects in mice models, complementary and even competitive to current histological and scanning microscopy approaches.

This methodolgy constitutes a new tool for the study of mineralization defects in genetically modified mice, providing rapid and high resolution imaging and quantitative correlates of the phenotypic alterations and comparison with wild type control mice.

Figure 1 :

 1 Figure 1: Method overview for qualitative and quantitative analysis of mandibles

Figure 2 :

 2 Figure 2 : Histograms of the whole mandible dataset before and after correction 2.5 Visualization with volume rendering Volume rendering is a well known technique used to visualize 3D discretely sampled datasets []. This technique is particularly suitable to render images acquired by a CT, MRI, or MicroCT []. In contrast to surface rendering, that need first primitive extraction (as for example an isosurface), volume rendering allows direct voxels visualization with opacity and color.

Figure 3 :

 3 Figure 3 : Transfer functions for opacity and HSV used for volume renderingTwo different views were computed : a sagittal orientation of the whole mandible and a thick region of interest in a frontal orientation located between the first molar and the second molar (Figure4).

Figure 4 :

 4 Figure 4 : Sagittal view of a mandible (left) and junction between the first and the second molar (right)

Fig 5 :

 5 Fig 5: Region growing segmentation of the pulp chamber of the second molar

Figure 6 :

 6 Figure 6 : Volume rendering visualization showing the different regions of a mandible

Fig 7 :

 7 Fig 7: Volume rendering visualization of the mandibles of a KOFM (left) and a control mouse (right). The KOFM mandible showed reduced mineralization in the condyl (posterior) part of the mandible, alveolar bone area (around the roots of the molars) and in the basal bone around the incisor

Fig 8 :

 8 Fig 8: Comparison of 2D slices across the mandible using scanning microscopy (left) and micro-CT (right).Micro-CT images allowed precise visualization of hypomineralized areas in KOFM slices, rivalizing with scanning microscopy

Fig 9 :

 9 Fig 9: Statistical results of pulp chamber volumes comparison (left) and 3D reconstruction (right) of pulp chamber volumes of KOFM (grey) and control mice (yellow)

ACKNOWLEDGEMENTS

This research was supported by the IRP Program of the NIDCR/NIH and by the Fondation de l'Avenir, Etude N° ETO-576.