Dahia Chibouti

Benoît Trouette

Eric Chénier
email: eric.chenier@univ-eiffel.fr

{dahia Chibouti;Benoit

Trouette;

Jump in the conduction heat flux at the gas/solid interface in micro-channels

Keywords: Heat transfer, Micro-channel, Microfluidics, Rarefied gas, Slip flow, Temperature jump, Heat flux jump

In micro-channel gas flows, velocity slip and temperature jump must be applied at the fluid/solid interface to take the thermodynamic non-equilibrium into account in the vicinity of the wall. In this work, simulations at molecular scale have been carried out and shown a jump in the conduction heat flux, in agreement with the boundary condition proposed by Maslen [10] in continuum mechanics.

Introduction

The understanding of fluid flows and heat transfer at scales of few microns is of paramount importance since the development of Micro-electro-mechanical Systems (MEMS). When the characteristic length scale of the flow decreases to reach about 10 times the free mean path λ of the gas molecules, the usual Navier-Stokes and energy equations remain valid in the core flow, but careful attention must be paid on the gas/solid boundary conditions. In a fluid layer of order λ close to the wall, the so-called Knudsen layer, the gas molecules are no more in a local thermodynamic equilibrium: the macroscopic model based on the continuity of the velocity, temperature and heat flux at the interface between the gas and the wall falls and slip or jump conditions have to be used, like those independently proposed by Navier [12] and later on by Maxwell [START_REF] Maxwell | On stresses in rarified gases arising from inequalities of temperature[END_REF]. Whereas the boundary conditions applied on the velocity and temperature are quite well es-tablished in the 'gaseous microfluidic' community, that associated to the jump in the conduction heat flux at the wall due to the viscous friction, and first proposed by Maslen [START_REF] Maslen | On Heat Transfer in Slip Flow[END_REF] in 1958, has often been neglected, forgotten or misunderstood, including in late papers. In the recent works by X. Nicolas et al. [START_REF] Nicolas | Thermal boundary conditions for convective heat transfer of dilute gases in slip flow regime[END_REF], a review of few papers using or ignoring the Maslen's flux condition is performed, followed by its proof which relies on the energy flux conservation principle at the macroscale. The same authors also showed that the use of the Maslen's flux condition explains the very small Nusselt number values in experiments performed on microfluidic gaseous flows [START_REF] Nicolas | Revisited analysis of gas convection and heat transfer in micro channels: Influence of viscous stress power at wall on Nusselt number[END_REF]. However, their mathematical approach is questionable since it relies on the assumption of the continuum mechanics which is valid far from the boundaries but seems no more applicable in the Knudsen layer. This short contribution aims to check if the Maslen's flux condition at the gas/solid interface is correct by getting rid of the continuum model and by simulat-ing the fluid flow and energy transfer at the molecular scale.

Continuum description

An established argon (Ar) fluid flow between two infinite (x, y)-parallel solid walls of platinum (Pt) is considered (Fig. 1). Walls are distant from H and their thickness is e. The flow is created by a volume force g = g x e x acting as a uniform pressure gradient in xdirection. The platinum wall starting at z = H/2 is thermally isolated, i.e. the conduction heat flux q z = 0, whereas the one ending at z = -H/2 is isotherm, at the temperature T 1 . For the compressible fluid flow in the micro-channel of length L, driven by the pressure difference p in -p out , the density decreases all along the channel. From the mass flow rate conservation, the velocity magnitude increases, and consequently, the viscous heating also. If the length scale associated with the pressure gradient is assumed much larger that the channel height, namely (p in + p out )/(p in -p out ) 2H/L, then the density variations over a length scale H of the channel can be disregarded. Therefore the axial dependence of all the macroscopic variables are locally negligible, at least far from the inlet and outlet regions [START_REF] Nicolas | Revisited analysis of gas convection and heat transfer in micro channels: Influence of viscous stress power at wall on Nusselt number[END_REF]. Thus, in what follows, the thermodynamics quantities are assumed to depend only on the transverse coordinate z. From a macroscopic point-of-view, and as discussed in [START_REF] Nicolas | Thermal boundary conditions for convective heat transfer of dilute gases in slip flow regime[END_REF], the total energy flux q (c) tot in the fluid region can be split into three contributions:

q (c) tot = q k + q µ + q h = -k f ∇T -( τ • u) + ρhu ( 1
)
where q k = -k f ∇T is the conduction heat flux with k f the fluid conductivity, q µ = -( τ • u) is the energy flux produced by the viscous friction and q h = ρhu is the convection enthalpy flux. Since the velocity simply writes u = u x (z)e x , the total energy flux in z-direction, q

(c)
tot , reduces into two terms in z-direction:

q (c) tot = -k f ∂T ∂z -u x µ ∂u x ∂z (2)
The Maslen's flux boundary condition [START_REF] Maslen | On Heat Transfer in Slip Flow[END_REF] is to equalize the total energy fluxes on both sides of the fluid/solid interface. In the wall region, the total energy flux reduces in its conduction contribution:

-k s ∂T ∂z H/2 = -k f ∂T ∂z H/2 -u x µ ∂u x ∂z H/2 (3)
where k s is the solid conductivity and µ the fluid dynamic viscosity. In order to check the validity of Eq. ( 3) in the Knudsen layer, simulations are carried out at the molecular scale.

Microscopic approach

Solid walls are composed of platinum atoms disposed on a FCC111 lattice. The dynamics of all atoms, either argon confined between the solid walls or platinum, relies on the classical Newton's law. The trajectory of each atom i, or particle, is governed by F i = m i ri where ri is the acceleration vector, m i is the mass and F i is the resultant of external forces acting on particle i. With r i is the particle position, g is the volume force which drives the fluid flow, and assuming only binary interactions given by the pairwise potential U α/β , the resultant F i is developed as

F i = -∇ g • (r i -r i ) + j U α/β (r ij
) for an argon particle and is reduced to

F i = -∇ j U α/β (r ij
) for a platinum one. The classical Lennard-Jones model [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF] is chosen to describe the argon/argon, platinum/platinum and argon/platinum interactions:

U α/β (r ij ) = 4ε α/β σ α/β r ij 12 - σ α/β r ij 6 (4)
where α and β denote the Ar-or Pt-atoms, r ij = r ij is the intermolecular distance between the i and j particles, r ij = r j -r i , ε α/β is the potential welldepth and σ α/β is the distance where the potential cancels out. The choice of argon as a fluid is justified by the binary interactions and spherical approximation induced by the pairwise Lennard-Jones potential choice. This latter allows a good restitution of the dynamical and thermal properties of mono-atomic gases, and to some extend, of non-spherical molecules for dense cases [START_REF] Bugel | Thermal conductivity of the Lennard-Jones fluid: An empirical correlation[END_REF][START_REF] Galliéro | Molecular Dynamics Study of the Lennard-Jones Fluid Viscosity: Application to Real Fluids[END_REF]. Like in continuum mechanics (Sec. 2), the total energy flux q (m) tot can be defined from a microscopic description. Using the Irving-Kirkwood formulation [START_REF] Fan | Force and heat current formulas for many-body potentials in molecular dynamics simulations with applications to thermal conductivity calculations[END_REF][START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF], it reads:

q (m) tot = 1 2V N i m i ṙ2 i ṙi + 1 2V N i N j =i U α/β (r ij )ṙ i - 1 2V N i N j =i r ij [f ij • ṙi ] ( 5
)
where ṙi is the velocity of particle i,

f ij = -∇U α/β (r ij )
is the interaction force between particles i and j and N is the total number of particles in the volume V where q

(m)
tot is evaluated. The first and second terms of the r.h.s. are related to summations of kinetic and potential energy carried by a particle i, while the third term represents the energy transfer by the pressure work [START_REF] Maruyama | Molecular dynamics methods in microscale heat transfer[END_REF]. It can be noticed that no contribution involving the volume force g exists, simply because the velocity and position are independent variables.

In order to impose a zero conduction heat flux in the solid platinum wall ranging from z = H/2 to z = H/2 + e, the fluid region is replicated at the wall ending, at z = H/2 + e. That trick also allows keeping the periodicity of molecular simulations in all space directions. Since the potential interaction between two particles decreases quickly with their distance, to become negligible in the total amount of energy, the computation cost is highly reduced by introducing a truncation length r c . In practice, the interaction force derived from Eq. ( 4) is set to zero, if r ij > r c . The truncation radius r c is set to r c = 2.5 σ.

The simulations are carried out on boxes of dimensions

x × y × z . With the fluid zone replica, the length is finally z = 2(H + e). The lengths i , i = x, y, are multiples of the gaps of the FCC111 solid lattice. They also fulfil the condition i ≥ 3 × r c in order to manage efficiently the neighbourhood of each particle (Cell and Verlet lists [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF]). The height H is such that H ≥ C × min( x , z ) and is adjusted to fit the desired fluid density ρ in the volume defined by V = x × y × H. The constant C, initially fixed at C = 10, can be adapted to change the Kn number. In order the particles in both fluid regions do not interact between each other, the wall thickness e is such that e ≥ r c and e is a multiple of the vertical gap of the solid FCC111 lattice. The temperature of the left wall is maintained at T 1 using a Langevin thermostat [START_REF] Hünenberger | Thermostat Algorithms for Molecular Dynamics Simulations[END_REF], while that of the opposite wall is left free. The Velocity-Verlet algorithm [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF] is used to perform the time integration with a time step fixed at ∆t/t 0 = 5 × 10 -3 .

To calculate the macroscopic quantities ρ, u, T and q (m) tot , time and space averages are performed on n z uniform sub-domains covering the channel height H and over time intervals τ a . The converged estimate of the variance of the ensemble average is computed by a block average technique [START_REF] Flyvbjerg | Error estimates on averages of correlated data[END_REF][START_REF] Frenkel | Understanding Molecular Simulation[END_REF].

Results and discussions

The molecular dynamics simulations are based on a Lennard-Jones potential interactions. Therefore, the pressure, the temperature and the density are solutions of a specific equation of state [START_REF] Kolafa | The Lennard-Jones fluid: An accurate analytic and theoreticallybased equation of state[END_REF]. Likewise, the dynamic viscosity µ and the thermal conductivity k f are estimated through the correlations by Galliéro et al. [START_REF] Galliéro | Molecular Dynamics Study of the Lennard-Jones Fluid Viscosity: Application to Real Fluids[END_REF], and Bugel and Galliéro [START_REF] Bugel | Thermal conductivity of the Lennard-Jones fluid: An empirical correlation[END_REF]. The results are then presented as a function of the physical quantities, namely the mean density ρ and the volume force g x . The magnitude order of the slip velocity and jump temperature is controlled by the Knudsen number

Kn = λ 2H ( 6
)
with λ is the mean free path for molecules interacting through long-range potentials, which is expressed with the formulation by Cercignani [START_REF] Cercignani | Slow rarefied flows: Theory and application to Micro[END_REF]:

λ = µ ρ πm 2k B T ( 7
)
where m is the molecular mass, T is the temperature of the gas, and k B is the Boltzmann constant.

For the sake of simplicity, the physical parameters and the dynamic and thermal variables, averaged over the channel height, are gathered in Table 2 -8. The Knudsen number is Kn = 0.036, except the simulation with g x /g 0 = 0.002 where Kn = 0.0075 the Knudsen number Kn, the mean fluid density ρ and the volume force g x . The notation • denotes a spatial average over the channel height. All molecular dynamics simulations were carried out for a thermostated wall fixed at T 1 /T 0 = 1.5. The mean free path (Eq. 7), and therefore the Knudsen number (Eq. 6), is evaluated for a fluid flow without volume force g x = 0 (u x = 0), at the temperature equilibrium (T = T 1 ) and with a density measured in the fluid core, far from the boundaries.

Total energy flux analysis

The first simulation is carried out for a fluid density ρ/ρ 0 = 0.6 and a volume force g x /g 0 = 0.01. Based on the expression of the mean free path (Eq. 7), the Knudsen number ( 6) is Kn = 0.036. According to the fluid flow regimes, a velocity slip and a temperature jump are expected at the fluid/solid interface.

In order to compare the total energy flux at the molecular scale q (m) tot (Eq. 5) with its continuous counterpart q (c) tot (Eq. 2), the conduction heat flux q k and the viscous power density q µ have to be evaluated. To this end, the macroscopic velocity and temperature are first calculated by averaging the particle velocities and the local kinetic energy fluctuations (see Fig. 5 and Fig. 6 for the velocity and temperature profiles). Then, the macroscopic fluxes are approximated using the second order centred finite difference scheme to get the z-derivatives, combined with the dynamic viscosity and thermal conductivity provided by correlations [START_REF] Bugel | Thermal conductivity of the Lennard-Jones fluid: An empirical correlation[END_REF][START_REF] Galliéro | Molecular Dynamics Study of the Lennard-Jones Fluid Viscosity: Application to Real Fluids[END_REF]. In Fig. 2 are presented the different fluxes q

(m) tot , -2 -1.5 -1 -0.5 0 0.5 1 -1/2 -1/4 0 1/4 1/2 q / | q - k | z / H q tot (c) q tot (m) q k q µ
Figure 2: Normalized fluxes with respect to |q k |/q 0 = 2.74 × 10 -2 , the z-average of the conduction heat flux, for Kn = 0.036, ρ/ρ 0 = 0.6 and g x /g 0 = 0.01 q (c) tot , q k and q µ the viscous contribution to the energy equation. Each of them is normalized with respect to the absolute value of the mean conduction heat flux over the channel height. The error bars take into account the uncertainties on both, the calculated macroscopic data, the conductivity and viscosity correlations and the approximations of the velocity and temperature derivatives as well. The analysis of the curves confirms that the total energy flux computed with the molecular dynamics, q (m) tot , gets to zero at the fluid/solid interface z/H = 1/2, which is consistent with the conduction heat flux transmitted to the adiabatic solid wall. The reconstructed macroscopic flux q (c) tot exhibits a very good agreement with the molecular one q (m) tot in the core channel, and to a lesser extent, near the fluid/solid interfaces. In this problem, the viscous contribution q µ is far from being negligible, since it reaches almost half of the conduction heat flux. Its spatial distribution is anti-symmetric with respect to z = 0, which results from the symmetry of the nearly parabolic velocity profile (Fig. 5). Because the temperature profile is almost linear (Fig. 6), the conduction heat flux q k is nearly constant; the variations are mainly due to the conductivity dependence on the density and temperature. Finally, the macroscopic energy flux decomposition (2) is, as expected, accurate in the fluid core domain. On the other hand, differences are observed in the vicinity of the walls. At the adiabatic fluid/solid interface, the reconstructed flux is not as small as the molecular one, what may suggest that the macroscopic expression (2) fails to represent the right total energy flux, at least locally.

In order to examine thoroughly the departures between the total energy fluxes q (m) tot and q (c) tot , calculated by either the Irving-Kirkwood formulation (Eq. 5) or deduced from relations of the continuum mechanics (Eq. 2), their relative difference to q(m) tot is drawn for different densities, external forces and Knudsen numbers (Fig. 3). The reconstructed fluxes q z/H Kn = 0.0075, ρ = 0.6 , g x /g 0 = 0.002 Kn = 0.036 , ρ = 0.6 , g x /g 0 = 0.02 Kn = 0.036 , ρ = 0.6 , g x /g 0 = 0.01 Kn = 0.036 , ρ = 0.1 , g x /g 0 = 0.01 Kn = 0.036 , ρ = 0.05, g x /g 0 = 0.005 (m) tot (Eq. 5) and q (c) tot (Eq. 2) reconstructed with the molecular data. Filled (empty) symbols stand for points inside (outside) the Knudsen layer (|z|/H -1/2) < 2λ/H = 4Kn. For clarity reason, the error bars have been omitted. Reference values are given in Table 2 by less than 10% from the molecular ones q (m) tot in the bulk fluid domain. As discussed before, the discrepancies are mainly due to the reconstructed flux q (c) tot which relies on, both correlation laws for the thermal conductivity and dynamic viscosity evaluations, and approximations of the derivatives by the second order centred difference scheme. However, in a thickness region of about two mean free paths, the relative gap increases substantially, all the more that one gets closer to the wall. Looking at the density distributions (Fig. 4), growing oscillations develop towards the solid wall, approximately over the same length scale. These oscillations result from the clustering of the argon atoms close to the wall. Since the physical properties are density dependent, the dynamic viscosity µ and the thermal conductivity k f fluctuate in this zone as well.

This atomic patterning breaks the local thermodynamics equilibrium, which is required to get a macroscopic model for the fluid flow. The aim of the slip/jump relations is precisely to connect the macroscopic temperature and velocity variables on both side of the Knudsen layer, namely between the wall and the outer region of the Knudsen layer. The use of these boundary conditions, and more specifically the Maslen's flux condition, is the motivation of what follows.

Comparison between the molecular and continuum models

This section is devoted to the comparisons of the molecular dynamics results with those derived from the continuum model with constant physical properties, an incompressible fluid flow assumption and suitable boundary conditions. The simulations are carried out for volume force ranging from g x /g 0 = 0.0025 to 0.02 and a mean density ρ/ρ 0 = 0.6. The Knudsen number, evaluated at T 1 , is equal to Kn = 0.036.

The assumption of a constant density is acceptable in the core flow, where a weak decrease is observed due to the heating by the viscous friction (Fig. 4). On the 2). The number of bins for the molecular average calculation is n z = 100 other hand, at the vicinity of the boundaries, the density oscillates and departs up to about 75% of its bulk value. As already mentioned, this usual behaviour is related to the clustering of the argon atoms near solid wall due to the interaction parameters that make the Lennard-Jones potential attractive [START_REF] Nedea | Gas-surface interactions using accommodation coefficients for a dilute and a dense gas in a micro-or nanochannel: Heat flux predictions using combined molecular dynamics and Monte Carlo techniques[END_REF]. The oscillations are damped over a length of approximately 2λ that roughly corresponds to the Knudsen layer. The normalized velocity profiles, presented in Fig. 5, exhibit a Poiseuille like flow with slip at walls. Whereas the magnitude of the velocity and the dynamical slip depend on the applied volume force g x , they become quite identical once they are normalized by ūx , the averaged value of u x over the channel height (Table 2).

The superposition of the molecular data are fully consistent with the analytical function derived from the

0 0.2 0.4 0.6 0.8 1 1.2 1.4 -1/2 -1/4 0 1/4 1/2 u x / u - x z / H
g x /g 0 = 0.02 g x /g 0 = 0.01 g x /g 0 = 0.005 g x /g 0 = 0.0025 Eq. ( 8)

0.6 0.8 1 H/2-2λ H/2-λ H/2 z 0.6 0.8 1 H/2-2λ H/2-λ H/2
u x (z) ūx = 3(1 + 8σ µ Kn) -12(z/H) 2 2 + 24σ µ Kn (8)
with σ µ is the viscous slip coefficient classically used in the expression of the dynamical slip boundary condition:

u f -u w = σ µ λ (∇u • n wf )| f (9)
The indexes f and w indicate fluid and wall quantities evaluated at the fluid/solid interface and n wf its unit normal vector, directed from the wall to the fluid. In the present problem, u w = 0. Since the velocity profile is not constant in the channel section, the viscous friction leads to mechanical energy dissipation and a temperature increase by getting away from the lower wall at T 1 . This heating and the temperature jump are even more pronounced as the external force g x increases. However, like for the velocity, the dimensionless reduced temperature θ/θ b , with θ = T -T 1 and θ b = ρu x θ dz/ ρu x dz the dimensionless bulk temperature (Table 2), turns out to be independent on the external force g x (Fig. 6). Such a result is recovered by integrating the energy equation with the velocity field (8), the temperature jump at the lower boundary z/H = -1/2

T f -T w = ξλ (∇T • n wf )| f ( 10
)
with ξ the thermal jump coefficient, and the Maslen's flux condition (3) at the adiabatic wall z/H = 1/2. The analytical reduced temperature θ(z)/θ b then reads:

θ(z) θ b = A + B z H + C z H 4 (11)
with

A = a 0 + a 1 Kn
a 0 = 5, a 1 = 48σ µ + 32ξ, a 2 = 192σ µ ξ, d 0 = 43 20160 , d 1 = ξ 72 + 11σ µ 240 , d 2 = σ µ 4 (ξ + σ µ ), d 3 = σ 2 µ ξ
This analytical solution, which is based on the Maslen's flux condition, fits perfectly the molecular data. This observation suggests that the flux boundary condition [START_REF] Fan | Force and heat current formulas for many-body potentials in molecular dynamics simulations with applications to thermal conductivity calculations[END_REF], which can be interpreted as a jump in the conduction heat flux, permits to connect correctly the thermal quantities on both part of the Knudsen layer.

As well as the velocity and temperature profiles, the total energy flux becomes independent of the volume force g x , once the suitable normalization q

(m) tot /|q (m)
tot | is considered (Fig. 7). An analytical approximation | z / H g x /g 0 = 0.02 g x /g 0 = 0.01 g x /g 0 = 0.005 g x /g 0 = 0.0025 Eq. ( 12)

-0.3 -0.2 -0.1 0 0.1 H/2-2λ H/2-λ H/2 z -0.3 -0.2 -0.1 0 0.1 H/2-2λ H/2-λ H/2
q (a)
tot (z) of the molecular data is developed by substituting the velocity [START_REF] Kolafa | The Lennard-Jones fluid: An accurate analytic and theoreticallybased equation of state[END_REF] and temperature [START_REF] Maxwell | On stresses in rarified gases arising from inequalities of temperature[END_REF] functions into the macroscopic expression of the total energy flux (2):

q (a) tot (z) |q (a) tot | = (24σ µ Kn + 3)(z/H) -4(z/H) 3 12σ µ Kn + 1 -1 (12)
Again, the analytical expression provides a very good approximation of the results coming from the molecular dynamics simulations. It is worth to noticed that, as for the density profiles (Fig. 4), flux oscillations are observed at the molecular scale, over a length of about 2λ.

To point out the robustness of the Maslen's boundary conditions, the total heat fluxes q (m) tot , achieved by the molecular simulations and presented in Fig. 3, are reconsidered and compared to the analytical solutions q (a) tot (z) (Fig. 8). Taking into account the z / H Kn = 0.0075, ρ/ρ 0 = 0.6 , g x /g 0 = 0.002 Kn = 0.036 , ρ/ρ 0 = 0.6 , g x /g 0 = 0.02 Kn = 0.036 , ρ/ρ 0 = 0.6 , g x /g 0 = 0.01 Kn = 0.036 , ρ/ρ 0 = 0.1 , g x /g 0 = 0.01 Kn = 0.036 , ρ/ρ 0 = 0.05, g x /g 0 = 0.005 Figure 8: Relative gap between the total energy fluxes stemmed from molecular dynamics simulations and the analytical expression (12). Filled symbols stand for points inside the Knudsen layer (|z|/H -1/2) < 2λ, while empty ones denote measurement points in the bulk. For clarity reason, the error bars have been omitted. Reference values are given in Table 2 approximation level used to develop the continuum model, the analytical function q (a) tot (z) gives satisfying results, with relative gaps less than 10%. The accuracy is much better close to the adiabatic wall z/H = 1/2 (< 2%), where the total energy flux is controlled by the Maslen's condition. The slight deviation observed by moving away from the adiabatic wall is probably the result of error accumulations induced by the assumptions used to develop the analytical solution. Again, the comparisons between the analytical and molecular total energy fluxes emphasize the robustness of the Maslen's flux boundary condition (3) to achieve a correct approximation.

From the dimensionless writing of the Navier-Stokes and energy equations with their boundary conditions, 6 dimensionless parameters appear: the heat conductivity ratio k s /k f , the ratio between the volume and viscous forces measured by the Jeffreys number Je = g x H 2 /(µU ), the energy flux ratio between the viscous friction and the heat conduction, namely the Brinkman number Br = µU 2 /(k f ∆T ), and obviously the above introduced Knudsen number Kn and the slip and jump coefficients σ µ and ξ. Setting Je = 1 leads to define the scale for the velocity: U = g x H 2 /µ. Likewise, the equilibrium between the viscous dissipation and the heat conduction allows to define the scale for the temperature variation: ∆T = µU 2 /k f = g 2

x H 4 /µk f . Table 3 8), [START_REF] Maxwell | On stresses in rarified gases arising from inequalities of temperature[END_REF] and (12).

the Jeffreys, Brinkman and Knudsen parameters, and the slip and jump coefficients as well. For each fixed Knudsen, the numerical values σ µ and ξ result from a minimization process between the equations ( 8), ( 11) and ( 12) and the averaged molecular data, using the non-linear least-squares Marquardt-Levenberg algorithm. Since the fluid properties results from the Lennard-Jones interaction potential, the ranges of variation of the the Jeffreys and Brinkman parameters are moderate: the ratio between the maximal and minimal values are close to 2.

The numerical parameters involved in the molecular dynamics simulations, such as the cut-off radius r c , the time step ∆t, the lengths x and y , the interaction potential of the platinum atoms of the solid wall or the (σ, ε)-values of the Ar/Pt interaction coefficients do not alter the analyses presented below. Obviously, the changes in the Lennard-Jones interaction parameters may modify slightly the fluid flow and heat transfer, especially the dynamical slip and temperature jump magnitudes, but the velocity, temperature, and flux profiles remain similar; the total energy fluxes computed by the molecular simulation and the analytical model, using the Maslen's flux boundary condition [START_REF] Fan | Force and heat current formulas for many-body potentials in molecular dynamics simulations with applications to thermal conductivity calculations[END_REF], always cancel at the adiabatic solid wall z = H/2.

Conclusion

Results of molecular dynamics simulation, obtained for a Lennard-Jones fluid, have shown that the Maslen's boundary condition [START_REF] Maslen | On Heat Transfer in Slip Flow[END_REF], demonstrated in the framework of the continuum mechanics (see [START_REF] Nicolas | Thermal boundary conditions for convective heat transfer of dilute gases in slip flow regime[END_REF]), although questionable in the Knudsen layer, is proved to be efficient to connect the conduction heat fluxes between the wall and the fluid domain. Therefore, the total energy flux conservation at the fluid/solid interface introduces effectively a jump in the conduction heat flux which is equal to the energy flux produced by the viscous friction.

Figure 1 :

 1 Figure 1: Physical configuration

Figure 3 :

 3 Figure 3: Relative difference between the total energy fluxes q

Figure 4 :

 4 Figure 4: Density ρ(z)/ρ for different volume forces g x /g 0 and Kn = 0.036 (see Table2). The number of bins for the molecular average calculation is n z = 100

Figure 5 :

 5 Figure 5: Velocity u x (z)/ū x for different volume forces g x /g 0 and Kn = 0.036. The average values are provided in Table 2. The parameters for the analytical function is σ µ = 2.46 and ξ = 36.0. The number of bins for the molecular average calculation is n z = 25, except the encapsulated figures which are based on n z = 100

Figure 6 :

 6 Figure 6: Same caption as Fig. 5 -Reduced temperature θ(z)/θ b

Figure 7 :

 7 Figure 7: Same caption as Fig. 5 -Total energy flux q (m) tot (z)/|q (m) tot |.

Table 1 :

 1 Physical scales, expressions and values used in the molecular dynamics simulations, with k B is the Boltzmann constant

	Taking argon as reference for the physical
	properties, the interaction coefficients write ε Ar/Ar =
	1 ε, σ Ar/Ar = 1 σ, m Ar/Ar = m, ε P t/P t = 31.3 ε,
	σ P t/P t = 0.73 σ, ε Ar/P t = 0.66 ε, σ Ar/P t = 0.865 σ,
	m P t = 4.8833 m with ε = 1.656 × 10 -21 J, σ = 3.405 Å
	and m = 6.633 × 10 -26 kg. The reference scale expres-
	sions and their derived quantities are given in Table 1.
	Scale	Expression	Value
	mass (Ar)		m	6.633 × 10 -26	kg
	length		σ	3.405 × 10 -10	m
	energy		ε	1.656 × 10 -21	J
	time, t 0		mσ 2 /ε	2.155 × 10 -12	s
	temperature, T 0		ε/k B	120 K
	velocity, v 0		ε/m	158.03 m/s
	density, ρ 0		m/σ 3	1680.31 kg/m 3
	force, g 0 viscosity, µ 0	√	ε/σ mε/σ 2	4.86 × 10 -12 9.0345 × 10 -5	N kg/m/s
	conductivity, k 0	k B ε/m/σ 2	0.0188 W/m/K
	heat flux, q 0	ε 3 /m/σ 3	6.62 × 10 9	W/m 2

Table 2 :

 2 Normalization parameters for Figs.[START_REF] Fan | Force and heat current formulas for many-body potentials in molecular dynamics simulations with applications to thermal conductivity calculations[END_REF]

	(m) tot |/q 0

as a function of

ρ/ρ 0 g x /g 0 H/σ ūx /v 0 θb /T 0 |q

  Kn + d 2 Kn 2 + d 3 Kn 3 ) (12σ µ Kn + 1)

		1.3	1.15 1.15
		1.2	1.1 1.1
			1.05 1.05
	b / θ -	1.1	H/2-2λ H/2-2λ	z H/2-λ H/2-λ	H/2 H/2
	θ	1
		0.9			g x /g 0 = 0.02 g x /g 0 = 0.01 g x /g 0 = 0.005 g x /g 0 = 0.0025
					Eq. (11)
		+ a 2 Kn 2 2304(d 0 + d 1 0.8 D D = -1/2 -1/4 0 1/4 1/2 z / H	, B =	96σ µ Kn + 8 D	, C =	-16 D	,

Table 3 :

 3 Dimensionless parameters for simulations presented in Figs.3 and 8; σ µ and ξ are used in the analytical functions (

						summarizes
	ρ/ρ 0	Kn	Je	Br	σ µ	ξ
	0.6 0.0075 9.6 0.13 2.63 36
	0.6	0.036 5.1 0.24 2.46 36
	0.1	0.036 5.5 0.21 4.13 37
	0.05	0.036 5.2 0.15 5.72 25


Pierrick Coupé

Boris Mansencal

Michaël Clément

Rémi Giraud

Baudouin Denis De Senneville

Vinh-Thong Ta

Vincent Lepetit

José V Manjon

AssemblyNet: A large ensemble of CNNs for 3D Whole Brain MRI Segmentation

Whole brain segmentation of fine-grained structures using deep learning (DL) is a very challenging task since the number of anatomical labels is very high compared to the number of available training images. To address this problem, previous DL methods proposed to use a single convolution neural network (CNN) or few independent CNNs. In this paper, we present a novel ensemble method based on a large number of CNNs processing different overlapping brain areas. Inspired by parliamentary decision-making systems, we propose a framework called AssemblyNet, made of two "assemblies" of U-Nets. Such a parliamentary system is capable of dealing with complex decisions, unseen problem and reaching a relevant consensus. AssemblyNet introduces sharing of knowledge among neighboring U-Nets, an "amendment" procedure made by the second assembly at higher-resolution to refine the decision taken by the first one, and a final decision obtained by majority voting. During our validation, AssemblyNet showed competitive performance compared to state-of-the-art methods such as U-Net, Joint label fusion and SLANT. Moreover, we investigated the scan-rescan consistency and the robustness to disease effects of our method. These experiences demonstrated the reliability of AssemblyNet. Finally, we showed the interest of using semi-supervised learning to improve the performance of our method.

Introduction

Quantitative brain analysis is crucial to better understand the human brain and to analyze different brain pathologies. However, whole brain segmentation is still a very challenging problem, mostly due to the high number of anatomical labels compared to the limited number of available training data, especially when considering fine-grained segmentation. Manual segmentation of the whole brain is indeed a very tedious and difficult task, preventing the production of large annotated datasets.

To address this question, several methods have been proposed in the past years. By extending the single-atlas method paradigm, the multi-atlas framework has been successfully applied to whole brain segmentation [START_REF] Heckemann | Automatic anatomical brain MRI segmentation combining label propagation and decision fusion[END_REF], [START_REF] Avants | A reproducible evaluation of ANTs similarity metric performance in brain image registration[END_REF]. In such approaches, labeled templates are first nonlinearly registered to the target image. Afterwards, the estimated deformations are applied to the manual segmentations before fusing them. This type of methods efficiently deals with limited training data; however, the required multiple nonlinear registrations can result in a huge computational time. Moreover, regularization involved in registration may prevent to accurately capture complex local anatomical patterns.

To reduce the computational time of multi-atlas methods and to better capture local anatomy, patch-based methods have been introduced [START_REF] Coupé | Patchbased segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF]. In such approaches, the label fusion step is based on the nonlocal patch-based estimator. These methods demonstrated state-of-the-art performance for whole brain segmentation [START_REF] Wang | Multi-atlas segmentation with joint label fusion and corrective learning-an open source implementation[END_REF]- [START_REF] Asman | Non-local statistical label fusion for multi-atlas segmentation[END_REF] [START_REF] Asman | Hierarchical performance estimation in the statistical label fusion framework[END_REF]. One of the main references in the domain is the patch-based joint label fusion (JLF) which won the MICCAI challenge in 2012 [START_REF] Wang | Multi-atlas segmentation with joint label fusion and corrective learning-an open source implementation[END_REF] and which is still considered as the state of the art for fine-grained whole brain segmentation. In patch-based methods, usual machine learning such as sparse coding [START_REF] Tong | Segmentation of MR images via discriminative dictionary learning and sparse coding: Application to hippocampus labeling[END_REF] or neural networks [START_REF] Sanroma | Learning non-linear patch embeddings with neural networks for label fusion[END_REF] has been used in place of the nonlocal estimator. Recently, a fast framework has been proposed [START_REF] Giraud | An Optimized PatchMatch for multi-scale and multi-feature label fusion[END_REF] to further reduce the computational time required by patch-based methods.

More recently, deep leaning (DL) methods have also been proposed for 3D brain segmentation. Most of these methods were dedicated to coarse segmentation considering only few structures (e.g., <35 structures). For instance QuickNat [START_REF] Roy | QuickNAT: A fully convolutional network for quick and accurate segmentation of neuroanatomy[END_REF] works on 27 structures, Bayesian QuickNat works on 33 structures [START_REF] Roy | Bayesian QuickNAT: Model uncertainty in deep whole-brain segmentation for structure-wise quality control[END_REF], 3DQ works on 28 structures [START_REF] Paschali | 3DQ: Compact Quantized Neural Networks for Volumetric Whole Brain Segmentation[END_REF], DeepNat works on 25 structures [START_REF] Wachinger | DeepNAT: Deep convolutional neural network for segmenting neuroanatomy[END_REF], the method proposed in [START_REF] Roy | Concurrent Spatial and Channel 'Squeeze & Excitation' in Fully Convolutional Networks[END_REF] works on 27 structures and the approach in [START_REF] Rickmann | Recalibrating 3D ConvNets with Project Excite[END_REF] works on 32 structures. The problem of whole brain segmentation considering fine-grained structures (i.e. >100 structures) is much more complex. Consequently, less works have been dedicated to this problem [START_REF] Henschel | FastSurfer-A fast and accurate deep learning based neuroimaging pipeline[END_REF]- [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF]. Moreover, most of the proposed methods were based on 2D frameworks. In fact, due to limited GPU memory, first attempts were based on patchwise strategies [START_REF] Wachinger | DeepNAT: Deep convolutional neural network for segmenting neuroanatomy[END_REF], [START_REF] De Brebisson | Deep Neural Networks for Anatomical Brain Segmentation[END_REF], [START_REF] Li | On the Compactness, Efficiency, and Representation of 3D Convolutional Networks: Brain Parcellation as a Pretext Task[END_REF] or 2D segmentation (slice by slice) [START_REF] Roy | QuickNAT: A fully convolutional network for quick and accurate segmentation of neuroanatomy[END_REF], [START_REF] Henschel | FastSurfer-A fast and accurate deep learning based neuroimaging pipeline[END_REF], [START_REF] Roy | Error Corrective Boosting for Learning Fully Convolutional Networks with Limited Data[END_REF]. Only recently, 3D fully convolutional network methods were proposed using reduced input data size (i.e., 128×128×128 voxels) [START_REF] Wong | 3D Segmentation with Exponential Logarithmic Loss for Highly Unbalanced Object Sizes[END_REF] or using an Spatially Localized Atlas Network Tiles (SLANT) strategy [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF]. This latter framework divides the whole volume into overlapping sub-volumes, each one being processed by a different U-Net [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] (e.g., 8 or 27). The ensemble SLANT strategy addresses the problem of limited GPU memory and simplifies the complex problem of fine-grained whole brain segmentation into simpler problems, better suited to limited training data.

In this paper, we present a new method able to deal with fine-grained whole brain segmentation at full resolution and based on 3D convolution neural networks (CNNs). To this end, we propose to extend the SLANT framework by using a much larger number of more compact 3D U-Nets (i.e., from 27 to 250) while keeping processing time similar. The main question to address is the optimal organization of this large ensemble of CNNs. To this end, we propose a new framework we call AssemblyNet. Inspired by the decision-making process developed by human societies to deal with complex problems, we decided to model a parliamentary system based on two separate assemblies. Such bicameral -meaning two chambers -parliament has been adopted by many countries around the world. A bicameral system is usually composed of an upper and a lower chamber, both having their own independency to ensure the balance of power. However, an assembly may communicate its vote to the other for amendment. Such parliamentary system is capable of dealing with complex decisions, unseen problem and reaching a relevant consensus. This study extends our conference paper [START_REF] Coupé | AssemblyNet: A Novel Deep Decision-Making Process for Whole Brain MRI Segmentation[END_REF] with more complete experiments investigating i) the impact of semi-supervised learning ii) the scan-rescan reliability of our method and iii) the robustness to disease effects of the proposed AssemblyNet. Moreover, we added additional ablation study and a new quality metric. Compared to [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF], our contributions are: i) the use of prior knowledge based on fast atlas registration, ii) a knowledge sharing between CNNs using nearest neighbor transfer learning, iii) iterative refinement process based on a multiscale cascade of assemblies and iv) the use of student-teacher semi-supervised learning based on a well-designed auxiliary dataset.

Materials and Methods

Method overview

In AssemblyNet, both assemblies are composed of 3D U-Nets considered as "assembly members" (see Figure 1). Each member represents one territory (i.e., brain area) in the final vote. To this end, we used spatially localized networks where each U-Net only processes a sub-volume of the global volume, as done in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF]. Sub-volumes overlap each other, so the final segmentation results from an overcomplete aggregation of local votes. A majority vote is used to obtain the global segmentation. Moreover, each member can share knowledge with their nearest neighbor in the assembly. In particular, we propose a novel nearest neighbor transfer learning strategy, where weights of the spatially nearest U-Net are used to initialize the next U-Net.

In addition, we also propose to use prior knowledge on the expected final decision which can be viewed as the bill (i.e., draft law) submitted to an assembly for consideration. As prior knowledge, we decided to use nonlinearly registered Atlas prior.

Finally, we also propose modeling communication between both assemblies using an innovative multiscale strategy. In AssemblyNet, we use a multiscale cascade of assemblies where the first assembly produces a coarse decision at 2×2×2 mm 3 . This coarse decision is then transmitted to the second assembly for analysis at 1×1×1 mm 3 . This amendment procedure is similar to an error correction or a refinement step. After consideration by both assemblies, the bill under consideration becomes a law which represents the final segmentation in our system.

Datasets

Training dataset: 45 T1w MRI from the OASIS dataset [START_REF] Marcus | Open Access Series of Imaging Studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF] manually labeled according to the BrainCOLOR protocol were used for training. The selected images were the same than the ones used in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF] . The age range was 18-96y for this dataset. This dataset, as provided by Neuromorphometrics, included several pre-processing steps. First, all the MRI scans were corrected for bias field inhomogeneity using [START_REF] Arnold | Qualitative and Quantitative Evaluation of Six Algorithms for Correcting Intensity Nonuniformity Effects[END_REF]. Second, all the scans were registered along the anterior commissure (AC) and the posterior commissure (PC) using anatomical landmarks. Therefore, the original space of images was not the native space. The OASIS scan resolution was 1×1×1 mm3 .

Testing dataset: 19 T1w MRI manually labeled according to the BrainCOLOR protocol were used for testing. These MRI came from three different datasets: 5 from the OASIS dataset, the Colin27 atlas [START_REF] Collins | Design and construction of a realistic digital brain phantom[END_REF] and 13 from the CANDI database [START_REF] Kennedy | CANDIShare: A Resource for Pediatric Neuroimaging Data[END_REF]. This testing dataset is the same one used in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF]. The age range was 20-89y for OASIS, 5-15y for CANDI and the age was 27y for Colin27. The OASIS and CANDI scans as provided by Neuromorphometrics included pre-processing (inhomogeneity correction; AC/PC registration). The resolution of CANDI scans was 0.94×1.5×0.94mm 3 . The Colin27 atlas was in the MNI space at 0.5×0.5×0.5 mm 3 (its original space) and included inhomogeneity correction, intensity normalization and averaging of multiple acquisitions.

Scan-Rescan dataset: 8 T1w MRI (from 4 subjects) manually labeled according to the BrainCOLOR protocol were used for scan-rescan experiment. The same expert segmented both the scan and the rescan image. These MRI came from 2 different datasets: 3 from the OASIS dataset (not included in the training) and one from a patient with Alzheimer's disease from the ADNI dataset [START_REF] Weiner | The Alzheimer's Disease Neuroimaging Initiative: A review of papers published since its inception[END_REF]. The OASIS included preprocessing and were in the AC/PC space while ADNI scans were in the native space at 1.2×0.94×0.94 mm 3 . Pathological dataset: 29 T1w MRI manually labeled according to the BrainCOLOR protocol were included in the pathological dataset. We did not use the rescan image of the AD patient described in the scan-rescan dataset. These MRI came from the ADNI dataset. There were 15 cognitively control subjects and 14 patients with Alzheimer's disease in this dataset. The age of the subjects varied from 62y to 88y. The ADNI scans did not include pre-processing and thus were in their native space.

During our experiments, we used the 132 anatomical labels consistent across subjects (see Tab 1 in supplementary material for the list). In addition to Neuromorphometrics datasets, we used external MRIs from several open access datasets for semisupervised learning.

Lifespan dataset: 360 unlabeled T1w MRI were randomly selected under constraints from the dataset used in our previous BigData studies [START_REF] Coupé | Lifespan Changes of the Human Brain In Alzheimer's Disease[END_REF], [START_REF] Coupé | Towards a unified analysis of brain maturation and aging across the entire lifespan: A MRI analysis[END_REF] to build the lifespan dataset. This dataset was based on 9 datasets publicly available (C-MIND1 , NDAR2 , ABIDE 3 , ICBM4 , IXI5 , OASIS6 , AIBL7 , ADNI1 and ADNI28 ). From 1y to 90y, we selected 2 females and 2 males for each age (i.e., 2F and 2M of 1 year old, 2F and 2M of 2 years old and so on). Therefore, we obtained a balanced group with 50% of each gender uniformly distributed from 1y to 90y. We made sure that none of the training or testing subjects were selected in this auxiliary dataset.

AssemblyNet framework

Preprocessing: To homogenize input orientations and intensities, all the images were first preprocessed using the volBrain pipeline9 [START_REF] Manjón | volBrain: An online MRI brain volumetry system[END_REF] with the following steps: i) denoising [START_REF] Manjón | Adaptive nonlocal means denoising of MR images with spatially varying noise levels[END_REF], ii) inhomogeneity correction [START_REF] Tustison | N4ITK: Improved N3 Bias Correction[END_REF], iii) affine registration into the MNI space (181×217×181 voxels at 1×1×1 mm 3 ) [START_REF] Avants | A reproducible evaluation of ANTs similarity metric performance in brain image registration[END_REF], iv) tissue-based intensity normalization [START_REF] Manjón | Robust MRI brain tissue parameter estimation by multistage outlier rejection[END_REF] and v) brain extraction [START_REF] Manjón | Nonlocal intracranial cavity extraction[END_REF]. Finally, image intensities were centralized and normalized within the brain mask and the background was set to zero.

Atlas priors:

To obtain priors knowledge on the expected results, we performed a nonlinear registration of an Atlas (based on the 45 training images) to the subject under consideration. To construct this atlas, we non-linearly registered the training cases to a reference -the MNI152 T1 template -using ANTS software [START_REF] Avants | A reproducible evaluation of ANTs similarity metric performance in brain image registration[END_REF]. The estimated nonlinear deformations were then applied to the manual segmentations. Finally, the warped manual segmentations were fused to construct the atlas labels by using a majority vote. The non-linear registration of the Atlas prior to the subject under study was done with an unsupervised deep learning framework similar to VoxelMorph trained on the lifespan dataset [START_REF] Balakrishnan | VoxelMorph: A Learning Framework for Deformable Medical Image Registration[END_REF]. There are three main differences between VoxelMorph and our non-linear registration network. First, our network works on deformation fields at 2 mm resolution, which are internally interpolated to 1 mm resolution by the network. This step assumes that the deformation field is spatially smooth. This approach has two 2 benefits, the interpolation process imposes an intrinsic regularization and the smaller size of the volume enables to increase the number of filters at each convolutional layer level giving us almost ten times more learnable parameters (2,165,955 compared to the 259,675 for VoxelMorph). Second, our network is trained using non-skull-striped images contrary to VoxelMorph that requires this preprocessing step. Finally, the used loss function is based on the mean absolute error of both image intensities and labels (last version of VoxelMorph can use also intensities and labels simultaneously, but this last version has been never released as far as we know).

Assembly description: Each assembly was composed of 125 3D U-Nets equally distributed in the MNI space along each axis (i.e., 5 along x, y and z). In the following, we use U(x, y, z) " x, y, z Î [1..5] as the position of the U-Net in the assembly. We experimentally found that 5×5×5 U-Nets produced the best compromise between segmentation accuracy and computational time (see Results section). Each 3D U-Net processed a sub-volume large enough to ensure at least 50% of overlap between subspaces. At the end, a majority vote was used to aggregate the local votes.

Nearest neighbor transfer learning:

To enable knowledge sharing between U-Nets within an assembly, we propose a new transfer learning where the weights of a trained U-Net are used to initialize nearest U-Nets in the assembly (i.e., U-Nets processing overlapping sub-volumes and thus dealing with similar anatomy). At the beginning, we trained a first U-Net from scratch in an image corner (e.g., the red U-Net at position U(1,1,1) in Figure 1). The weights of this first trained U-Net were copied to the next U-Net on the same line (e.g., the green U-Net at position U(2,1,1) in Figure 1). This second U-Net was then used to initialize the next network on the same line and so on. Once the first line was trained, each U-Net of the second line U(x,2,1) was initialized with the U-Net at the same position on the previous line U(x,1,1) and so on (e.g., the orange U-Net was initialized with weights of the red U-Net). Finally, Finally, once all the U-Nets on the same plan U(x,y,1) were trained, each U-Net of the next plan U(x,y,2) was initialized with the U-Net at the same position on the previous plan U(x,y,1) and so on. During the transfer learning, we copied only the weights of the descending/contraction path of the U-Net architecture.

Multiscale cascade of assemblies:

To make our decision-making system faster and more robust, we decided to use a multiscale framework. Consequently, the first assembly at 2×2×2mm 3 produced a coarse segmentation. Afterwards, an up-sampling to 1×1×1mm 3 of this segmentation was performed using nearest neighbor interpolation. The up-sampled segmentation at 1×1×1mm 3 of the first assembly was then used as an additional input in the second assembly. Consequently, the second assembly had three 3D sub-volumes as input (i.e., T1w, Atlas priors and up-sampled coarse segmentation all at 1×1×1mm 3 ).

Ensemble framework

Ensemble is a well-known paradigm in machine learning that is used to improve the global performance of a method. This improvement is obtained by training several models before fusing them. Over the past decades, ensemble learning has been extensively studied, specially to deal with small sample size and complex problems [START_REF] Wu | Brain segmentation based on multi-atlas guided 3D fully convolutional network ensembles[END_REF]. In brain segmentation, most of the recent DL methods are based on multi-view ensemble [START_REF] Roy | QuickNAT: A fully convolutional network for quick and accurate segmentation of neuroanatomy[END_REF], [START_REF] Henschel | FastSurfer-A fast and accurate deep learning based neuroimaging pipeline[END_REF], [START_REF] Mehta | BrainSegNet: a convolutional neural network architecture for automated segmentation of human brain structures[END_REF]. In such framework, 2D models along axial, coronal and sagittal view are trained before fusing their outputs to enforce 3D consistency. A variant consists in fusing 2D models and 3D models [START_REF] Mehta | BrainSegNet: a convolutional neural network architecture for automated segmentation of human brain structures[END_REF], [START_REF] Zheng | A new ensemble learning framework for 3D biomedical image segmentation[END_REF]. It is also possible to fuse predictions of models working at different scales [START_REF] Moeskops | Automatic Segmentation of MR Brain Images With a Convolutional Neural Network[END_REF]. Moreover, there exist frameworks based on aggregation of predictions from multiple models with different architectures [START_REF] Kamnitsas | Ensembles of Multiple Models and Architectures for Robust Brain Tumour Segmentation[END_REF] or trained with different subset of the training dataset [START_REF] Dolz | Deep CNN ensembles and suggestive annotations for infant brain MRI segmentation[END_REF]. In addition, dropout has been proposed to generate several instances of a model at test time to reduce over-fitting [START_REF] Wu | Brain segmentation based on multi-atlas guided 3D fully convolutional network ensembles[END_REF]. However, all these methods are based on few models (typically <10 networks). The investigation of using a larger number of networks is recent in DL for brain segmentation. In [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF], SLANT strategy proposes to train more models (>10) on different areas of the brain that enforces model diversity before fusing them. The use of a large ensemble of CNNs has yielded state-of-the-art results for fine-grained whole brain segmentation.

In this paper, the proposed AssemblyNet includes almost all these strategies in a single framework. First, AssemblyNet includes a multiscale cascade of assemblies that allows to take advantage of models trained at different resolutions. Second, each U-Net is trained with a different set for training and validation to take advantage of all the available training data. Third, each U-Net is trained on a different but overlapping brain area that ensures models diversity as in SLANT. In addition, we use dropout at test time to simulate several instances of each U-Net that reduces over-fitting. We also performed temporal averaging of model weights based on snapshot ensembles [START_REF] Laine | Temporal Ensembling for Semi-Supervised Learning[END_REF]. Finally, compared to [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF], we propose to use a much larger number (from 27 to 250) of more compact U-Nets to better deal with limited training data.

Semi-supervised learning framework

In this study, we also investigated the use of semi-supervised learning (SSL) to further improve segmentation accuracy. SSL aims at using a small amount of labeled data in combination with a large number of unlabeled images to achieve higher performance. In medical image analysis, these techniques are particularly interesting to overcome the limited amount of training data and the complexity of the labeling process.

Over the past years, several strategies have been proposed to make SSL efficient within DL framework mainly based on the teacher-student paradigm. Within this paradigm, SSL methods can integrate a consistency term on the predictions of unlabeled samples to force the student network not to diverge [START_REF] Laine | Temporal Ensembling for Semi-Supervised Learning[END_REF], [START_REF] Tarvainen | Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results[END_REF]. This idea is usually based on a "Mean teacher" by using exponential moving averages over the predictions or the model parameters at different steps during training [START_REF] Huang | Snapshot Ensembles: Train 1, get M for free[END_REF]. Another teacher-student strategy consists in using the available labeled data to generate weak labels (also called pseudo labels) for unlabeled examples using a teacher model. This weak labels are in turn used as training dataset to improve the robustness of the student network [START_REF] Yalniz | Billion-scale semi-supervised learning for image classification[END_REF]. Afterwards, the student network is fine-tuned on the original labeled data to avoid error propagation. SSL strategies have been successfully applied to improve segmentation results in different medical applications [START_REF] Zhou | Collaborative Learning of Semi-Supervised Segmentation and Classification for Medical Images[END_REF]- [START_REF] Chen | Multitask Attention-Based Semi-supervised Learning for Medical Image Segmentation[END_REF]. In our context of whole brain segmentation, the authors of [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF], [START_REF] Roy | Error Corrective Boosting for Learning Fully Convolutional Networks with Limited Data[END_REF] proposed to use auxiliary datasets segmented with traditional tools such as Freesurfer [START_REF] Fischl | FreeSurfer[END_REF] or non-local spatial staple label fusion (NLSS) [START_REF] Asman | Non-local statistical label fusion for multi-atlas segmentation[END_REF] to improve their segmentation framework based on deep learning. However, such methods can require impractical computational burden (e.g., 21 CPU years in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF]) and classical tools may provide suboptimal auxiliary segmentations.

Here, we take inspiration of the teacher-student paradigm from [START_REF] Yalniz | Billion-scale semi-supervised learning for image classification[END_REF] to leverage the fast processing capabilities and high segmentation accuracy of the proposed AssemblyNet. In our SSL framework, an AssemblyNet teacher -trained on the 45 training images -was first used to segment unlabeled images (i.e., the 360 images of the lifespan dataset). Then, these 360 pseudo-labeled images were used to train from scratch an AssemblyNet student. At the end, the AssemblyNet student was fine-tuned on the 45 manual segmentations of the training dataset. As shown in [START_REF] Yalniz | Billion-scale semi-supervised learning for image classification[END_REF], this finetuning step is able to limit error propagation within SSL framework. During our experiments, we investigated the iteration of this procedure considering that the obtained AssemblyNet student could be a good teacher for a second student generation. In our SSL framework, we took care to build the unlabeled images dataset balanced in age and gender in order to limit bias introduction in the pseudo-labeled population.

Implementation details

Data augmentation: First, the images of the training and lifespan datasets were flipped along mid sagittal plane in the MNI space. Then, we used MixUp data augmentation during training to minimize overfitting problems [START_REF] Zhang | mixup: Beyond Empirical Risk Minimization[END_REF]. This method performs a linear interpolation of a random pair of training examples and their corresponding labels.

Training framework: For all the networks, we used the 3D U-Net architecture proposed in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF], but with a lower number of filters. Instead of using a basis of 32 filters of 3×3×3 -32 for the first layer, 64 for the second and so on -we selected a basis of 24 filters of 3×3×3 to reduce by 25% the network size. We experimentally found that this setting reduced memory consumption without impacting performance. The used architecture is presented in Figure 2. Each block was composed of batch normalization, convolution and ReLU activation. The skip connections between encoder and decoder were based on concatenation. In addition, dropout was done between each level of the encoder. Moreover, upsampling in the decoder was based on trilinear interpolation. Finally, a SoftMax was done before performing argmax to obtain the final label for each voxel.

For all the U-Nets, we used the same parameters: batch size = 1, optimizer = Adam, epoch = 100, loss = Dice and dropout = 0.5 after each block of the descending path. For the U-Nets of the first assembly at 2×2×2 mm 3 , we used input resolution = 32×48×32 voxels and input channel = 2 (i.e., T1w and Atlas priors). For the U-Nets of the second assembly at 1×1×1 mm 3 , we used input resolution = 64×72×64 voxels and input channel = 3 (i.e., T1w, Atlas priors and up-sampled coarse segmentation). In addition, to compensate for the small batch size, we performed temporal averaging of model weights based on snapshot ensembles [START_REF] Laine | Temporal Ensembling for Semi-Supervised Learning[END_REF]. At the end of the 100 epochs, we performed additional 20 epochs where the model estimated at each epoch is averaged with previous ones using a moving average. Such average of model weights along the optimization trajectory leads to better generalization than usual training [START_REF] Izmailov | Averaging Weights Leads to Wider Optima and Better Generalization[END_REF]. For the SSL step, we used only 20 epochs for normal optimization and 10 epochs for moving average. Finally, we also performed dropout at test time [START_REF] Gal | A Theoretically Grounded Application of Dropout in Recurrent Neural Networks[END_REF]. For each U-Net, we generated 3 different outputs before averaging them (with dropout layer active). Such method helps reducing variance of the networks. As in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF], the experiments were done with an NVIDIA Titan Xp with 12 GB memory and thus processing times are comparable.

Computational time:

The preprocessing steps take around 90 seconds. The nonlinear registration of the atlas takes less than 5 seconds. The first assembly at 2×2×2 mm 3 requires 3 minutes to segment an image while the second assembly at 1×1×1 mm 3 requires 5 minutes. At the end, the final segmentation is registered back to the original space using the inverse affine transform estimated during preprocessing. This interpolation takes around 30 seconds. Therefore, the full AssemblyNet process takes around 10 minutes including preprocessing, segmentation, and inverse registration back to the original space.

Validation framework

First, for each testing subject, we estimated the average Dice coefficient on the 132 considered anatomical labels (without background) in the original space. Afterwards, we estimated the global mean Dice in % over the 19 images of the testing dataset. In this experiment, we compared AssemblyNet with several state-of-the-art methods. First, the patch-based joint label fusion (JLF) [START_REF] Wang | Multi-atlas segmentation with joint label fusion and corrective learning-an open source implementation[END_REF] was used as reference. In addition, we included U-Net [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF], SLANT-8 and SLANT-27 methods as proposed in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF]. SLANT-8 is based on 8 U-Nets processing non-overlapping sub-volumes of 86×110×78 voxels while SLANT-27 is based on 27 U-Nets processing overlapping sub-volumes 96×128×88 voxels. All these methods were trained on the same 45 training images described in the section datasets. Moreover, all these methods used the following pre-processing steps: i) affine registration to MNI space using NiftyReg [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF], ii) bias filed correction using N4 [START_REF] Tustison | N4ITK: Improved N3 Bias Correction[END_REF] and iii) intensity harmonization using regression-based normalization. In addition, we included SLANT-27 FT trained on 5111 auxiliary images segmented using NLSS [START_REF] Asman | Hierarchical performance estimation in the statistical label fusion framework[END_REF] and fined tuned on the 45 training images. These are the best published results for fine-grained whole brain segmentation to our knowledge. Afterwards, to show the impact of affine registration to MNI space, we also presented the results of "naïve U-Net" working directly in their original space. For all these methods, we report the results published in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF]. We also used the docker implementation of SLANT-27 10 to produce segmentations of all the considered datasets. Finally, we also compared mean surface distance between AssemblyNet, AssemblyNet SSL (with semi-supervised learning) and SLANT-27 FT.

For the scan-rescan reliability experiment, we rigidly registered the re-scan image into the original space of the scan images [START_REF] Avants | A reproducible evaluation of ANTs similarity metric performance in brain image registration[END_REF]. We then interpolated the re-scan segmentations into the original space of the scan images using the estimated transformation matrix. By this way, we estimated the Dice coefficients between both manual segmentations (i.e., intra-rater consistency) and both automatic segmentations (i.e., intra-method consistency). Moreover, we estimated the methodexpert consistency as the Dice coefficients between the automatic segmentation of the rescan images and the manual segmentation of the scan image.

For the experiment on the robustness to disease effects, we first computed the Dice on the 29 ADNI subjects. We then compared the Dice coefficients obtained for cognitively normal (CN) subjects and patients with Alzheimer's Disease (AD) to study the impact of pathology on segmentation accuracy.

For all these experiments, we used one-sided non-parametric Wilcoxon signed-rank test at 95% of confidence to assess the significance of Dice improvement as in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF]. Moreover, we used one-sided Mann-Whitney rank test at 95% of confidence to assess the significance of Dice decrease between the AD group compared to the CN group.

Results

AssemblyNet performance

First, we evaluated the proposed contributions (see Table 1). Compared to baseline results at 2×2×2 mm 3 (Dice=67.4%), the use of Atlas prior provided a gain of 0.3 percentage point (pp) in terms of mean Dice. Moreover, the combination of Atlas prior and transfer learning improved by 0.5 pp the baseline mean Dice. In addition, multiscale cascade of assemblies increased by 1.1 pp the mean Dice obtained with Assembly at 1×1×1 mm 3 without multiscale cascade (Dice=72.2%). Finally, AssemblyNet outperformed by 5.9 pp the mean Dice obtained with baseline Assembly at 2×2×2 mm 3 . The Dice coefficients produced by AssemblyNet were significantly better than the Dice coefficients produced by all the considered alternatives. Note these results were obtained using only 45 training cases.

Table 1: Evaluation of the proposed contributions. The mean Dice (std) is evaluated on the 19 images of the test dataset in the original space for the 132 considered labels (without background). Testing time includes image preprocessing and registration back to the original space. * indicates a significant lower Dice compared to AssemblyNet using a Wilcoxon test.

Methods

Atlas prior

Transfer learning

Multiscale

Dice in % (std)

Training time

Testing time

Assembly at 2×2×2 mm Afterwards, we studied the impact of the number of U-nets on the performance of the Assembly at 2×2×2 mm 3 with atlas prior and transfer learning (see Table 2). During this experience, the accuracy reached a plateau from 125 U-Nets. Using more networks did not provide additional improvements while increasing computational time. Therefore, in the following, 125 was used as the default number of networks in each assembly.

Impact of semi-supervised learning

Second, we evaluated the impact of the proposed teacher-student semi-supervised learning (SSL) framework (see Table 3). The used of the lifespan dataset -labeled with teacher -in the training of the student lead to an improvement of 0.6pp of mean Dice compared to the mean Dice of the teacher (Dice=73.3%). The fine-tuning (FT) step further increased the mean Dice by 0.3 pp. These results are in line with previous literature [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF], [START_REF] Yalniz | Billion-scale semi-supervised learning for image classification[END_REF] on the role of the FT step to limit error propagation within semisupervised learning framework. Finally, the second iteration after FT produced marginal improvement and lead to a mean Dice of 74.0%. This improvement was not significant. Therefore, in the following, we used the first student generation since the time required for the second iteration is not justified by the performance improvement.

Comparison with state-of-the-art methods

We compared AssemblyNet with state-of-the-art methods (see ). This number of 21 CPU years is reported in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF].

In addition, compared to SLANT-27 FT, AssemblyNet provided better results without library extension while being faster to train and to execute. Using our SSL framework based on 360+45 images, AssemblyNet SSL obtained a gain of 1pp compared to SLANT-27 FT trained over 5111+45 images. According to [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF], their library extension required 21 CPU years to be completed. Consequently, such an approach is impractical or very costly using a cloud-based solution. The proposed SSL framework is more practical in term of time and resources. Finally, our AssemblyNet SSL was significantly better than AssemblyNet and SLANT-27 FT (docker). In our framework, SLANT-27 FT (docker) obtained slightly lower results than the ones published in the original paper [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF]. This may come from hardware and environment differences.

In addition, we analyzed the performance of the methods according to the dataset. Mean Dice coefficients obtained on each testing dataset (i.e., OASIS, CANDI and Colin27) are provided in As expected, all the methods performed better on adult scans from the OASIS dataset since the training dataset comes from the same cohort. Moreover, all the images were acquired with the same protocol on the same scanner and provided in the same space. First, we can note the good performance of JLF compared to U-Net and SLANT-8. Moreover, when considering methods trained with 45 training, AssemblyNet outperformed U-Net by 8.2 pp, JLF by 4.4 pp and SLANT-27 by 1.2 pp of mean Dice. When considering all the methods, AssemblyNet SSL obtained significantly better Dice than SLANT-27 FT (docker). It has to be noted that in [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF], the authors have shown that the SLANT-27 FT significantly outperformed U-Net, JLF and SLANT-8. Finally, on the OASIS images, using SSL did not significantly improve the AssemblyNet results.

On child scans from the CANDI dataset acquired with a different protocol, we can first note a dramatic drop in performance for all the methods except for AssemblyNet, AssemblyNet SSL and SLANT-27 FT. Moreover, when considering methods trained with 45 training, AssemblyNet outperformed U-Net by 19.7 pp, JLF by 19.2 pp and SLANT-27 by 9 pp of mean Dice. When considering all the methods, AssemblyNet SSL obtained significantly better Dice than AssemblyNet and SLANT-27 FT (docker).

On the high-resolution Colin27 image, we also observed an important decrease of performance for all the methods except for AssemblyNet, AssemblyNet SSL and SLANT-27 FT. As for CANDI dataset, AssemblyNet obtained the best segmentation accuracy with or without SSL on this dataset.

The comparison between Naïve U-Net (working in the original space) and U-Net (working in the MNI space) showed that performing an affine registration to MNI produced a gain of 16 pp (see Table 4). On the OASIS dataset (in the AC/PC space), the training and testing spaces were similar and thus the naïve U-Net obtained descent results (see Table 5). However, on the Colin27 atlas in the MNI space at 0.5×0.5×0.5mm 3 , the Naïve U-Net obtained a Dice=0 since this Atlas is not in the training space. While slower, using an additional affine registration to MNI space allowed to improve performance and to be robust to image space and resolution.

Finally, we compared mean surface distance (MSD) from manual segmentations to automatic segmentations on AssemblyNet, AssemblyNet SSL and SLANT-27 FT (see As an illustration, Figure 3 shows the segmentations of the central slides in the original space obtained by SLANT-27 FT (docker) and AssemblyNet SSL on the first subject of the testing dataset (ID=1120_3). Both methods provided good segmentations although AssemblyNet SSL segmentation was less smooth especially around sulci (e.g., cerebellum -see red ellipses). Moreover, we can observe an over segmentation of cortical gray matter in SLANT-27 FT segmentation as visible in the error map where structures appeared (see green ellipses). Finally, this figure shows the staircasing artifacts present in the human segmentation (e.g., pallidum -see the pink ellipses) while automatic methods were more regular and consistent.

Scan-rescan consistency

The study of segmentation reproducibility produced by a segmentation method is also highly important especially in medical imaging. Therefore, we carried out a scanrescan experiment to investigate the consistency and reliability of the proposed method. The results on the 4 images of the scan-rescan dataset are provided in Table 7. It has to be highlighted that there are variations between both acquisitions due to patient's motion, distortion, inhomogeneity and noise. Therefore, agreement between both acquisitions is not expected to be perfect but higher Dice indicates better method stability, consistency and reliability. First, we estimated the consistency of the segmentations provided by SLANT-27 FT (docker), AssemblyNet, AssemblyNet SSL and the expert on the scan and the rescan images. As expected, we can note that automatic methods were much more consistent (>90% of Dice) than the expert who obtained a scan-rescan segmentation consistency of 76.8% of Dice. Moreover, we can see that AssemblyNet was more consistent than SLANT-27 FT especially when using the proposed teacher-student SSL framework. The method consistency of AssemblyNet SSL was significantly higher than the SLANT-27 FT (docker) and AssemblyNet.

In addition, we estimated the Expert-Method consistency as the mean Dice coefficients between the automatic segmentation on the rescan image and the manual segmentation of the scan image. The Expert-Method consistency of AssemblyNet SSL was significantly higher than the consistency of SLANT-27 FT (docker) but not than the consistency of AssemblyNet. Finally, the Expert-Method consistencies obtained by automatic methods were not significantly lower than the intra-expert consistency although this difference was almost significant for SLANT-27 FT docker (p=0.36 for AssemblyNet SSL, p=0.23 for AssemblyNet and p=0.07 for SLANT-27 FT).

Robustness to disease effects

The last part of our validation is dedicated to robustness to disease effects. To this end, we compared SLANT-27 FT (docker), AssemblyNet and AssemblyNet SSL on the pathological dataset composed of CN subjects and AD patients.

First, we estimated the mean Dice for both groups (see Table 8). For the CN group, AssemblyNet SSL obtained a significantly better Dice than both other methods. For the AD group, we obtained similar results although the improvement obtained by AssemblyNet SSL was higher for AD group (2.2 pp) than for the CN group (1.7 pp) compared to SLANT-27 FT (docker).

In addition, we compared the accuracy between CN and AD groups for the three methods. We found no significant differences between groups for all the methods although this difference was almost significant for SLANT-27 FT docker (p=0.31 for AssemblyNet SSL, p=0.18 for AssemblyNet and p=0.06 for SLANT-27 FT).

Finally, AssemblyNet obtained a global Dice 73.1% without SSL and 73.6% with SSL, while SLANT-27 FT obtained 71.6%. The results of AssemblyNet SSL were significantly better than the results obtained with both other methods.

Discussion

In this work, we presented a novel whole brain segmentation framework based on a large number of 3D CNN (i.e., 250 U-Nets) called AssemblyNet. First, we showed that the use of Atlas prior, nearest neighbor transfer learning and multiscale cascade of Assemblies enable to improve global segmentation accuracy. In further work alternative options could be investigated. First, atlas prior could be replaced by fast multi-atlas prior. Thanks to nonlinear registration methods based on deep learning, such prior is no more too expensive. Moreover, more advanced communication between assembly members should be investigated. Recent advances in multi-agent reinforcement learning seem a promising way [START_REF] Omidshafiei | Deep Decentralized Multitask Multi-Agent Reinforcement Learning under Partial Observability[END_REF]. Finally, in this paper, we focused only on the optimal organization of a large group of CNNs without studying the optimal assembly composition. Additional works should investigate this point, for instance by introducing model diversity in the assembly.

Second, we studied the impact of the proposed SSL based on a teacher-student paradigm. We showed that using few hundreds of well-balanced unlabeled data could significantly improve the results of AssemblyNet in all the cases (i.e., unseen acquisition protocol, age period and pathology). Compared to previous methods using larger auxiliary datasets labelled with classical tools [START_REF] Huo | 3D whole brain segmentation using spatially localized atlas network tiles[END_REF], [START_REF] Roy | Error Corrective Boosting for Learning Fully Convolutional Networks with Limited Data[END_REF], the proposed SSL framework is more practical in terms of computational time and resources. However, SSL is currently receiving special attention in deep learning community. Consequently new paradigms should be considered [START_REF] Ren | Meta-Learning for Semi-Supervised Few-Shot Classification[END_REF], [START_REF] Luo | Smooth Neighbors on Teacher Graphs for Semi-supervised Learning[END_REF].

Afterwards, we compared our AssemblyNet with state-of-the-art methods. We demonstrated the high performance of our method in terms of segmentation accuracy and computational time. First, these experiments demonstrated the advantage of using several CNNs to segment the whole brain since SLANT-27 and AssemblyNet clearly outperformed the use of a single U-Net. Moreover, these results showed that using a larger number of simpler CNNs within a multiscale framework is an efficient strategy.

While for OASIS we obtained Dice higher than intra-expert consistency, the accuracy on the CANDI dataset is still limited. This can come from several factors such as the lower image quality (e.g., more motion artifacts in child images) or the larger distance between adult training dataset and this child dataset. These points should be deeper investigated in future whole brain segmentation methods. In terms of computational time, our method could be further improved by using several GPUs. At training time, once the first U-Net is trained, several following U-Nets can be trained in parallel despite our transfer learning strategy. At testing time, AssemblyNet can be fully parallelized and thus the processing time could be drastically reduced using multiple GPUs.

In addition, we investigated the scan-rescan consistency of the proposed method. We showed that the intra-method consistency of our method reached 92.8% while intrarater consistency was limited to 76.8%. This result clearly demonstrates that our automatic method segments the whole brain in a more consistent manner than a human expert. Moreover, we show that the expert-method consistency obtained with our method is not significantly lower than the intra-rater consistency, which is an encouraging result. However, these results were obtained using only 4 scan-rescan subjects. In addition, the Neuromorphometrics dataset does not contain material to evaluate method reproducibility (same image segmented twice by the same expert). Such data would have been useful to evaluate if automatic methods had reached human variability. Finally, these results raise the question of using human segmentation as "gold standard" with a consistency lower than 80%. Semi-manual segmentations could be considered in order to reduce intra-rater variability.

Finally, we studied the robustness of AssemblyNet to pathology. To this end, we compared its accuracy on CN and AD groups. We observed a small but non-significant decrease of Dice for the AD group compared to CN group. Moreover, compared to SLANT-27, AssemblyNet was less impacted by the presence of the pathology. This is a first step towards a more extensive validation with other pathologies.

Conclusion

In this paper, we proposed to use a large number of CNNs to perform whole brain segmentation. We investigated how to organize this large ensemble of CNNs to accurately segment the brain. To this end, we designed a novel deep decision-making process called AssemblyNet based on two assemblies of 125 3D U-Nets. Our validation showed the very competitive results of AssemblyNet compared to state-ofthe-art methods. We also demonstrated that AssemblyNet is very efficient to deal with limited training data and to accurately achieve segmentation in a practical training and testing times. Finally, we demonstrated the interest of semi-supervised learning to improve the performance of our method on unseen acquisition protocol, age period and pathology.

Figure 1 :

 1 Figure 1: Illustration of the proposed AssemblyNet framework. Our method is based on two assemblies of 125 3D U-Nets integrated into a multiscale framework. The first assembly (in yellow) provides a coarse segmentation at 2×2×2 mm 3 . The second assembly (in purple) refines this coarse segmentation to produce the final segmentation at 1×1×1 mm 3 . Each 3D U-Net processes a different but overlapping area of the brain. The U-Nets in red in both Assemblies process the area indicated by red rectangles in the input images. The U-Nets in green process the area indicated by green rectangles. During training, the U-Nets in green and orange are initialized using the weights of the U-Nets in red by transfer learning. The output segmentations for each assembly are obtained by majority voting of the 125 3D U-Nets.

Figure 2 :

 2 Figure 2: Illustration of the used U-Net architecture. The number of input channels (NC) depends on the considered assembly (i.e., NC=2 at 2×2×2 mm 3 and NC=3 at 1×1×1 mm 3 ). Each block is composed of batch normalization (BN), convolution and ReLU activation. The number of 3x3x3 filters is indicated on the top of each block. The final output size depends on the number of structures (NS) in the considered sub-volume.

Figure 3 :

 3 Figure 3: Example of segmentations in the original space for the first testing subject (ID=1120_3). First rows: sagittal; axial and coronal views for the T1w MRI. Second row: manual segmentation produced by the expert. Third row: segmentation obtained by our AssemblyNet SSL. Fourth row: segmentation obtained by SLANT-27 FT (docker). Fifth row: binary difference between manual and AssemblyNet SSL segmentations. Last row: binary difference between manual and SLANT-27 FT (docker) segmentations. Colored ellipses indicate areas of interest

Table 2 :

 2 Impact of the number of U-nets on the Assembly at 2×2×2 mm 3 with atlas prior and transfer learning. The mean Dice (std) is evaluated on the 19 images of the test dataset in the original space for the 132 considered labels (without background). Testing time includes image preprocessing and registration back to the original space. * indicates a significant lower Dice compared to AssemblyNet based on 343 U-Nets using a Wilcoxon test.

	Methods	Number of U-Nets	Dice in % (std)	Training time	Testing time
	Assembly at 2×2×2 mm 3	27 (3×3×3)	66.1 (3.4)*	6h	3min
	Assembly at 2×2×2 mm 3	64 (4×4×4)	67.6 (3.4)*	15h	3min
	Assembly at 2×2×2 mm 3	125 (5×5×5)	67.9 (3.3)	29h	5min
	Assembly at 2×2×2 mm 3	216 (6×6×6)	67.9 (3.3)	2 days	7min
	Assembly at 2×2×2 mm 3	343 (7×7×7)	67.9 (3.3)	3 days	10min

Table 3 :

 3 Impact

	Methods	Training	Dice in %	Training	Library
		images	(std)	time	extension time
	AssemblyNet Teacher	45	73.3 (4.2)*	7 days	0s
	AssemblyNet SSL First student generation	360	73.6 (4.1)*	12 days	2.5 days
	AssemblyNet SSL First student generation	360 + FT 45	73.9 (4.0)	14 days	2.5 days
	AssemblyNet SSL Second student generation	360	73.9 (4.0)	26 days	5 days
	AssemblyNet SSL Second student generation	360 + FT 45	74.0 (3.9)	28 days	5 days

of the proposed semi-supervised learning framework on the 19 images of the testing dataset. The mean Dice (std) is evaluated on the 132 considered labels (without background) in the original space. * indicates a significant lower Dice compared to second generation of AssemblyNet SSL using a Wilcoxon test.

Table 4

 4

). When considering only methods trained with 45 images, AssemblyNet improved mean Dice obtained with U-Net and SLANT-8 by 16.3 pp, JLF by 9.9 pp and SLANT-27 by 7.2 pp. AssemblyNet was also efficient in terms of training and testing times compared to SLANT-based methods. It has to be noted that the Assembly at 2×2×2 mm outperformed all the methods using 45 training images (except AssemblyNet) while working at low resolution.

Table 4 :

 4 Comparison with state-of-the-art methods on the 19 images of the testing dataset. The mean Dice is evaluated on the 132 considered labels (without background) in the original space. * indicates a significant lower Dice compared to AssemblyNet SSL using a Wilcoxon test when compared to SLANT-27 FT (docker) AssemblyNet and the Assembly at 2×2×2 mm 3 .

	Methods	Training	Dice in %	Training	Testing	Library
		images	(std)	time	time	extension time
	Naïve U-Net [24]	45	41.0	33h	1min	0s
	U-Net [20]	45	57.0	33h	8min	0s
	SLANT-8 [20]	45	57.0	11 days	10min	0s
	JLF [4]	45	63.4	0s	34h	0s
	SLANT-27 [20]	45	66.1	42 days	15min	0s
	SLANT-27 FT [20]	5111 + FT 45	72.9	27 days	15min	21 years a
	Assembly at 2×2×2 mm 3	45	67.9 (3.3)*	29h	5min	0s
	SLANT-27 FT (docker)	5111 + FT 45	72.6 (2.8)*	27 days	15min	21 years a
	AssemblyNet	45	73.3 (4.2)*	7 days	10min	0s
	AssemblyNet SSL	360 + FT 45	73.9 (4.0)	14 days	10min	2.5 days

a Library extension time represents the CPU time required to segment 5111 MRI using NLSS (i.e., 34h×5111

Table 5

 5

(see Fig 1-4 in supplementary material for boxplots of Dice distributions).

Table 5 :

 5 Comparison with state-of-the-art methods using Dice on the different testing datasets (5 adult scans from OASIS, 13 child scans from CANDI child and the high-resolution Colin27 image based on scans average). Using Wilcoxon tests * indicates a significant lower Dice compared to AssemblyNet SSL when compared to SLANT-27 FT (docker) and AssemblyNet.

	Methods	Training	OASIS	CANDI	Colin27
		images	Dice in %	Dice in %	Dice in %
			(std)	(std)
	Naïve U-Net [24]	45	60.6 (0.6)	37.5 (4.3)	0.0
	U-Net [20]	45	70.6 (0.9)	51.4 (8.1)	62.1
	SLANT-8 [20]	45	69.9 (1.4)	51.9 (7.0)	59.7
	JLF [4]	45	74.6 (0.9)	59.0 (3.3)	64.6
	SLANT-27 [20]	45	76.6 (0.8)	62.1 (6.2)	66.5
	SLANT-27 FT [20]	5111 + FT 45	77.6 (1.2)	71.1 (2.3)	73.2
	SLANT-27 FT (docker)	5111 + FT 45	75.9 (1.7)*	71.3 (2.2)*	73.5
	AssemblyNet	45	78.8 (1.7)	71.1 (2.9)*	74.2
	AssemblyNet SSL	360 + FT 45	79.0 (2.0)	71.9 (2.9)	75.0

Table 6 and

 6  in supplementary material). Average mean surface distance showed similar trends than Dice scores. AssemblyNet SSL produced significant lower MSD in all the considered cases except for OASIS dataset compared to AssemblyNet.

Table 6 :

 6 Comparison with state-of-the-art methods using mean surface distance (MSD) on the different testing datasets (5 adult scans from OASIS, 13 child scans from CANDI child and the high-resolution Colin27 image based on scans average). Using Wilcoxon tests * indicates a significant greater MSD compared to AssemblyNet SSL.

	Methods	OASIS	CANDI	Colin27	Global
		MSD in mm	MSD in mm	MSD in mm	MSD in mm
	SLANT-27 FT (docker)	0.699 (0.076)*	1.062 (0.098)*	0.242	0.923 (0.243)*
	AssemblyNet	0.550 (0.058)	1.028 (0.153)*	0.212	0.859 (0.289)*
	AssemblyNet SSL	0.553 (0.067)	0.996 (0.134)	0.206	0.838 (0.270)

Table 7 :

 7 Reliability study on the scan-rescan datasets (3 adult scans from OASIS and one scan of a patient with AD from ADNI). The intra-method consistency is the mean Dice between the automatic segmentations obtained on the scan and rescan images. The Expert-Method consistency is the mean Dice between the automatic segmentation on the rescan image and the manual segmentation of the scan image. The Dice coefficients obtained on the scan and the rescan images are averaged. The intra-rater consistency is the mean Dice between the manual segmentations obtained on the scan and rescan images. Using Wilcoxon tests * indicates a significant lower Dice compared to AssemblyNet SSL.

	Methods	Training	Intra-method	Intra-rater	Expert-Method
		images	consistency	consistency	consistency
			Dice in % (std)	Dice in %	Dice in % (std)
				(std)
	SLANT-27 FT (docker)	5111 + FT 45	91.7 (1.5)*		72.9 (1.7)*
	AssemblyNet	45	92.0 (1.8)*	76.8 (4.1)	75.6 (1.7)
	AssemblyNet SSL	360 + FT 45	92.8 (1.8)		75.8 (1.9)

Table 8 :

 8 Methods comparison on the pathological dataset (29 scans from the ADNI dataset including 15 CN and 14 patients with AD). * indicates a significant lower Dice compared to AssemblyNet SSL using Wilcoxon tests.

	Methods	Training	CN	AD	ADNI
		images	Dice in %	Dice in %	Dice in %
			(std)	(std)	(std)
	SLANT-27 FT (docker)	5111 + FT 45	72.3 (1.6)*	71.0 (2.6)*	71.6 (2.2)*
	AssemblyNet	45	73.6 (1.6)*	72.6 (2.6)*	73.1 (2.2)*
	AssemblyNet SSL	360 + FT 45	74.0 (1.5)	73.2 (2.5)	73.6 (2.1)

https://research.cchmc.org/c-mind/

https://ndar.nih.gov

http://fcon_1000.projects.nitrc.org/indi/abide/

http://www.loni.usc.edu/ICBM/

http://brain-development.org/ixi-dataset/

http://www.oasis-brains.org

http://adni.loni.ucla.edu/research/protocols/mri-protocols

www.loni.usc.edu

https://www.volbrain.upv.es

https://github.com/MASILab/SLANTbrainSeg

Acknowledgements

This work benefited from the support of the project DeepvolBrain of the French National Research Agency (ANR-18-CE45-0013). This study was achieved within the context of the Laboratory of Excellence TRAIL ANR-10-LABX-57 for the BigDataBrain project. Moreover, we thank the Investments for the future Program IdEx Bordeaux (ANR-10-IDEX-03-02, HL-MRI Project), Cluster of excellence CPU and the CNRS/INSERM for the DeepMultiBrain project. This study has been also supported by the DPI2017-87743-R grant from the Spanish Ministerio de Economia, Industria Competitividad. The authors gratefully acknowledge the support of NVIDIA Corporation with their donation of the TITAN Xp GPU used in this research.

Moreover, this work is based on multiple samples. We wish to thank all investigators of these projects who collected these datasets and made them freely accessible. The C-MIND data used in the preparation of this article were obtained from the C-MIND Data Repository (accessed in Feb 2015) created by the C-MIND study of Normal Brain Development. This is a multisite, longitudinal study of typically developing children from ages newborn through young adulthood conducted by Cincinnati Children's Hospital Medical Center and UCLA and supported by the National Institute of Child Health and Human Development (Contract #s HHSN275200900018C). A listing of the participating sites and a complete listing of the study investigators can be found at https://research.cchmc.org/c-mind. The NDAR data used in the preparation of this manuscript were obtained from the NIH-supported National Database for Autism Research (NDAR). NDAR is a collaborative informatics system created by the National Institutes of Health to provide a national resource to support and accelerate research in autism. The NDAR dataset includes data from the NIH Pediatric MRI Data Repository created by the NIH MRI Study of Normal Brain Development. This is a multisite, longitudinal study of typically developing children from ages newborn through young adulthood conducted by the Brain Development Cooperative Group and supported by the National Institute of Child Health and Human Development, the National Institute on Drug Abuse, the National Institute of Mental Health, and the National Institute of Neurological Disorders and Stroke (Contract #s N01-HD02-3343, N01-MH9-0002, and N01-NS-9-2314, -2315, -2316, -2317, -2319 and -2320). A listing of the participating sites and a complete listing of the study investigators can be found at http://pediatricmri.nih.gov/nihpd/info/participating_centers.html. The ADNI data used in the preparation of this manuscript were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904). The ADNI is funded by the National Institute on Aging and the National Institute of Biomedical Imaging and Bioengineering and through generous contributions from the following: Abbott, AstraZeneca AB, Bayer Schering Pharma AG, Bristol-Myers Squibb, Eisai Global Clinical Development, Elan Corporation, Genentech, GE Healthcare, GlaxoSmithKline, Innogenetics NV, Johnson & Johnson, Eli Lilly and Co., Medpace, Inc., Merck and Co., Inc., Novartis AG, Pfizer Inc., F. Hoffmann-La Roche, Schering-Plough, Synarc Inc., as well as nonprofit partners, the Alzheimer's Association and Alzheimer's Drug Discovery Foundation, with participation from the U.S. Food and Drug Administration. Private sector contributions to the ADNI are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study was coordinated by the Alzheimer's Disease Cooperative Study at the University of California, San Diego. ADNI data are disseminated by the Laboratory for NeuroImaging at the University of California, Los Angeles. This research was also supported by NIH grants P30AG010129, K01 AG030514 and the Dana Foundation. The OASIS data used in the preparation of this manuscript were obtained from the OASIS project funded by grants P50 AG05681, P01 AG03991, R01 AG021910, P50 MH071616, U24 RR021382, R01 MH56584. See http://www.oasis-brains.org/ for more details. The AIBL data used in the preparation of this manuscript were obtained from the AIBL study of ageing funded by the Common-wealth Scientific Industrial Research Organization (CSIRO; a publicly funded government research organization), Science Industry Endowment Fund, National Health and Medical Research Council of Australia (project grant 1011689), Alzheimer's Association, Alzheimer's Drug Discovery Foundation, and an anonymous foundation. See www.aibl.csiro.au for further details. The ICBM data used in the preparation of this manuscript were supported by Human Brain Project grant PO1MHO52176-11 (ICBM, P.I. Dr John Mazziotta) and Canadian Institutes of Health Research grant MOP-34996. The IXI data used in the preparation of this manuscript were supported by the U.K. Engineering and Physical Sciences Research Council (EPSRC) GR/S21533/02 -http://www.brain-development.org/. The ABIDE data used in the preparation of this manuscript were supported by ABIDE funding resources listed at http://fcon_1000.projects.nitrc.org/indi/abide/. ABIDE primary support for the work by Adriana Di Martino was provided by the NIMH (K23MH087770) and the Leon Levy Foundation. Primary support for the work by Michael P. Milham and the INDI team was provided by gifts from Joseph P. Healy and the Stavros Niarchos Foundation to the Child Mind Institute, as well as by an NIMH award to MPM (R03MH096321). http://fcon_1000.projects.nitrc.org/indi/abide/ This manuscript reflects the views of the authors and may not reflect the opinions or views of the database providers.