Anda Diemer

Johan Ondo

Carlier

Christophe Collewet

Regulation of spatially developing free shear flow in steady desired state by closed-loop control

Entrainment and mixing processes in the turbulent layer between two adjacent fluids are commonly involved in industrial applications. In case of a spatially developing plane shear layer, induced by two parallel incident streams with velocity and fluid differences, the flow is convectively unstable and a typical noise amplifier, leading to a very high sensitivity to upstream boundary conditions. Open-or closed-flow control by acting on its upstream conditions is therefore promising to enhance or reduce mixing due to turbulence.

Coherent structures of the turbulence in plane free shear layer can be view in many flow visualization experiments. It is shown in Ref. 1 and 2 that turbulence is conducted by large spanwize roller vortices surrounded by smaller quasi-streamwize rib vortices. Roller vortices come from the development of perturbation due to Kelvin-Helmholtz instability leading to vortex shedding at Strouhal frequency and near the inflection point of the velocity profile. 3, 4 Rib vortices come from three-dimensional shear instabilities in the thin vortical braid region between the rollers. 5 Together with spanwise oscillations of rollers, they speed up the transition to fully three-dimensional turbulence.

Self similar behavior of the mixing layer occuring downstream to the transition region is inescapable. Self similar behavior means geometric affinity of the turbulence statistic profiles and linear growth of the mixing layer thickness. Nevertheless, upstream conditions influence the mixing layer development in the transition and the self-similar regions, as for example the length of the transition region and the constant spreading rate in the self-similar region.

Putting aside effects of velocity and density ratios between the two streams, many studies showed periodic excitations act significantly on the mixing layer development. They trigg the roller vortex shedding and therefore modify the structure of the flow in the mixing layer and its downstream evolution, at least in the transition region. Open-loop control was achieved in Ref. 6,7,8 using a predeterminated periodic excitations. Closed-loop control of the mixing layer, more robust and efficient due to the use of measures in the loop which allow to adapt and check the control, has been recently investigated in Ref. 9, 10, 11. In Ref. 9 and 10, a temporally developing mixing layer was used to make easier the design of a predeterminated control optimizing the thickness of the mixing layer, this predeterminated control being applied in open-loop on the spatially developing mixing layer. In Ref. 11, different closed-loop approaches (extremum-seeking adaptive control, POD mode feedback control, machine learning control) were achieved by finding the optimal frequency actuation online (rather than offline as previous studies using open-loop control) that enhanced turbulence in the mixing layer.

In the present study, we achieved a classical control scheme from control theory, i.e. state feedback by Linear-Quadratic Regulator (LQR), able to maintain mixing-layer in a desired state and to reject upstream perturbations. Note that we just considered here a steady desired state as a first step in the way to cope with an unsteady desired state as a particular space-time organization of the flow knowing the benefits for industrial applications. The paper is structured as follow. Section II is devoted to fluid flow modeling and 1 of

the closed-loop control design, including an original model reduction based on simple fluid flow assumption. Section III presents the behaviour of the linearized model and the control law applied to a Navier-Stokes 3D solver. The ability of the control law to reject disturbances is also discussed.

II. Control design

II.A. Governing equations

In incompressible flow and Newtonian fluid, the motion of viscous fluid is described by the Navier-Stokes equations

∂ t u + (u • ∇)u = -1 ρ ∇p + 1 Re ∆u, ∇ • u = 0, ( 1
)
where u is the velocity field vector defined as [u v w] T and p the pressure, both being non-dimensionalized using the velocity difference between the two free stream U 0 = U 2 -U 1 and the initial thickness of the shear layer δ 0 . In this system, the dynamic viscosity of the fluid appears as the Reynolds number Re = U0δ0 ν . The flow was assumed to be near the desired profile U b which corresponded to the steady 2D incompressible base-flow solution so that

(U b • ∇)U b = -1 ρ ∇p + 1 Re ∆U b , ∇ • U b = 0. ( 2
)
In this study, the Blasius solution of this system with free shear layer boundary conditions was approximed by

U b = U 1 + 1 2
(tanh(y -y 0 ) + 1) , 0, 0

T . (3)
It is of note that this desired profile is homogeneous in the streamwise direction x, and then does not take into account the slight physical extension of the shear layer downstream. For the given profile

U b = [U b 0 0] T disturbed by u' = [u v w ]
T and p with u' U b , the Navier-Stokes equations (1) become

         ∂ t u + U b ∂ x u + v ∂ y U b = -∂ x p + 1 Re ∆u , ∂ t v + U b ∂ x v = -∂ y p + 1 Re ∆v , ∂ t w + U b ∂ x w = -∂ z p + 1 Re ∆w , ∂ x u + ∂ y v + ∂ z w = 0. (4)
The state space model obtained from equations (4) leads to a descriptor system (due to the presence of the pressure p in the equations). Even there exists a considerable study on descriptor systems, we chose to determinate a standard linear system using vorticity-stream function formulation, see Ref. 12, 13. The stream function ψ(x, y, t) is defined with 2D and incompressible flow hypothesis such that

u = ∂ y ψ (5) v = -∂ x ψ (6)
ω z = -∇ 2 ψ, (7)
where ω z is the vorticity in spanwise direction z. This formulation allows to eliminate the pressure in the equations and will allow to design a standard state space model. The vorticity-stream function formulation is given by

(∂ t + U b ∂ x )∆ -U b ∂ x - 1 Re ∆ 2 ψ = 0. ( 8
)
where U b (y) is the second derivative of U b (y) with respect to y. It is of note that plane wave perturbation leads eq. ( 8) to the well-know Rayleigh or Orr-Sommerfeld equations. 12

II.B. Boundary conditions

The extent of the considered domain was L x × L y = 64 × 16. Dirichlet boundary conditions are used at the top and bottom boundaries with u' (x, y -y 0 = ±8, t) = 0.

Input and output boundaries were defined separately. The output boundary was used to preserve the continuity at x = 64, so that

∂ x u' (x = 64, y, t) = 0. ( 10
)
and the input boundary was modeled to take in account the control. The input boundary can be modeled as a function χ (x, y, t) to take into account the overall effect of downstream actuators (jet, plasma, wind blower. . . ). Without loss of generality, this function was decomposed as a separable function such that

χ(x, y, t) = f (x)Ω(y)Γ(t), (11)
where

• f (x) represents the effects of the input control in streamwise direction and is defined such that for x = 0, χ = Ω(y)Γ(t), (12a)

for x = 0, ∂ x χ = 0, (12b)
for x = 64, χ = 0, (12c)

for x = 64, ∂ x χ = 0. ( 12d
)
These conditions lead to the function f (x) = 2( x 64 ) 3 -3( x 64 ) 2 + 1, see Ref. 13; • Ω (y) represents the control of two wind blowers used to get the two co-flowing streams U 1 and U 2 .

Comparing to the literature (see Ref.

3), Ω(y) is often given by Λ exp -σ|y-y0| .

In this study, we kept U1+U2 2 constant so that we chose Ω (y) = tanh(y -y 0 );

• Γ (t) represents the amplitude of the control.

Taking into account the input boundary, we converted the homogeneous equation ( 8) with unhomogeneous conditions (eq. in streamwise direction) into an unhomogeneous system with homogeneous conditions, see Ref. 13, 14, 15. For this purpose, the stream function was broken down into

ψ = ψ h + χ, (14)
with χ the change applied at the boundary and ψ h the corresponding inside stream function. Substituting ψ in the linearized Navier-Stokes equation ( 8), this equation becomes

     F (ψ h ) = -F (χ), ψ h (x = 0, y, t) = 0, ψ(x, y = ±1) = 0, (15)
where F (ψ) is the left side of eq. ( 8).

II.C. Model reduction

The model we propose in this paper is obtained in two steps:

• a spacial discretization (in x and y directions) of the linearized Navier-Stokes equation ( 8) using a finite differences scheme to obtain a finite model. We used n x and n y as the number of discretized sites in x and y direction, respectively;

• an extraction of the equation in y 0 (centerline of the shear layer).

Applying the finite differences scheme to (15) (simply a 2nd-order centred scheme in this study) and taking into account the input boundary, we got a reduced system

A 1 ∂ t X r + A 2 X r = B 1 ∂ t Γ(t) + B 2 Γ(t). ( 16
)
where X r is the discretized ψ h and corresponds to the state vector of dimension n xy , A 1 and A 2 are matrices of dimension n xy × n xy , B 1 and B 2 are matrices of dimension n xy × 1 and n xy = n x • n y is the total number of discretized sites.

The matrix A 1 is not singular, so that the eq. ( 16) is rewriten such that

∂ t X r = -A -1 1 A 2 X r + A -1 1 B 1 ∂ t Γ + B 2 Γ. ( 17
)
As the dynamic occuring at centerline of the shear layer lead and express the behavior of the overall shear layer, we extracted the n x rows corresponding to the dynamic in y 0 (y (1-ny)/2 = -8 and y (ny-1)/2 = 8), so we can write

∂ t X 0 = (ny-1)/2 i=(1-ny)/2 A 0i X i + B 01 ∂ t Γ + B 02 Γ ( 18
)
where X i is the discretized ψ h in y i . Matrices A 0i , B 01 and B 02 allow to capture the effects of X i , ∂ t Γ and Γ, respectively, on the dynamic of X 0 . X i need to be expressed as function of X 0 to close the system. This was done using a numerical computation of eq. ( 17) under Matlab2012 with Γ a sinusoidal function with a subharmonic. The scenario consisted to excitate the system with Γ and to measure the nodes in normal direction (the effects on streamwise direction are less important) the mean profile can be obtained for a given time. This obtained mean profile could be approximated using a Gaussian function such that

X i = e -d(yi-y0) 2 X 0 , (19)
with d 1/(2 * 0.119). Using eq. ( 18) and eq. ( 19), we can define a reduced controlled state space model given by

∂ t X(t) = AX(t) + Bq(t), (20)
with the control vector q = ∂ t Γ, the considered state vector X = X T 0 Γ T T and dynamic and control matrices are

A = (ny-1)/2 i=(1-ny)/2 A 0i e -d(yi-y0) 2 B 02 0 0 (21)
and

B = B 01 1 . ( 22
)
The properties of this linearized state is discussed in section III.

II.D. Control law

Let consider the system (20), we want to design a linear quadratic control law and therefore consider the following cost function to minimize

J = ∞ 0 (X T (t)QX(t) + q T (t)Rq(t))dt (23)
where Q, a semi-definite positive matrix (Q = Q T ≥ 0), and R, a definite positive matrix (R = R T > 0), are weighted matrices.

It is of note that the cost function was based on the the full state space, as in Ref. [START_REF] Tatsambon Fomena | Fluid Flow Control: a Vision-Based Approach[END_REF][START_REF] Gautier | Feed-forward control of a perturbed backward-facing step flow[END_REF]. Indeed, the use of image sensor in fluid mechanic enables the measure in an efficient way of the flow state, see Ref. [START_REF] Tropea | Springer Handbook of Experimental Fluid Mechanics[END_REF].

The control is computed using a feedback gain such that q = -KX (24)

with K = R -1 B T P and P is the solution of the well known algebraic Riccati equation

A T P + PA -PBR -1 B T P + Q = 0. ( 25
)
The weighted matrices Q and R are chosen to ensure a correct behavior of the closed-loop system. The matrix Q is such that X T QX is related to the turbulent kinetic energy of the system, see Ref. [START_REF] Mckernan | Control of Plane Poiseuille Flow: A theoretical and Computationnal Investigation[END_REF]. The matrix R was chosen such that R = r 2 I where r was adjusted depending on the ability of the designed control law to lead the system through the desired profile on the linearized system. In our case, I = 1 as we have a scalar control (q).

III. Results

In this section, we discuss the characterization of the linearized state space model derived from the navier-Stokes equations and the obtained results on Incompact3D with the control law eq. ( 24) obtained from the state space model.

III.A. Behavior of the state space model

The search of a control for the mixing layer assumes that the system is controllable, that means there exists a control q such that the solution X(t) → X f when t → t f . Using the Popov-Belevitch-Hautus (PBH) criteria, see, 19 it can be seen that the state space model is stabilizable.

In the particular case of the linearized system about a desired state, it's obvious that the stabilizability will be equivalent to the controllability if the final state is the desired profile (despite the numerical problem which can be encountered from sparse matrices). We note also that the state space model A is stable as expected for a noise amplifier flow. Indeed, all its eigenvalues lie in the left half complexe plane except one which is null, but their real part are very closed of the imaginary axis. The objective of our proposed control law in (24) will be to improve the stability of the closed-loop system such that we ensure the minimal kinetic energy in the system whatever disturbances are.

The stability of the state space means that in the absence of permanent disturbance, the system will converge to the desired profile U b . In the next, we look at what happens if there is permanent disturbance to its input boundary. We consider the system given by

∂ t X(t) = AX(t) + Bq(t) + Bw(t), (26)
where w(t) is the disturbance to be rejected and we assumed that disturbances and control acted on the system in the same way. We consider that w perturbs velocities U 1 and U 2 and we use the control law to compensate the effect of the disturbances. In that configuration, the disturbance acted on the input boundary and was seen as a perturbation on the control. The disturbances w(t) was formulated as a temporal sinusoidal excitation (to model fluctuations of the blowers) and was given by

w(t) = 1 k=0 a k sin(2πf k t), (27)
where f 0 is the fundamental frequency (most amplified periodic perturbation, see 3 ) and f 1 corresponds to its subharmonic (f 1 = f0 2 ), a 0 = 5e -3 and a 1 = a0 2 . We present on Fig. 1, the evolution of the input of the system (control law and disturbance) and the norm of the state vector in presence of the disturbances, with and without control. As can be seen, even if the system without control is stable, the flow is perturbated by the disturbance. The effects of this disturbance is far much lower when using our closed loop control scheme. Indeed, we remark that the obtained control was very closed to the opposite of the applied disturbance.

III.B. Behavior on the nonlinear model

To completely validated our control law we have now to apply it on the system obtained from the resolution of the 3D Navier-Stokes equations. This is done thanks to Incompact3D solver, see 20, 21 First, let us exam the behavior of the system without control and considering a perturbation w(t) described by ( 27) with a 0 = 5e -3 and f 0 = 0.08. We represent on Fig. 2 the vorticity map. As can be seen we observe the classical Kelvin-Helmholtz instabilities and pairing due to subharmonics. On the top, we represent the region of interest.The behavior we can expect from a classical excitation in terms of properties of the system (Kelvin Helmholtz instabilities, pairing due to subharmonic, . . . ) is exactly the one we obtained with our excitation. The idea here was to verify that whatever the configuration or the excitation/disturbances we encountered, we could lead the system to the desired profile.

We have to mention that we did not measure the whole system, but only data for y = 128 (corresponds to the center of the ) in the region of interest. The considered measure was done using the number of considered sites in the linearized model (n x = 17).

In this simulation and the one before, we started from the same saved configuration (simulation used to avoid difference induced by numerical initialization). The simulation proceeded such that the control was first applied alone (approximately for 3s) and the disturbance was applied when it reached 0 (near 58s in global time). We can see on Fig. 3 the effect of the proposed control law when the flow is perturbated by w(t) according to (27). The reject of the disturbance is estimated to more than 95% excepted for some periods where we just observed a rejection of 74% (the minimal rejection obtained near 184s).

Fig. 4 depicts the vorticity map. Since the control law compensates the disturbance (see Fig. 3) the observation of the first Kelvin-Helmholtz instabilities is delayed. Moreover, as can be seen, the control law avoids the pairing of vortices and even avoids the construction of vortices in the region of interest. Fig. 4 has directly to be compared with Fig. 2.

IV. Conclusion

In this paper, we showed the ability of feedback control to regulate spatially developing mixing-layer flow around a steady desired state, and to reject relatively slow disturbances compared to the acquisition rate. This was performed using a linear control law designed from a reduced linearized state space model of the Navier-Stokes equations. The reduction was based on basic assumption allowing to deal with the dynamic of the mixing layer only on the centerline. As the measure for the control loop are distributed along the centerline of the mixing layer, non-intrusive measurements method have to be used to prevent additional disturbances and wake effects on downstream measures. Optical method such as particule image velocimetry or optical flow are very appropriate, particularly if the analysis is restricted to a limited number of lines in the image. Further works will concern the case of an unsteady desired state and the validation on experiments in a dedicated wind tunnel, see Ref. 22.

Figure 1 .

 1 Figure1. Evolution of the system: on top, the control law and the disturbance; on bottom the norm of the state vector in presence of permanent disturbances.

Figure 2 .

 2 Figure 2. Instantaneous vorticity maps in the model plant in presence of permanent disturbances and without control after 202s of simulation: on top, a zoom on the region of interest; on bottom, the whole map.

Figure 3 .

 3 Figure 3. The evolution of the control in presence of permanent disturbances.

Figure 4 .

 4 Figure 4. Instataneous vorticity maps in the model plant in presence of permanent disturbances and with control after 202.5s of simulation: on top a zoom on the region of interest and on bottom the whole map.


Richard Desper
email: desper@ncbi.nlm.nih.gov

Olivier Gascuel
email: gascuel@lirmm.fr

Fast and Accurate Phylogeny Reconstruction Algorithms Based on the Minimum-Evolution Principle

 of the minimum evolution principle. For the standard version, we provide a greedy construction algorithm (GME) which produces a starting topology in O(n 2 ) time, and a tree swapping algorithm (FASTNNI) that searches for the best topology using nearest neighbor interchanges (NNIs), where the cost of doing p NNIs is O(n 2 + pn), i.e. O(n 2 ) in practice because p is always much smaller than n. The combination of GME and FASTNNI produces trees which are fairly close to Neighbor Joining (NJ) trees in terms of topological accuracy, especially with large trees. We also provide two closely related algorithms for the balanced version, called BME and BNNI, respectively. BME requires O(n 2 × diam(T )) operations to build the starting tree, and BNNI is in O(n 2 + pn × diam(T )), where diam(T ) is the topological diameter of the output tree. In the usual Yule-Harding distribution on phylogenetic trees, the diameter expectation is in log(n), so our algorithms are in practice faster that NJ. Furthermore, BNNI combined with BME (or GME) has better topological accuracy than NJ, BIONJ and WEIGHBOR (all three in O(n 3 )), especially with large trees.

Minimum evolution (ME) was proposed by several authors as a basic principle of phylogenetic inference [START_REF] Swofford | Phylogenetic inference[END_REF]. Given the matrix of pairwise evolutionary distances between the taxa being studied, this principle involves first estimating the length of any given topology and then selecting the topology with shortest length. Numerous variants exist, depending on how the edge lengths are estimated and how the tree length is calculated from these edge lengths [START_REF] Gascuel | Strengths and limitations of the minimum evolution principle[END_REF]. Several definitions of tree length have been proposed, differing from one another by the treatment of negative edge lengths. The most common solution [START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF][START_REF] Saitou | The neighbor-joining method: A new method for reconstructing phylogenetic trees[END_REF] simply defines the tree length as the sum of all edge lengths, regardless of whether they are positive or negative. Edge lengths are usually estimated within the least-squares framework. When all distance estimates are supposed to be independent and to have the same variance, we use ordinary least-squares (OLS).

The weighted least-squares corresponds to the case were distance estimates are independent but (possibly) with different variances, while the generalized leastsquares approach does not impose any restriction and is able to take benefit from the covariances of the distance estimates. Distance estimates obtained from sequences do not have the same variance, large distances are much more variable than the shortest ones [START_REF] Fitch | Construction of phylogenetic trees[END_REF], and are mutually dependent when they share a common history (or path) in the true phylogeny [START_REF] Nei | Variances of the average numbers of nucleotide substitutions within and between populations[END_REF]. Therefore, to estimate edge lengths from evolutionary distances, using generalized least-squares is theoretically superior to using weighted least-squares, which is in turn more appropriate than ordinary least-squares [START_REF] Bulmer | Use of the method of generalized least squares in reconstructing phylogenies from sequence data[END_REF].

The ME principle has been shown to be statistically consistent when combined with ordinary least-squares [START_REF] Denis | On the consistency of the minimum evolution principle of phylogenetic inference[END_REF][START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF]. However, surprisingly (and unfortunately) we recently demonstrated [START_REF] Gascuel | Strengths and limitations of the minimum evolution principle[END_REF] that the combination of this principle with weighted least-squares might be inconsistent in some (likely unrealistic) cases, while the combination with generalized least-squares is deeply inconsistent and cannot be used for phylogenetic inference.

The ME principle, combined with a form of OLS length estimation, constitutes the basis of the popular Saitou and Nei's [START_REF] Saitou | The neighbor-joining method: A new method for reconstructing phylogenetic trees[END_REF] Neighbor Joining (NJ) algorithm, which greedily refines a star tree by agglomerating taxa pairs in order to minimize the tree length at each step. The traditional approach to using ME is then to start with the NJ tree, and then do a topological search from that starting point. The first stage requires O(n 3 ), where n is the number of taxa, while the current implementations [START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF] of the second are in O(pn 3 ), where p is the number of swaps performed by the program. Another approach, available in FITCH 3.6 [START_REF] Felsenstein | PHYLIP -Phylogeny Inference Package (Version 3.2)[END_REF], is to greedily construct a tree by iterative addition of taxa to a growing tree, and then to perform swaps from this initial tree. But the time complexity is O(n 4 ) or more.

This paper further investigates the ME principle. First, we demonstrate that its combination with OLS, even when not fully optimal in terms of topological accuracy, has the great advantage to lead to very fast O(n 2 ) algorithms, able to build huge trees as envisaged in biodiversity studies. Second, we show that a new version of this principle, first introduced by Pauplin [START_REF] Pauplin | Direct calculation of a tree length using a distance matrix[END_REF] to simplify tree length computation, is more appropriate than the OLS version. In this new version, sibling subtrees have equal weight, as opposed to the standard unweighted OLS, where all taxa have the same weight so that the weight of a subtree is equal to the number of its taxa. Algorithms to deal with this new "balanced" version are faster than NJ, although not as fast as their OLS counter-part, and have a better topological accuracy than NJ, but also than BIONJ [START_REF] Gascuel | BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data[END_REF] and WEIGHBOR [START_REF] Bruno | Weighted neighbor joining: A likelihoodbased approach to distance-based phylogeny reconstruction[END_REF]. The rest of this paper is organized as follows: we first provide the notation, definitions and formulae, we describe the tree swapping algorithms for both OLS and balanced ME versions, we explain how these algorithms are modified to greedily construct a tree by iterative taxa addition, we provide simulation results to illustrate the gain in topological accuracy and run times, and we then conclude by a brief discussion. We have omitted some of the details of the algorithms and some mathematical proofs, which will appear in the journal version of the paper.

Notation, Definitions, and Formulae

A tree is made of nodes (or vertices) and of edges. Among the nodes we distinguish the internal (or ancestral) nodes and the leaves (or taxa), The leaves are denoted as i, j or k, the internal nodes as u, v or w, while an edge e is defined by a pair of nodes and a length l(e). We shall be considering various length assignments of the same underlying shape. In this case, we shall use the word "topology", while "tree" will be reserved for an instance of a topology with given edge lengths associated. We use the letter T to refer to a topology and T to refer to a tree. A tree is also made of subtrees typically denoted as A, B, C or D. For the sake of simplicity, we shall use the same notation for the subtrees and for the sets of taxa they contain. Accordingly, T also represents the set of taxa being studied, and n is the number of these taxa. Moreover, we shall use lowercase letters, e.g. a, b, c or d, to represent the subtree roots. If A and B are two disjoint subtrees, with roots a and b respectively, we'll say that A and B are distant-k subtrees if there are k edges in the path from a to b.

∆ is the matrix of pairwise evolutionary distance estimates, with ∆ ij being the distance between taxa i and j. Let A and B be two non-intersecting subtrees from a tree T . We define the average distance between A and B as:

∆ A|B = 1 |A||B| i∈A,j∈B ∆ ij .
(

∆ A|B may also be defined recursively as:

-If A and B are singleton sets, i.e. A = {a} and B = {b}, then ∆ A|B = ∆ ab , -Else, without loss of generality let B = B 1 ∪ B 2 , so we then have

∆ A|B = |B 1 | |B| ∆ A|B1 + |B 2 | |B| ∆ A|B2 . ( 2
)
It is easily seen that Equations ( 1) and ( 2) are equivalent. Equation (2) follows the notion that the weight of a subtree is proportional to the number of its taxa. Thus every taxon has the same weight, and the same holds for the distances as shown by Equation (1). Thus, this average is said to be unweighted [START_REF] Sneath | Numerical Taxonomy[END_REF]. It must be noted that the unweighted average distance between subtrees does not depend on their topologies, but only of the taxa they contain.

∆ T is the distance induced by the tree T , i.e., ∆ T ij is equal to the length of the path connecting i to j in T , for every taxon pair (i, j). Given a topology T and a distance matrix ∆, the OLS edge length estimation produces the tree T with topology T minimizing the sum of squares i,j∈T Vach [START_REF] Vach | Least squares approximation of addititve trees[END_REF], Rzhetsky and Nei [START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF] and others showed analytical formulae for the proper OLS edge length estimation, as functions of the average distances. Suppose e is an internal edge of T , with the four subtrees A, B, C and D defined as depicted in Figure 1(a). Then, the OLS length estimate of e is equal to:

(∆ T ij -∆ ij ) 2 .
l(e) = 1 2 λ(∆ A|C + ∆ B|D ) + (1 -λ)(∆ A|D + ∆ B|C ) -(∆ A|B + ∆ C|D ) , (3)
where

λ = |A||D| + |B||C| (|A| + |B|)(|C| + |D|)
.

Suppose e is an external edge, with i, A and B as represented in Figure 1(b).

Then we have:

l(e) = 1 2 (∆ A|i + ∆ B|i -∆ A|B ). (4)
Following Saitou and Nei [START_REF] Saitou | The neighbor-joining method: A new method for reconstructing phylogenetic trees[END_REF] and Rzhetsky and Nei [START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF], we define the tree length l(T ) of T to be the sum of the edge lengths of T . The OLS minimum evolution tree is then that tree with topology T minimizing l(T ), where T has the OLS edge length estimates for T , and T ranges over all possible tree topologies for the taxa being studied. Now, suppose that we are interested in the length of the tree T shown in Figure 1(a), depending on the configuration of the corner subtrees. We then have (proof deferred until journal publication):

l(T ) = 1 2 λ(∆ A|C + ∆ B|D ) + (1 -λ)(∆ A|D + ∆ B|C ) + ∆ A|B + ∆ C|D (5) +l(A) + l(B) + l(C) + l(D) -∆ a|A -∆ b|B -∆ c|C -∆ d|D ,
where λ is defined as in Equation (3). The advantage of Equation ( 5) is that the lengths l(A), l(B), l(C) and l(D) of the corner subtrees as well as the average root/leaf distances, ∆ a|A , ∆ b|B , ∆ c|C and ∆ d|D , do not depend of the configuration of A, B, C and D around e. Exchanging B and C or B and D might change the length of the five edges shown in Figure 2(a), and then the length of T , but not the lengths of A, B, C and D. This simply comes from the fact that the edge e is within the corner subtrees associated to any of the edges of A, B, C and D.

As we shall see in the next section, this property is of great help in designing fast OLS tree swapping algorithms.

Let us now turn our attention toward the balanced version of minimum evolution, as defined by Pauplin [START_REF] Pauplin | Direct calculation of a tree length using a distance matrix[END_REF]. The tree length definition is the same. Formulae for edge length estimates are identical to Equations ( 3) and ( 4), with λ replaced by 1/2, but using a different definition of the average distance between subtrees, a definition which depends on the topology under consideration. Letting T be this topology, the balanced average distance between two non-intersecting subtrees A and B is then recursively defined by:

-If A and B are singleton sets, i.e. A = {a} and B = {b} , then ∆ T A|B = ∆ ab , -Else, without loss of generality let B = B 1 ∪ B 2 . We then have

∆ T A|B = 1 2 (∆ T A|B1 + ∆ T A|B2 ). ( 6
)
The change from Equation ( 2) is that the sibling subtrees B 1 and B 2 now have equal weight, regardless of the number of taxa they contain. Thus taxa do not have the same influence depending on whether they belong to a large clade or are isolated, which can be seen as consistent in the phylogenetic inference context [START_REF] Sneath | Numerical Taxonomy[END_REF]. Moreover, a computer simulation comparing variants of the NJ algorithm [START_REF] Gascuel | On the optimization principle in phylogenetic analysis and the minimum-evolution criterion[END_REF] showed that this "balanced" (or "weighted") approach is more appropriate than the unweighted one for reconstructing phylogenies with evolutionary distances estimated from sequences. Therefore, it was tempting to test the performance of this balanced version of the minimum evolution principle.

Unfortunately, this new version does not have all of the good properties as the OLS version: the edge length estimates given by Equations ( 3) and ( 4) now depend on the topology of the corner subtrees, simply because the weighted average distances between these subtrees depend on their topologies. As we shall see, this makes the algorithms more complicated and more expensive in computing time than with OLS.

However, the same tree length Equation ( 5) holds with ∆ being replaced by ∆ T and λ by 1/2, and, fortunately, we still have the good property that tree lengths l(A), l(B), l(C) and l(D) as well as average root/leaves distances ∆ T a|A , ∆ T b|B , ∆ T c|C and ∆ T d|D remain unchanged when B and C or B and D are swapped. Edge lengths within the corner subtrees may change when performing a swap, but their (balanced) sums remain identical (proof deferred until journal publication).

FASTNNI & BNNI Tree Swapping Algorithms

This paper presents two algorithmic schemes for phylogenetic inference. The first constructs an initial tree by the stepwise addition of taxa to a growing tree, while the second improves this tree by performing local rearrangements (or swapping) of subtrees. Both follow a greedy approach and seek, at each step, to improve the tree according to the minimum evolution criterion. This approach does not guarantee that the global optimum will be reached, but only a local optimum. However, this kind of approach has proven to be effective in many optimization problems [START_REF] Cormen | Introduction to algorithms[END_REF], and we shall see that further optimizing the minimum evolution criterion would not yield significant improvement in terms of topological accuracy. Moreover, such a combination of heuristic and optimization algorithms is used in numerous phylogenetic reconstruction methods, for example those from the PHYLIP package [START_REF] Felsenstein | PHYLIP -Phylogeny Inference Package (Version 3.2)[END_REF].

In this section, we consider tree swapping to minimize the minimum evolution criterion, when edge lengths are assigned to the topology either by OLS or the balanced version. The algorithm iteratively exchanges subtrees of an initial tree, in order to minimize the length of the tree as described in Equation [START_REF] Cormen | Introduction to algorithms[END_REF]. There are many possible definitions of "subtree exchange". Since the number of combinations is high, we usually only consider exchanging neighboring subtrees, and, at least initially, we restrict ourselves to the exchange of subtrees separated by three edges, for example B and C in Figure 1(a). Such a procedure is called a "nearest neighbor interchange", since exchanging subtrees separated by one or two edges does not yield any modification of the initial tree. We adopted this approach because it allows a fast algorithm, and it is sufficient to reach a good topological accuracy.

The FASTNNI Tree Swapping Algorithm

Assume that OLS edge lengths are used, and that the swap between B and C in Figure 1(a) is considered. Let T be the initial tree and T the swapped tree. According to Equation (5) and our remarks, we have:

L(T ) -L(T ) = 1 2 (λ -1)(∆ A|C + ∆ B|D ) -(λ -1)(∆ A|B + ∆ C|D ) -(λ -λ )(∆ A|D + ∆ B|C ) , ( 7
)
where λ is as defined in Section 1, and

λ = |A||D| + |B||C| (|A| + |C|)(|B| + |D|) .
If L(T ) -L(T ) > 0, T is not optimal, and we greedily choose the tree swap which corresponds to the largest difference between L(T ) and L(T ) (there are 2 such possible swaps for each internal edge of T ). Moreover, assuming that the average distances between the corner subtrees have already been computed, L(T ) -L(T ) can be obtained in O(1) time via Equation [START_REF] Denis | On the consistency of the minimum evolution principle of phylogenetic inference[END_REF]. Instead of computing the average distances between the corner subtrees (which change when swaps are realized), we compute the average distances between every pair of nonintersecting subtrees. This takes place before evaluating the swaps and requires O(n 2 ) time, using an algorithm that is based on the recursive Equation ( 2) and is completely explained in the journal version of the paper. The whole algorithm can be summarized as follows: 1. Precompute the average distances between non-intersecting subtrees; 2. Run over all internal edges and select the best swap using Equation ( 7); 3. If the best swap is not relevant (L(T ) -L(T ) ≤ 0) stop and return the current tree, else perform the swap, compute the average distances between the newly created subtrees (A ∪ C and B ∪ D in our example above) and the other non-intersecting subtrees using Equation ( 2), and go to Step 2.

Step Rzhetsky and Nei [START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF] describe a procedure that requires O(n 2 ) to compute every edge length. In one NNI, five edge lengths are changed, so evaluating a single swap is in O(n 2 ), searching for the best swap in O(n 3 ), and their whole procedure in O(pn 3 ). This can be improved using Bryant and Waddell [START_REF] Bryant | Rapid evaluation of least-squares and minimum-evolution criteria on phylogenetic trees[END_REF] results, but the implementation in the PAUP environment of these ideas is still in progress (David Bryant, personal communication). In any case, our O(n 2 ) complexity is optimal since it is identical to the data size.

The BNNI Tree Swapping Algorithm

The BNNI algorithm is modeled after FASTNNI, with a simple change to each step. Step 1 is changed by replacing Equation ( 2) with Equation ( 6) to calculate subtree average distances. Step 2 is changed in the balanced setting by using the equation:

L(T ) -L(T ) = 1 4 (∆ T A|B + ∆ T C|D ) -(∆ T A|C + ∆ T B|D )
instead of Equation [START_REF] Denis | On the consistency of the minimum evolution principle of phylogenetic inference[END_REF]. There is a major difference within Step 3, where average distances between subtrees are updated. Consider the swap of subtrees B and C as depicted in Figure 2. To maintain an accurate table of average distances, BNNI must recalculate ∆ T X|Y . Suppose there are l -2 edges in the path from y to u in T . We then compute

∆ T X|Y = ∆ T X|Y + 2 -l ∆ T X|C -2 -l ∆ T X|B .
We'll need to do such a recomputation for every pair (x, y) in the same subtree (here A, see Figure 1(a) and Figure 2) of the tree swapping step. Since there are at most n × diam(T ) such pairs for any tree T , the time complexity of the updating step is O(n×diam(T )). The upper bound is worst when T is a chain. In this case, the diameter is proportional to n and both the number of distances to update and the bound are proportional to n 2 . However, the diameter is usually much lower. Assuming, as usual, a Yule-Harding speciation process [START_REF] Harding | The probabilities of rooted tree-shapes generated by random bifurcation[END_REF][START_REF] Yule | A mathematical theory of evolution, based on the conclusions of Dr[END_REF], the expected diameter is O(log(n)) [START_REF] Erdös | A few logs suffice to build (almost) all trees: Part II[END_REF], which implies an average complexity of the updating step in O(n log(n)). Other (e.g. uniform) distributions on phylogenetic trees are discussed in [START_REF] Aldous | Stochastic models and descriptive statistics for phylogenetic trees from Yule to today[END_REF][START_REF] Mckenzie | Distributions of cherries for two models of trees[END_REF], with expected diameter at most O( √ n). Therefore, the total time complexity of the algorithm is O(pn 2 ) in the worst case, and

O(n 2 + pn log(n)) (or O(n 2 + pn √ n)) in practice.
As we shall see, BNNI is a bit slower than FASTNNI, but is still able to improve thousand-taxa trees.

Neither FASTNNI nor BNNI explicitly computes the edge lengths until all tree swaps are completed. This is done in O(n) time by using the edge length formulae (3) and ( 4) for FASTNNI, while for BNNI, we use

l(e) = ∆ T A∪B|C∪D - 1 2 (∆ T A|B + ∆ T C|D ) and l(e) = 1 2 (∆ T i|A + ∆ T i|B -∆ T A|B ),
for the external and internal edges, respectively (see Figure 1 for the notation).

The GME and BME Addition Algorithm

The tree swapping algorithms of Section 2 are very fast at moving from an initial topology to a local minimum for their respective minimum evolution criterion. Rzhetsky and Nei [START_REF] Rzhetsky | A simple method for estimating and testing minimumevolution trees[END_REF] proposed the Neighbor-Joining tree as a good initial starting point. However, Neighbor-Joining requires O(n 3 ) computations, more than FASTNNI or BNNI would require. In this section, we will consider an insertion algorithmic scheme, which iteratively adds taxa to partial trees by using tree swapping techniques similar to those from Section 2 to find the optimal insertion point. This simple method provides an O(n 2 ) (resp. O(n 2 × diam(T )) algorithm for finding a good initial topology for the OLS (resp. balanced) version of minimum evolution.

The GME Greedy Addition Algorithm

Given an ordering on the taxa, denoted as (1, 2, 3, . . . , n) , for k = 4 to n we create a tree T k on the taxa set (1, 2, 3, . . . , k). We do this by testing each edge For each edge e ∈ T k-1 , define T e to be tree created by inserting k along e, and set c(e) to be the OLS size of T e . Computing c(e) for each e would seem to be computationally expensive. However, a much simpler approach using tree swapping exists to determine the best insertion point.

Consider the tree T of Figure 3, where k is inserted between subtrees C and A ∪ B, and assume that we have the length l(T ) = L of this new tree. If we knew the size of T , we could quickly calculate the size of T , using a special version of Equation [START_REF] Cormen | Introduction to algorithms[END_REF]:

l(T ) = L + 1 2 (λ -λ )(∆ k|A + ∆ B|C ) + (λ -1)(∆ A|B + ∆ k|C ) +(1 -λ)(∆ A|C + ∆ k|B ) . ( 8
)
In this situation,

λ = |A| + |B||C| (|A| + |B|)(|C| + 1)
,

and λ = |A| + |B||C| (|A| + |C|)(|B| + 1)
.

Suppose we were to first root the tree at some leaf r, let e 0 = (r, s) be the edge incident to r, and L to be c(e 0 ). We could then calculate c(e 1 ) and c(e 2 ) where e 1 and e 2 are the two edges incident to e 0 , and in fact we could find c(e) for each edge e ∈ T k-1 merely by repeated applications of Equation [START_REF] Erdös | A few logs suffice to build (almost) all trees: Part II[END_REF].

To simplify matters, we could write c(e) = L + f (e). Because we only seek to determine the best insertion point, we need not calculate the actual value of L, as it is sufficient to minimize f (e) with f (e 0 ) = 0.

To carry out the (2k -5) calculations of Equation ( 8), we need only precalculate 1. all average distances ∆ k|S between k and any subtree S from T k-1 ; 2. all average distances between subtrees of T k-1 separated by two edges, for example A and B in Figure 3; 3. the number of leaves of every subtree.

The algorithm can be summarized as follows:

-For k = 3, initialize the matrix of average distances between distant-2 subtrees and the array counting the number of taxa per subtree. Form T 3 with leaf set {1, 2, 3}. -For k = 4 to n,

1. Compute all ∆ k|S average distances; 2. Starting from an initial edge e 0 of T k-1 , set f (e 0 ) = 0, and recursively search each edge e to obtain f (e) from Equation (8). 3. Select the best edge by minimizing f , insert k on that edge to form T k , and update the average distance between every pair of distant-2 subtrees as well as the number of taxa per subtree. -Return T n .

To achieve Step 1, we recursively apply Equation ( 2), which requires O(k) computing time. Step 2 is also done in O(k) time, as explained above. Finally, to update the average distance between any pair A, B of distant-2 subtrees, if k is inserted in the subtree A,

∆ {k}∪A|B = 1 1 + |A| ∆ k|B + |A| 1 + |A| ∆ A|B . ( 9
)
Step 3 is also done in O(k) time, because there are O(k) pairs of distant-2 subtrees, and all the quantities in the right-hand side of Equation ( 9) have already been computed. So we build T k from T k-1 in O(k) computational time, and thus the entire computational cost of the construction of T , as we sum over k, is O(n 2 ) . This is much faster than NJ-like algorithms which require O(n 3 ) operations, and the FITCH [START_REF] Felsenstein | An alternating least-squares approach to inferring phylogenies from pairwise distances[END_REF] program which requires O(n 4 ) operations. As we shall see, this allows trees with 4,000 taxa to be constructed in few minutes, while this task requires more than an hour for NJ, and is essentially impossible for FITCH on any existing and imaginable machine. This algorithm is called GME (Greedy Minimum Evolution algorithm).

The BME Addition Algorithm

This algorithm is very similar to GME. Equation (8) simplifies into:

L(T ) = L + 1 4 (∆ T A|C + ∆ T k|B ) -(∆ T A|B + ∆ T k|C ) . ( 10
)
Step 1 is identical and provides all ∆ T k|S distances by recursively applying Equation (6). Step 2 simply uses Equation (10) instead of Equation ( 8) to calculate the relative cost of each possible insertion point. The main difference is the updating performed in Step 3. Equation [START_REF] Felsenstein | PHYLIP -Phylogeny Inference Package (Version 3.2)[END_REF] only requires the balanced average distances between distant-2 subtrees, but to iteratively update these distances, we use (and update) a data structure that contains all distances between every pair of non-intersecting subtrees (as for FASTNNI).

When k is inserted into T k-1 , we must calculate the average ∆ T A|B for any subtree A of T k which contains k, and any subtree B disjoint from A. We can enumerate all such pairs by considering their respective roots. Let a be the root of A and b the root of B. Regardless of the position of k, any node of T k could serve as the root b of B. Then, considering a fixed b, any node in the path from b to k could serve as the root a. Thus there are O(k × diam(T k )) such pairs.

Given such a pair, A, B, let us consider how we may quickly calculate ∆ T A|B from known quantities. Consider the situation as depicted in Figure 4. Suppose k is inserted by creating a new node w which pushes the subtree A 1 farther away from B. Suppose there are (l -1) edges in the path from w to a, and the subtrees branching off this path are, in order from w to a, A 2 , A 3 , . . . , A l . Then

∆ T A∪k|B = 2 -l (∆ T k|B + ∆ T A1|B ) + l i=2 2 -(l+1-i) ∆ T Ai|B .
However, we already know the value of

∆ T A|B = 2 -(l-1) ∆ T A1|B + l i=2 2 -(l+1-i) ∆ T Ai|B . Thus ∆ T A∪k|B = ∆ T A|B + 2 -l (∆ T k|B -∆ T A1|B
). The upper bound on the number of pairs is worst when T k is a chain, with k inserted at one end. In this case, the diameter is proportional to k and both the number of distances to update and the bound are proportional to k 2 . However, the diameter is usually much lower. We repeat the arguments of Section 2.2: if we assume a Yule-Harding speciation process, the expected diameter is O(log(k)) [START_REF] Erdös | A few logs suffice to build (almost) all trees: Part II[END_REF], which implies an average complexity of the updating step in O(k log(k)), while a uniform distribution on phylogenetic trees yields an expected diameter of O( √ k).

Therefore, the complexity of the whole insertion algorithm is O(n 3 ) in the worst case, and O(n 2 log(n)) (or O(n 2 √ n)) in practice. As we shall see in the next section, this is still less than NJ and allows trees with thousands of taxa to be constructed within a reasonable amount of time.

Results

Protocol

We used simulations based on random trees. Parameter values were chosen so as to cover the features of most real data sets, as revealed by the compilation of the phylogenies published in the journal Systematic Biology during the last few years. This approach induces much smaller contrasts between the tested methods than those based on model trees (e.g., [START_REF] Bruno | Weighted neighbor joining: A likelihoodbased approach to distance-based phylogeny reconstruction[END_REF][START_REF] Gascuel | BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data[END_REF]). Indeed, model trees are generally used to emphasize a given property of the studied methods, for example their performance when the molecular clock is strongly violated. Thus, model trees are often extreme and their use tends to produce strong and possibly misleading differences between the tested methods. On the other hand, random trees allow comparisons with a large variety of tree shapes and evolutionary rates, and provide a synthetic and more realistic view of the average performances.

We used 24-and 96-taxon trees, and 2000 trees per size. For each of these trees, a true phylogeny, denoted as T , was first generated using the stochastic speciation process described by Kuhner and Felsenstein [START_REF] Kuhner | A simulation comparison of phylogeny algorithms under equal and unequal rates[END_REF], which corresponds to the usual Yule-Harding distribution on trees [START_REF] Harding | The probabilities of rooted tree-shapes generated by random bifurcation[END_REF][START_REF] Yule | A mathematical theory of evolution, based on the conclusions of Dr[END_REF]. Using this generating process makes T ultrametric (or molecular clock-like). This hypothesis does not hold in most biological data sets, so we created a deviation from the molecular clock. Each edge length of T was multiplied by 1.0 + µX, where X followed the standard exponential distribution (P (X > η) = e -η ) and µ was a tuning factor to adjust the deviation from molecular clock; µ was set to 0.8 with 24 taxa and to 0.6 with 96 taxa. The average ratio between the mutation rate in the fastest evolving lineage and the rate in the slowest evolving lineage was then equal to about 2.0 with both tree sizes. With 24 taxa, the smallest value (among 2000) of this ratio was equal to about 1.2 and the largest to 5.0 (1.0 corresponds to the strict molecular clock), while the standard deviation was approximately 0.5. With 96 taxa, the extreme values became 1.3 and 3.5, while the standard deviation was 0.33. These (2 × 2000) trees were then rescaled to obtain "slow", "moderate" and "fast" evolutionary rates. With 24 taxa, the branch length expectation was set to 0.03, 0.06 and 0.15 mutations per site, for the slow, moderate and fast conditions, respectively. With 96 taxa, we had 0.02, 0.04 and 0.10 mutations per site, respectively. For both tree sizes, the average maximum pairwise divergence was of about 0.2, 0.4 and 1.0 substitutions per site, with a standard deviation of about 0.07, 0.14 and 0.35 for 24 taxa, and of about 0.04, 0.08 and 0.20 for 96 taxa. These values are in good accordance with real data sets. The maximum pairwise divergence is rarely above 1.0 due to the fact that multiple alignment from highly divergent sequences is simply impossible. Moreover, with such distance value, any correction formula, e.g. [START_REF] Kimura | A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences[END_REF], becomes very doubtful, due to our ignorance of the real substitution process and to the fact that the larger the distance the higher the gap between estimates obtained from different formulae. The medium condition (∼0.4) corresponds to the most favorable practical setting, while in the slow condition (from ∼0.1 to ∼0.3) the phylogenetic signal is only slightly perturbed by multiple substitutions, but it can be too low, with some short branches being not supported by any substitution.

SeqGen [START_REF] Rambaut | Seq-Gen: An application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees[END_REF] was used to generate the sequences. For each tree T (among 3 × 2 × 2000), these sequences were obtained by simulating an evolving process along T according to the Kimura [START_REF] Kimura | A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences[END_REF] two parameter model with a transition/transversion ratio of 2.0. The sequence length was set to 500 sites. Finally, DNADIST from the PHYLIP package [START_REF] Felsenstein | PHYLIP -Phylogeny Inference Package (Version 3.2)[END_REF] was used compute the pairwise distance matrices, assuming the Kimura model with known transition/transversion ratio. The data files are available on our web page (see title page).

Every inferred tree, denoted as T , was compared to the true phylogeny T (i.e., that used to generate the sequences and then the distance matrix) with a topological distance equivalent to Robinson and Foulds' [START_REF] Robinson | Comparison of phylogenetic trees[END_REF]. This distance is defined by the proportion of internal edges (or bipartitions) which are found in one tree and not in the other one. This distance varies between 0.0 (both topologies are identical) and 1.0 (they do not share any internal edge). The results were averaged over the 2000 test sets for each tree size and evolutionary rate. Finally, to compare the various methods to NJ, we measured the relative error reductions (P M -P NJ )/P NJ , where M is any tested method different from NJ and P X is the average topological distance between T and T when using method X. This measure highlights the performance difference between the related, distance-based, tested methods.

Phylogeny Estimation Algorithm Comparison

We used a variety of different algorithms to try to reconstruct the original tree, given the matrix of estimated distances. We used the NJ [START_REF] Saitou | The neighbor-joining method: A new method for reconstructing phylogenetic trees[END_REF] program from the PAUP package [START_REF] Swofford | PAUP-Phylogenetic Analysis Using Parsimony (and other methods)[END_REF], with and without the BIONJ [START_REF] Gascuel | BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data[END_REF] flag; WEIGHBOR version 1.2 [START_REF] Bruno | Weighted neighbor joining: A likelihoodbased approach to distance-based phylogeny reconstruction[END_REF], available at http://www.t10.lanl.gov/billb/weighbor/, the FITCH [START_REF] Felsenstein | An alternating least-squares approach to inferring phylogenies from pairwise distances[END_REF] program from the PHYLIP package [START_REF] Felsenstein | PHYLIP -Phylogeny Inference Package (Version 3.2)[END_REF]; the Harmonic Greedy Triplet (HGT/FP) program, provided by Miklos Csűrös [START_REF] Csűrós | Fast recovery of evolutionary trees with thousands of nodes[END_REF]; GME and BME. Also, we used output topologies from other programs as input for FASTNNI and BNNI. GME, BME, FASTNNI, and BNNI will be available in the near future at our web page, or via ftp at ftp://ncbi.nlm.nih.gov/∼pub/desper/ME.

We also measured how far from the true phylogeny one gets with NNIs. This served as a measure of the limitation of each of the minimum evolution frameworks, as well as a performance index for evaluating our algorithms. Ordinarily, the (OLS or balanced) minimum evolution criterion will not, in fact, observe a minimum value at the true phylogeny. So, starting from the true phylogeny and running FASTNNI or BNNI, we end up with a tree with a significant proportion of false edges. When this proportion is high, the corresponding criterion can be seen as poor regarding topological accuracy. This proportion represents the best possible topological accuracy that can be achieved by optimizing the considered criterion, since we would not expect any algorithm optimizing this criterion in the whole tree space to find a tree closer from the true phylogeny than the tree that is obtained by "optimizing" the true phylogeny itself.

Results are displayed in Table 1 and Table 2. The first number indicates the average topological distance between the inferred tree and the true phylogeny.

The second number provides the percentage of relative difference in topological distance between the method considered and NJ (see above); the more negative this percentage, the better the method was relative to NJ.

The performances of the basic algorithms (1st columns) are strongly correlated with the number of computations that they perform. Both O(n 2 ) algorithms are clearly worse than NJ. BME (in O(n 2 log(n))) is still worse than NJ, but becomes very close with 96 taxa, which indicates the strength of the balanced minimum evolution framework. BIONJ, in O(n 3 ) like NJ and having identical computing times, is slightly better than NJ in all conditions, while WEIGH-BOR, also in O(n 3 ) but requiring complex numerical calculations, is better than BIONJ. Finally, FITCH, which is in O(n 4 ), is the best with 24 taxa, but was simply impossible to evaluate with 96 taxa (see the computing times below).

After FASTNNI (2nd columns), we observe that the output topology does not depend much on the input topology. Even the deeply flawed HGT becomes close to NJ, which indicates the strength of NNIs. However, except for the true phylogeny in some cases, this output is worse than NJ. This confirms our previous results [START_REF] Gascuel | On the optimization principle in phylogenetic analysis and the minimum-evolution criterion[END_REF], which indicated that the OLS version of minimum evolution is reasonable but not excellent for phylogenetic inference. But this phenomenon is much more visible with 24 than with 96 taxa, so we expect the very fast GME+FASTNNI O(n 2 ) combination to be equivalent to NJ with large n.

After BNNI (3rd columns), all initial topologies become better than NJ. With 24 taxa, the performance is equivalent to that of FITCH, while with 96 taxa the results are far better than those of WEIGHBOR, especially in the Fast condition where BNNI improves NJ by 21%, against 10% for WEIGHBOR. These results are somewhat unexpected, since BNNI has a low O(n 2 + pn log(n)) average time complexity. Moreover, as explained above, we do not expect very high contrast between any inference method and NJ, due to the fact that our data sets represent a large variety of realistic trees and conditions, but do not contain extreme trees, notably concerning the maximum pairwise divergence. First experiments indicate that maximum parsimony is close to BNNI in the Slow and Moderate conditions, but worse in the Fast condition, the contrast between these methods being in the same range as those reported in Table 1 and2. The combination of BNNI with BME or GME can thus be seen as remarkably efficient and accurate. Moreover, regarding the true phylogeny after BNNI, it appears that little gain could be expected by further optimizing this criterion, since our simple optimization approach is already close to these upper values.

Computing Times

To compare the actual running speeds of these algorithms, we tested them on a Sun Enterprise E4500/E5500, with ten 400 MHz processors and 7 GB of memory, running the Solaris 8 operating system. Table 3 summarizes the average computational times (in hours:minutes:seconds) required by the various programs to build phylogenetic trees. The leftmost two columns were averaged over two thousand 24-and 96-taxon trees, the third over ten 1000-taxon trees, and the final column over four 4000-taxon trees. Stars indicate entries where the algorithm was deemed to be feasible FITCH would take approximately 25 minutes to make one 96-taxon tree, a speed which makes it impractical for real applications, because they often include bootstrapping which requires the construction of a large number of trees. As WEIGHBOR's running time increased more than 60-fold when moving from 24-taxon trees to 96-taxon trees, we judged it infeasible to run WEIGHBOR on even one 1000-taxon tree.

Table 3 confirms that all of our algorithms are faster than NJ, with the improvement becoming notably impressive when thousands of taxa are dealt with. Furthermore, we suspect that implementation refinements such as those used in PAUP's NJ could be used to make our algorithms still much faster.

Table 4 contains the number of NNIs performed by each of the three combinations which appear in Table 3. Not surprisingly, the largest number of NNIs was consistently required when the initial topology was made to minimize the OLS ME criterion, but NNIs were chosen to minimize the balanced ME criterion (i.e., GME + BNNI). This table shows the overall superiority of the BME tree over the GME tree, when combined with BNNI. In all of the cases considered, the average number of NNIs considered for each value of n was considerably less than n itself.

Discussion

We have presented a new distance-based approach to phylogeny estimation which is faster than most distance algorithms currently in use. The current most popu-lar fast algorithm is Neighbor-Joining, an O(n 3 ) algorithm. Our greedy ordinary least-squares minimum evolution tree construction algorithm (GME) runs at O(n 2 ), the size of the input matrix. Although the GME tree is not quite as accurate as the NJ tree, it is a good starting point for nearest neighbor interchanges (NNIs). The combination of GME and FASTNNI, which achieves NNIs according to the ordinary least-squares criterion, also in O(n 2 ) time, has a topological accuracy close to that of NJ, especially with large numbers of taxa and moderate or slow rates.

However, the balanced minimum evolution framework appears much more appropriate for phylogenetic inference than the ordinary least-squares version. This is likely due to the fact that it gives less weight to the topologically long distances, i.e. those containing numerous edges, while the ordinary least-squares method puts the same confidence on each distance, regardless of its length. Even when the usual and topological lengths are different, they are strongly correlated. The balanced minimum evolution framework is thus conceptually closer to weighted least-squares [START_REF] Fitch | Construction of phylogenetic trees[END_REF], which is more appropriate than ordinary least-squares for evolutionary distances estimated from sequences.

The balanced NNI algorithm (BNNI) achieves a good level of performance, superior to that of NJ, BIONJ and WEIGHBOR. BNNI is an O(n 2 + np × diam(T )) algorithm, where diam(T ) is the diameter of the inferred tree, and p the number of swaps performed. With p × diam(T ) = O(n) for most data sets, the combination of GME and BNNI effectively gives us an O(n 2 ) algorithm with high topological accuracy.

Fig. 1 .

 1 Fig. 1. Corner subtrees used to estimate the length of e, (a) for e an internal edge; (b) for e an external edge

Fig. 2 .

 2 Fig. 2. T is obtained from T by swapping B and C

Fig. 3 .

 3 Fig. 3. T is obtained from T by swapping A and k

Fig. 4 .

 4 Fig. 4. Calculating balanced average ∆ TA|B when k is inserted into A

Table 1 .

 1 Topological accuracy for 24-taxon trees at various rates of evolution

		slow rate		moderate rate
		w/o NNIs FASTNNI BNNI w/o NNIs FASTNNI BNNI
	True Tree	.109 -1.6% .104 -6.2%	.092 3.7% .083 -5.8%
	FITCH .109 -1.9% .113 2.0% .107 -3.4% .085 -4.9% .094 6.0% .085 -4.0%
	WEIGHBOR .109 -1.8% .112 1.7% .107 -3.0% .085 -4.3% .094 6.2% .085 -4.0%
	BIONJ .111 -0.3% .113 2.0% .107 -3.6% .087 -2.0% .094 6.5% .085 -4.2%
	NJ .111 0% .113 2.0% .107 -3.5% .088 0% .094 6.6% .085 -4.0%
	BME .118 7.1% .113 1.9% .107 -2.8% .100 13% .094 6.3% .084 -4.9%
	GME .122 10% .113 2.1% .107 -3.4% .107 21% .095 7.1% .084 -4.8%
	HGT/FP .334 202% .112 1.1% .107 -2.9% .326 268% .095 7.5% .088 -0.2%
			fast rate
		w/o NNIs FASTNNI BNNI
		True Tree	.088 6.5% .076 -8.3%
		FITCH .076 -7.8% .090 8.8% .077 -7.0%
		WEIGHBOR .077 -6.8% .089 8.2% .077 -6.9%
		BIONJ .079 -3.6% .090 9.0% .077 -6.6%
		NJ .082 0% .090 9.1% .077 -6.9%
		BME .098 19% .090 9.1% .076 -7.1%
		GME .105 28% .090 9.8% .076 -7.6%
		HGT/FP .329 300% .090 9.8% .083 0.8%

Table 2 .

 2 Topological accuracy for 96-taxon trees at various rates of evolution

		slow rate		moderate rate
		w/o NNIs FASTNNI BNNI w/o NNIs FASTNNI	BNNI
	True Tree	.172 -5.6% .167 -8.8%	.132 -3.0% .115 -15.4%
	WEIGHBOR .178 -2.5% .181 -0.7% .173 -5.2% .129 -5.4% .137 0.5% .118 -13.0%
	BIONJ .180 -0.9% .182 -0.3% .173 -5.1% .134 -1.9% .138 1.3% .118 -13.0%
	NJ .183 0% .182 -0.2% .173 -5.2% .136 0% .139 1.8% .119 -12.9%
	BME .186 1.9% .181 -0.6% .173 -5.3% .137 1.0% .138 1.1% .118 -13.2%
	GME .199 8.8% .183 0.3% .173 -5.3% .158 16% .140 2.7% .118 -13.2%
	HGT/FP .512 185% .185 1.5% .175 -4.3% .480 253% .143 5.2% .123 -.9.3%
			fast rate
		w/o NNIs FASTNNI	BNNI
		True Tree	.115 0.6% .088 -23.4%
		WEIGHBOR .103 -10% .119 3.8% .091 -21.0%
		BIONJ .112 -2.5% .121 5.1% .090 -21.7%
		NJ .115 0% .121 5.5% .090 -21.3%
		BME .117 1.8% .120 4.4% .090 -21.4%
		GME .144 25% .122 6.3% .091 -21.1%
		HGT/FP .465 306% .126 9.4% .098 -14.7%

Table 3 .

 3 Average computational time (HH:MM:SS) for reconstruction algorithms

	Algorithm	24 taxa 96 taxa 1000 taxa 4000 taxa
	GME	.02522 .07088	8.366	3:11.01
	GME + FASTNNI .02590 .08268	9.855	3:58.79
	GME + BNNI	.02625 .08416	11.339	6:02.10
	HGT/FP	.02518 .13491	13.808	3:33.14
	BME	.02534 .08475	19.231 12:14.99
	BME + BNNI	.02557 .08809	19.784 12:47.45
	NJ	.06304 .16278	21.250 20:55.89
	BIONJ	.06528 .16287	21.440 20:55.64
	WEIGHBOR	.42439 26.88176	*****	*****
	FITCH	4.37446	*****	*****	*****

Table 4 .

 4 Average number of NNIs performed

	Algorithm	24 taxa 96 taxa 1000 taxa 4000 taxa
	GME + FASTNNI 1.244 8.446	44.9	336.50
	GME + BNNI	1.446 11.177	59.1	343.75
	BME + BNNI	1.070 6.933	29.1	116.25

Acknowledgements

Special thanks go to Stephane Guindon, who generated the simulated data sets used in Section 4.


Radajewski Dimitri

Pierre Roblin

Patrice Bacchin

Martine Meireles

Yannick Hallez

In situ structural analysis with a SAXS laboratory beamline directly in a microfluidic chip

Recent advances have been made in coupling microfluidic chips with X-ray equipment enabling structural analysis of samples directly in microfluidic devices.This important step mainly took place on powerful synchrotron facilities because of the need of a beam reduced in size to fit the microfluidic channel dimensions but still intense. In this work, we discuss how improvements of an X-ray laboratory beamline and an optimal design of a microfluidic device allow to obtain reliable structural information without the need of a synchrotron. We evaluate the potential of these new developments by probing several well known dispersions. These include dense inorganic gold and silica nanoparticles that scatter photons quite intensely, the bovine serum albumin (BSA) macromolecule, with more moderate contrast, to highlight possible applications in biology, and latex nanospheres with only weak contrast with the solvent to show the limits of the setup. We established a proof of concept for a versatile setup that will open the way for more complex lab-on-chip devices suitable for in situ and operando structural analysis by small angle X-ray scattering analysis without the necessity of a synchrotron source.

Introduction

Small-angle X-ray scattering (SAXS) has now been proven to be a powerful technique to investigate the structures and interactions of hard and soft matter systems. [START_REF] Feigin | Structure analysis by small-angle X-ray and neutron scattering[END_REF][START_REF] Glatter | Small angle X-ray scattering[END_REF][START_REF] Guinier | Small-angle scattering of X-rays[END_REF] It can provide a wide range of information and has been used extensively over the years to determine molecular interactions [START_REF] Ducruix | [END_REF] and molecular weights, 5 , investigate structural conformations, 6,7 formation of nanoparticles [8][9][10] or folding dynamics of proteins in solution. 11,12 Lab-on-chip devices have been developed more recently and offer the possibility to perform various operations at the microscale, using only small quantities of sample. 13,14 They are found in diverse scientific fields, such as biology, biochemistry, chemistry or physics. Among many examples, they can be used to mix fluids, 15,[START_REF] Ward | structures, devices, and systems[END_REF] perform filtration [17][18][19] or compression 20 operations, observe crystallization phenomena, [21][22][23] separate blood cells 24,25 or carry out high-throughput screening. [26][27][28] Coupling these two techniques has been recently successfully used to study nanoparticle formation, [29][30][31][32][33][34] protein folding and unfolding, [35][36][37][38] self-assembly of biomacromolecules 39,40 and lipid nanoparticles 41 or protein crystallization. 42,43 Microfluidic systems allow to perform unit operations with great control and reproducibility, with the possibility for segmented-flow configurations to quickly screen numerous experimental conditions and access to structural information in the compartmentalized Laboratoire de Génie Chimique, Université de Toulouse, CNRS, INPT, UPS, Toulouse, France. E-mail: dimitri.radajewski@gmail.com, yannick.hallez@univ-tlse3.fr units. Flowing samples ensures a constant renewal of the exposed sample and therefore avoids radiation damage. 44 Measuring at different positions in compartmentalized units allows in situ time-resolved analysis since the position in the channel and the timescale investigated are related by the flowrate. However, since the X-ray beam size has to fit the small size of the channel and at the same time suitable signal-to-noise ratio is needed, such experiments are commonly carried out at facilities equipped with high-brilliance synchrotron sources. Since access to them is limited, conducting long-term studies or performing experimental campaigns on demand is difficult.

At the same time, major technical developments concerning the performances of X-ray sources, detectors, X-ray optics and scatterless slits for laboratory X-ray beamlines 45,46 have opened the field for the use of microfocused beams at the laboratory scale. Microfluidic systems thus began to be used in conjunction with laboratory X-ray equipment to investigate the nucleation and growth of zinc oxide 47 and gold nanoparticles. 48,49 However, in these experiments, the microfluidic chip is basically used to synthesize the nanoparticles while the X-ray scattering is actually performed on a millimetric capillary glued at the exit of the chip. To our knowledge, there is a single report of a microfluidic device (made of PMMA, PTFE and Polyimide) directly mounted in a chamber of a laboratory X-ray instrument. 50 In this work, concentrated and dilute inorganic slurries (10 wt% silica nanoparticles vs. 0.5 wt% CaCO3 crystals) have been probed in segmented water-in-oil flows. The authors reported short time exposure SAXS images for concentrated silica nanoparticles contained in aqueous droplets flowing along a 300×300 µm2 channel. The quality of produced scattering pattern allowed a good agreement with a model for spherical particles. However for dilute calcite crystals almost zero scattering intensity was observed with small exposure times.

The advantages of using microfluidics are linked to the small lateral dimensions and depth of the channels. This imposes two major limitations when performing SAXS on a microfluidic chip: a small beam size to fit the lateral microfluidic channels dimensions, and a limited number of potential scatterers due to the small channel depth. This is not a problem when working with high brilliance synchrotron sources because the high flux will generate lots of scattered photons anyway, but it is a major issue limiting the joint use of laboratory beamlines with microfluidics. On laboratory beamlines, it is therefore very important to optimize X-ray optics (photon flux, beam size, smearing at small angles) together with microfluidic chip characteristics (channel geometry, materials, wall thickness) to minimize absorption and background scattering from the chip.

In this work, we report how a microfluidic device was developed and how a state-of-the-art laboratory beamline was optimized to perform in situ, operando SAXS measurements without the need for a synchrotron source. We evaluate the potential of these new developments by probing several well known types colloidal dispersions. These include dense inorganic gold and silica nanoparticles that scatter photons quite intensely, the bovine serum albumin (BSA) macromolecule, with more moderate contrast, to highlight possible applications in biology, and latex nanospheres with only weak contrast with the solvent to show the limits of the setup. This discussion is complementary to the one proposed in Ref. 50, in which additional constraints on the acquisition times allowed to have exploitable signals only for inorganic nanoparticles.

Material and methods

Microfluidic system

The microfluidic chip used in this work was fabricated using soft lithography and injection molding techniques.

First, the 2D design of the chip was produced with the Computer Aided Design software AutoCAD. It consists of four parallel, 400 µm wide, 3 cm long channels connected to inlet and outlet tubing (see Figure 1). The channel width was chosen to accommodate a relatively large beam size of about 250 µm, which allows to maximize the photon flux on the sample. Black circles are on the mask to form pillars that will favor a homogeneous distribution of OSTEMER during injection molding.

The next step is the fabrication of a master mold. Several layers of inexpensive multilevel negative tone photoresist dry film (WBR2000 series, DuPont, France) are first laminated on a glass substrate to obtain the required 400 µm thickness (Thermo Scientific Menzel-Glaser, Germany). The 2D chip design is then patterned on the dry film by UV exposure (UV-KUB2, Kloé, Montpellier, France) through a low cost emulsion mask, and 3D structures are obtained by developing the non-exposed photoresist using a 1%w solution of sodium carbonate (Na 2 CO 3 ). After rinsing by tap water, the master is dried and hydrophobized with NOVEC TM 1720 (3M, France).

A hollow PDMS (Sylgard 184 elastomer kit, Dow Corning, USA) mold is then produced by pouring liquid PDMS on the surface of the dry film master, curing it for few hours at 70°C before gently and easily peeling it off from the hydrophobic master mold. The procedure to elaborate the dry film master mold and the PDMS mold is described in more details in Ref. 42 for the interested reader.

To assemble the final microfluidic chip, two OSTEMER (Mercene Labs, Sweden 51 ) stickers have to be prepared. First, a transparent mylar film is placed on top of the hollow PDMS mold and liquid OSTEMER is injected with a needle in the gap between the mold and the film (see Figure 1). OSTEMER is then UV-cured, becoming solid, but still flexible and sticky, so the structure can be peeled off from the PDMS mold. Second, another OSTEMER sticker is obtained by injecting and curing it in a flat 50 µm deep PDMS mold on top of which a 25 µm polyimide Kapton film is positioned. Finally, the two OSTEMER stickers are assembled together and a second curing of the whole chip is made by placing it in an oven at 90°C overnight. The mylar film used as a fabrication support can finally be removed to get the final chip made of only OSTEMER supported by a Kapton film.

PEEK capillaries 360 µm OD 150 µm ID (Trajan Scientific, United Kingdom) are used as connectors and epoxy-glued to the inlet and outlets of the chip. They are connected to PTFE tubing 1/16" OD with specific adaptors and colloidal dispersions are manually injected in the device with a syringe.

As this work is concerned by the optimization of a microfluidic chip for a use on a laboratory X-ray beamline with low photon flux compared to a synchrotron, the important point is that X-rays will pass through 25 µm of Kapton, two 50 µm thick layers of OSTE-MER, and 400 µm of sample. Maximizing the sample thickness compared to the wall thickness is of course essential. Increasing the sample thickness further would be possible but the device would become closer to a standard capillary than to a microfluidic chip. Polyimide (Kapton) films and OSTEMER are known for their very good compatibility with SAXS measurements in particular due to their low attenuation and low background scattering. 31 The Kapton layer is present to provide some mechanical rigidity to the chip, and decreasing its thickness is not an option. Cutting a small window in this layer (or another support material) so the X-rays pass only through OSTEMER and the sample is however possible. Trying to reduce the OSTEMER thickness below 50 µm led to fluid leakage through microcracks, so the present thickness seems to be optimal to maximize the photon flux and minimize background scattering. Small-Angle X-ray scattering SAXS experiments were performed on a Xeuss 2.0 laboratory beamline from Xenocs equipped with a microfocus GeniX 3D Cu microsource and a Pilatus3 1M HPC detector from DECTRIS. Before modification, the X-ray beam was collimated with FOX3D 12_INF single reflection multilayer optics, and cut vertically and horizontally by two scatterless slits to minimize beam divergence. Two modifications of this setup have been performed to increase the photon flux while minimizing the beam size on the sample. First, the optics were replaced by a FOX3D HFVL mirror which allowed to increase the flux from 80 to 450 million photons per A similar issue was mentioned by Levenstein and coworkers 50 using the same Xeuss 2.0, with the same beam size, although their more powerful liquid-metal-jet source allowed them to maintain a flux of 3.7 × 10 6 ph s -1 . By inserting a pinhole located very close to the sample in our setup, the slits could be opened more than on the original setup and the flux on the sample could be increased back to about 9.5 × 10 6 ph s -1 . The flux on the sample is thus maximized and possible artifacts due beam divergence are minimized. With this setup, the divergence of the beam after the pinhole is about 1 mrad. The microfluidic chip is mounted on a support, fixed on a precise motorized translation stage such that it can be readily moved in horizontal and vertical positions with respect to the X-ray beam in order to position the investigated channel in front of the beam.

The sample to detector distance was 1.2195 m and the wavelength is λ = 1.541 Å giving a q-range extending from 0.005 Å -1 to 0.5 Å -1 . All measurements were performed at atmospheric pres-sure, 21 °C, and ambient humidity. For the relatively low scattering materials (proteins and latex), 6 acquisitions of 10 minutes were performed and the raw 2D scattered image was obtained by averaging the 6 recorded frames. For the most scattering materials (gold and silica nanoparticles), single acquisitions of 10 minutes were sufficient to have a decent scattering pattern. The recorded images were corrected for detector efficiency and distortion, scaled by the acquisition time, sample thickness, and sample transmission factor, and averaged in the azimuthal direction while applying a mask to remove faulty regions of the image to finally obtain the intensity profiles I s (q).

The intensity q of the scattering vector is

q = 4π λ sin θ 2 ( 1
)
where θ is the scattering angle.

For each experiment, a background scattering profile I b (q) was also obtained from the buffer, using the same standardized scaling procedure. The final intensity profiles reported hereafter are defined as I(q) = I s (q) -I b (q). For single scattering they can be written as

I(q) = Aφ P(q)S m (q) ( 2
)
where A is a q-independent constant, φ is the colloid volume fraction, P(q) is the form factor and S m (q) is the so-called measurable structure factor. At large dilution S m (q) → 1 so a measurement in these conditions gives I 0 (q) ≡ Aφ 0 P(q). A measurement at finite concentration then allows to define S m (q) and A from (2) and the requirement S m (q → ∞) = 1. Note that S m (q) is the true structure factor only for monodisperse systems without correlations between positions and orientations. For polydisperse or anisotropic systems, S m (q) is merely an easily measurable quantity with no clear thermodynamic meaning.

Preparation of solutions

Gold nanoparticles were synthesized by collaborators (see Acknowledgements) using an adapted protocol from Ref. 52. A mixture of oleylamine and HAuCl4xH2O in pentane was stirred at room temperature for 5 min. Afterward, a solution of tertbutylamine borane in presence of oleylamine in pentane was added to reduce the gold precursor. The color of the solution immediately turned to brown. After stirring for 60 min at room temperature, the nanoparticles were purified by precipitating with ethanol and centrifuged. The supernatant was discarded and the precipitated NPs were redispersed in pentane at a volume fraction of 1.3 × 10 -4 . The aqueous dispersion of Ludox HS-40 (Sigma-Aldrich, France) colloidal silica 40 wt. % was filtered through a 5 µm teflon filter and then dialyzed for 15 days in a 7 kDa cellulose membrane bag (Sigma-Aldrich, France) against a 1mM potassium chloride solution. Subsequent dilutions of the washed dispersion were obtained with the 1mM KCl solution to reach volume fractions φ = 0.00234, 0.00496, 0.00826, 0.0144, 0.0287, 0.0420 and 0.0907. The volume fraction was determined by weighting a 2 g sample before and after drying at 120°C. For this calculation we assumed a silica density of 2200 kg m -3 .

Protein Bovine Serum Albumin (ThermoFisher Scientific, France) was used as a model macromolecule and was solubilized at a stock concentration of 100 mg/mL in phosphate-buffered saline PBS (Fisher BioReagents, USA). Subsequent dilutions led The aqueous dispersions of polystyrene latex beads with sulfate charges (Invitrogen, ThermoFisher Scientific, USA) were prepared using a 1 mM NaCl solution to reach volume fractions φ = 0.034, 0.07, 0.105, and 0.14. The average radius announced by the manufacturer was 13.5 nm.

Results and discussion

Gold nanoparticles

The first samples used to test the setup were dispersions of gold nanoparticles in pentane. These colloids are dense, spherical, and quite monodisperse so they are expected to have a very clear SAXS signature even at very low volume fraction.

Results are presented in Figure 3 together with model spectra corresponding to sphere suspensions in which particle radii follow a Schulz distribution 53 with average 2.75 nm and 17 % polydispersity. This is in good agreement with the specification from the synthesis protocol. 52 . The absolute intensity of the model (coefficient A in (2)) has been adjusted only once on the data corresponding to the highest concentration. The model spectra for the three other samples then follow eq. ( 2) without need for further adjustments.

As expected for these dense colloids, the signal to noise ratio is quite good, and sufficient to fit a form factor with confidence. It remains however too limited to use Porod's law 2,54 based on the high-q limit of Iq 4 quantitatively. For this, and for the present dispersion, the background noise would need to be negligible up to q ≃ 0.4 Å -1 , which would require an increase of the scattered signal intensity by about two orders of magnitude. The Q invariant (integral of Iq 2 ) 2,54 could be calculated with reasonable accuracy in particular because the minimum q value has been made low enough with the use of the pinhole.

Silica nanoparticles and interactions

The test with gold nanospheres allowed to discuss pure form factor measurements. The next step is to test the possibility to mea-sure the mesostructure of a dispersion, and therefore indirectly colloidal interactions, directly inside a microfluidic chip. For this purpose we used Ludox HS40 aqueous silica dispersions. These particles scatter X-rays quite intensely due to their large contrast, which makes SAXS measurements relatively easy. They have already been investigated using synchrotron sources (see e.g. Refs. 55-57) so their SAXS signature is quite well documented, and their surface chemistry is also well known. We are concerned here with the following questions: can the present setup yield spectra with good enough signal-to-noise ratios to fit the first structure factor peak? The second one? Can we access small enough angles to measure the plateau S(q → 0), which is important for thermodynamics?

The most dilute sample (φ =0.00234) was used to measure the form factor, which is reported in the inset of Figure 4. It could be modeled using a Schulz distribution for particle radii with an average of 8 nm and 13 % size polydispersity, in line with frequently reported values. [56][57][58] The measurable structure factors of the more concentrated samples are reported in Figure 4. The apparition of a marked structure peak and the decrease of S(q → 0) with increasing particle concentration could be observed clearly enough to compare them to theoretical model and simulation predictions. These features, due to long-range electrostatic interactions, could be computed theoretically by solving the Ornstein-Zernike (OZ) equation with the Rogers-Young closure and a Hard-Sphere-Yukawa effective potential (see Ref. 59 and references therein for more details). Polydispersity was accounted for using the decoupling approximation, shown to be quite accurate for the present moderate polydispersity. 53,60 The effective colloidal charge and screening length have been calculated using the Extrapolated Point Charge renormalization method, 61 accounting for a volume-fraction-dependent charge regulation mechanism involving the surface reaction

Si-OH --⇀ ↽--Si-O -+ H + (3)
with a dissociation constant pK = 7.5, a surface concentration of ionizable sites Γ 0 = 8 nm -2 , and a Stern capacitance C S = 2.9 F.m -2 . These values have been measured independently on large silica beads by Trefalt and coworkers. 62 The pH value of 6.9 used in the regulation model was the one directly measured in the dispersions. Considering there are several approximations in this theoretical model, and the fact that we did not use any fitting parameter, the agreement observed in Figure 4 is generally very good.

We also performed simulations of the two most concentrated systems using a special version of Brownian Dynamics that does not rely on the definition of effective pair potentials or on any decoupling approximation to handle polydispersity. In these simulations, the 3D electric potential field is computed throughout the suspension at every time step by solving the non-linear Poisson-Boltzmann equation. The excess osmotic stress tensor field is then calculated and integrated to obtain the many-body forces exerted on every colloid accurately. [63][64][65][66] The structure factors obtained from these simulations and from the theoretical model described above are in very good agreement (see open circles in Figure 4). The slightly excessively marked peak in the theory for the sample at φ =0.042 is thus not due to the approximate treatment of polydispersity or a possible poor definition of the effective potential. As the experimental data seems close to that of a system with a smaller surface charge, a possible explanation could be an excessive acidification of the sample due to a prolonged exposition to ambient air.

Most importantly, the set of measurements on silica colloids and the comparison with theoretical and numerical models show that the meso-structure of interacting colloids can be measured quantitatively in situ and operando in a microfluidic chip, with a laboratory beamline. The photon flux is not really an issue here as silica is quite dense. The critical point for this kind of studies is the access to low wavenumbers since the value S(q → 0) is actually a direct measure of osmotic compressibility at least for monodisperse systems. Here we see that with the present setup where q min ≃ 0.005 Å -1 , it can be measured for volume fractions above 5%. Note that increasing the sample-to-detector distance would not be an option here because the divergence of the beam actually sets the lowest observable angle. Improvements of the beam size and divergence would be very interesting to reach lower q values presently accessible only on synchrotrons.

Macromolecule in solution

The BSA model protein was used to evaluate the possibility to obtain good scattering data from macromolecules with the present chip and beamline setup. The corresponding scattering curves are presented in Figure 5. No interactions are supposed to appear for the most dilute sample. The intensity recorded for this sample is compared to two models in the inset of Figure 5. The first one is the form factor obtained from the known atomic structure of crystallized BSA from the Protein Data Bank (3v03.pdb). This form Even for the most concentrated case, the structure factor presents only a weak, hardly visible peak around 0.08 Å -1 due to strong screening by the 0.15 mM PBS buffer. Interactions need to be accounted for, however, as the value S(q → 0) can be as low as about 0.5 for a BSA concentration of 100 g/L. Here, we found that using a unique value of Z = 16 unit charges per BSA molecule, with an effective radius a = 3.4 nm allowed to match the experimental data reasonably well. This is consistent with values found previously by similar fitting procedures, 68,69 and by ab initio calculations. 70 Note that the equation of state of BSA (osmotic pressure versus concentration) could also be predicted satisfactorily with these values. 70,71 The bump visible around q = 0.2 Å -1 especially for the most concentrated case is not an artifact. It cannot be reproduced with the ellipsoid model, but it is actually present in the structure factor extracted from the structure of crystallized BSA. Note that, however, using the latter form factor, the intensity at low q is significantly overestimated by the model. Determining the best coarse-graining strategy to model BSA is not the point of this work though, so we leave this issue open for future research.

Here we have provided, using BSA, a proof of concept highlighting the potential of the present mirofluidic/SAXS setup to study the structure and interactions of macromolecules in solution.

Latex nanoparticles

The three types of colloids mentioned above and used for benchmarking the setup all have a significant contrast with their solvent, thus allowing rather straightforward measurements. The last type of sample considered is an aqueous dispersion of nega-tively charged sulfate latex beads in a 1mM NaCl buffer. These colloids were chosen as they are well characterized, but also have a low electronic contrast with water so their SAXS signature is among the most difficult to obtain quantitatively. Data were acquired for volume fractions φ =0.035, 0.07, 0.105, and 0.14. The acquisition time was set to one hour for each sample. As expected, the scattering curves presented in Figure 6 contain a significant amount of noise compared to the samples investigated above. The first structure peak is well resolved for the three most concentrated samples and its location follows the expected φ 1/3 scaling due to long-range repulsions. The signal of the most dilute sample is quite difficult to fit with a form factor model because it is quite close to the background, and because of artifacts visible for q < 0.01 Å -1 . The latter are present because the slits were opened further to maximize the photon flux for this experiment. The same measurements were also performed in 1.5 mm diameter quartz capillaries in order to increase the signal to noise ratio, and the form factor corresponding to a polydisperse sphere model with radius 15 nm and polydispersity 16% was shown to fit the data quite well. This is consistent with the 13.5 nm average radius given by the manufacturer. In Figure 6, this form factor was used and structure has been calculated for an effective surface charge density σ =0.02 e/nm 2 . This value is about one tenth of the bare surface charge density measured by titration, which is reasonable. The use of latex nanospheres reveals the limits of laboratory SAXS measurements in microfluidic chips. The meso-structure of dispersions of small objects with low contrast such as polystyrene latex beads can be measured to some extent with acquisition times of the order of one hour but the poor signal-to-noise ratio makes attempts at modeling the data quite difficult. In order to improve this with the current setup, the pinhole could be positioned directly on the microfluidic chip, which would allow to open the slits further. However, since the pinhole is only at 2 cm from the sample now, this should not lead to a dramatic improvement of measurements. Another possible improvement of the setup would be to replace the solid anode source by a liquidmetal-jet or a rotating anode source, while keeping the use of the scatterless pinhole. This would probably allow an increase of at least one order of magnitude on the photon flux for the same small beam size compatible with on-chip measurements.

Conclusion

In summary, we established a proof-of-concept for a versatile microfluidic system coupled to a laboratory beamline for X-ray scattering experiments. The microfluidic chip was made of OS-TEMER and polyimide, with very thin walls, to maximize the signal-to-noise ratio by minimizing absorption and background scattering by the device itself. This system shows great potential to study various dispersions and was successfully used to obtain structural information on both inorganic nanoparticles and bio-macromolecules in suspension. The form factor and the first structure peak of concentrated dispersions of repulsive colloids could be measured accurately. The second peak was more difficult to distinguish from the background noise. The zero wavenumber limit of the structure factor could be observed in dense systems with the addition of a pinhole close to the sample. It allowed to open the slits of the beamline more than without pinhole, so the photon flux could be kept high enough despite the reduced beam size. A limitation of the device appeared for low-contrast materials such as polystyrene spheres suspended in water. The low quality of the data made quantitative fitting procedures difficult, although qualitative structural information could still be obtained. There is no doubt that using carefully chosen materials for microfluidic chips and state-of-the art laboratory beamlines with limited modifications now allows some in situ, operando SAXS measurements in the laboratory on a wide range of colloidal dispersions. The limited photon flux delivered by laboratory X-ray sources imposes significant acquisition times, ranging from seconds to tens of minutes. To perform time-resolved measurements in rapidly evolving systems, keeping a steady-state flow allows to associate a point in space to a time so long acquisitions can still be performed to observe fast transformations, although this may require large volumes of solutions. Time-resolved measurements on static samples can be undertaken for slowly evolving systems, as for example slow reactions, colloidal aggregation, rearrangements of the structure of metallic nanoparticles after synthesis, or processes involving slow solvent transfers such as osmotic compression 20 or solvent exchange. The simple setup investigated here opens the way to more sophisticated lab-on-chip devices in which in-situ scattering can be performed without the absolute need for a synchrotron facility.

Fig. 1

 1 Fig. 1 Top left panel: design of the microfluidic chip consisting in four parallel 400 µm wide, 3 cm long channels. Top right panel: picture of the final chip with inlet/outlet capillaries. Bottom panel: fabrication process of an OSTEMER chip: first the two OSTEMER stickers are prepared by injection and curing in hollow PDMS molds, then they are assembled together and fully cured. For SAXS optimal acquisition, in this set-up, X-rays pass through 25 µm of Kapton, 50 µm of OSTEMER, 400 µm of sample and 50 µm of OSTEMER

Fig. 2

 2 Fig. 2 Top: microfluidic chip on the Xeuss 2.0 SAXS beamline. The pinhole is located on the black motorized stage added to the standard setup. Bottom: schematic representation of the beamline showing the source emitting a wide beam (diameter about 1 mm), which is cut by slits 1 (0.7×0.7 mm 2 ) and 2 (0.3×0.3 mm 2 ), and by the pinhole (0.25×0.25 mm 2 ) located at 2 cm from the sample.

Fig. 3

 3 Fig. 3 Scattering data for gold nanoparticles and model for polydisperse spheres with radii following a Schulz distribution with average 2.75 nm and 17 % polydispersity. From top to bottom: dispersion as prepared, and diluted 2, 4, and 8 times.

to 5 ,

 5 12.5, 25 and 50 mg/mL solutions used in these experiments.

Fig. 4

 4 Fig. 4 Scattering data for Ludox HS40 silica nanoparticles. Main figure: structure factors shifted by steps of one for volume fractions φ = 0.00234, 0.0050, 0.00826, 0.0144, 0.0287, 0.0420 and 0.0907, from bottom to top. × and + symbols: experiments in two different chips. • symbls: simulations (see main text). Continuous lines: theory (see main text). Inset: intensity in dilute conditions (arbitrary units) and polydisperse sphere model for the form factor based on a Schulz distribution of radii with average 8 nm and 13 % polydispersity.

Fig. 5

 5 Fig.5Scattering data for different BSA concentrations: 100, 50, 25 and 12.5 mg/mL top to bottom (symbols), theory using the classical ellipsoid form factor (continuous lines) and using the form factor fitted from the BSA atomic structure (dashed lines). The inset shows data and theory for the most dilute sample: 5mg/mL.

Fig. 6

 6 Fig. 6 Scattering data for latex nanoparticles (symbols) and Hard-Sphere-Yukawa theory (continuous lines). From top to bottom: initial dispersion φ = 0.14 and dilution factors 3/4, 1/2, and 1/4. The overall quality of the data is lower for low-contrast material

| 1-8

1-8 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Acknowledgements

We thank J.-B. Salmon, I. Rodriguez-Ruiz, and S. Teychené for discussions concerning microfluidic chips development, G. Nägele and M. Heinen for discussions concerning the modelling of the structure of polydisperse systems and for lending their OZ code to us, and Sébastien Lledos, Simon Tricard, Lise-Marie Lacroix and Guillaume Viau for providing the gold nanoparticles. We acknowledge the ANR program Grant No. ANR-18-CE06-0021 for financial support. This work was performed using HPC resources from CALMIP (Grant 2023A-P21039).

Notes and references

Author Contributions

Conflicts of interest

There are no conflicts to declare.