Vianney Bruned

André Mas
email: andre.mas@umontpellier.fr

Sylvain Wlodarczyk

Weak convergence of particle swarm optimization

Keywords: Particle swarm optimization,, Convergence,, Central limit theorem

Particle swarm optimization algorithm is a stochastic meta-heuristic solving global optimization problems appreciated for its efficiency and simplicity. It consists in a swarm of interacting particles exploring a domain and searching a global optimum. The trajectory of the particles has been well-studied in a deterministic case. More recently authors shed some light on the stochastic approach to PSO, considering particles as random variables and studying them with probabilistic and statistical tools. These works paved the way to the present article. We focus here on weak convergence, the kind of stochastic convergence that appears in the Central Limit Theorem. We obtain three main results related to three different frameworks. These three settings depend either on the regime of the particles (oscillation or fast convergence) or on the sampling strategy (along the trajectory or in the swarm). The main application of these results is the construction of confidence intervals around the local optimum found by the swarm. The theorems are illustrated by a simulation study.

1 Introduction [START_REF] Eberhart | A new optimizer using particle swarm theory[END_REF] introduced the particle swarm optimization algorithm (PSO) based on social interactions (behaviors of birds). Since then PSO has known a great popularity in many domains and gave birth to many variants of the original algorithm (see [START_REF] Zhang | A comprehensive survey on particle swarm optimization algorithm and its applications[END_REF] for a survey of variants and applications). PSO is a stochastic meta heuristic solving an optimization problem without any evaluation of the gradient. The algorithm explores the search space in an intelligent way thanks to a population of particles interacting with each other and updated at each step their position and their velocity. The dynamic of the particles relies on two attractors: their personal best position (historical best position of the particle denoted p s n below), and the neighborhood best position (corresponding to the social component of the particles, denoted g s n ). In the dynamic equation, the attractors are linked with a stochastic process in order to explore the search space. Algorithm 1 refers to the classical version of PSO with S particles and N iterations.

Algorithm 1 Classical PSO

Initialize the swarm of S particles with random positions x s 0 and velocities v s 0 over the search space.

for n = 1 to N do Evaluate the optimization fitness function for each particle. Update p s n (personal best position) and g s n (neighborhood best position). Change velocity (v s n ) and position (x s n ) according to the dynamic equation. end for

The convergence and stability analysis of PSO are important matters. In the literature, there are two kinds of convergence:

the convergence of the particles towards a local or global optimum. This convergence is not obtained with the classical version of PSO. [START_REF] Van Den Bergh | A convergence proof for the particle swarm optimiser[END_REF] and [START_REF] Schmitt | Particle swarm optimization almost surely finds local optima[END_REF] proposed a modified version of PSO to obtain the convergence. the convergence of each particle to a point (e.g. [START_REF] Poli | Mean and variance of the sampling distribution of particle swarm optimizers during stagnation[END_REF]).

If we focus on the convergence of each particle to a point, a prerequisite is the stability of the trajectory of the particles. In a deterministic case, [START_REF] Clerc | The particle swarm-explosion, stability, and convergence in a multidimensional complex space[END_REF] dealt with the stability of the particles with some conditions on the parametrization of PSO. Later, [START_REF] Kadirkamanathan | Stability analysis of the particle dynamics in particle swarm optimizer[END_REF] used the Lyapunov stability theorem to study the stability. About the convergence of PSO, [START_REF] Van Den Bergh | A study of particle swarm optimization particle trajectories[END_REF] looked at the trajectories of the particles and proved that each particle converges to a stable point (deterministic analysis). Under stagnation hypotheses (no improvement of the personal and neighborhood best positions), [START_REF] Poli | Mean and variance of the sampling distribution of particle swarm optimizers during stagnation[END_REF] gives the exact formula of the second moment. More recently, [START_REF] Bonyadi | Stability analysis of the particle swarm optimization without stagnation assumption[END_REF] or [START_REF] Cleghorn | Particle swarm stability: a theoretical extension using the non-stagnate distribution assumption[END_REF] provided results for the order-1 and order-2 stabilities with respectively stagnant and non-stagnant distribution assumptions (both weaker than the stagnation hypotheses).

Let us introduce some notations. We consider here a cost function f : R d → R + that should be minimized on a compact set Ω. Consequently the particles evolve in Ω ⊂ R d .

Let x s n ∈ R d 1 ≤ s ≤ S denote the position of particle number s in the swarm at step n. Let (r j,n ) j=1,...,n≥1 be sequences of independent random vectors in R d whose margins obey a uniform distribution over [0, 1] and denote by ω, c 1 and c 2 three positive constants which will be discussed later. Then the PSO algorithm considered in the sequel is defined by the two following equations (or dynamic equation, [START_REF] Poli | Mean and variance of the sampling distribution of particle swarm optimizers during stagnation[END_REF]:

v s n+1 = ω • v s n + c 1 r 1,n (p s n -x s n ) + c 2 r 2,n (g s n -x s n ) , x s n+1 = x s n + v s n+1 .
(1)

where stands for the Hadamard product:

u v = (u 1 v 1 , ..., u d v d ) .
and p s n (resp. g s n ) is the best personal position (resp. the best neighborhood position of the particle s) :

p s n = argmin t∈{x s 0 ,...,x s n } f (t) , g s n = argmin t∈{p s n :s ∈V(s)} f (t) .
with V(s) the neighborhood of particle s. This neighborhood depends on the swarm's topology: if the topology is called global (all the particles communicate between each other) then g s n = g n = argmin t∈{p 1 n ,...,p S n } f (t) (see [START_REF] Lane | Particle swarm optimization with spatially meaningful neighbours[END_REF]).

Our main objective is to provide (asymptotic) confidence sets for the global or local optimum of the cost function f . If g = argmin t∈Ω f (t) for some domain Ω, a confidence region for g is a random set Λ such that :

P (g ∈ Λ) ≥ 1 -α
for some small α ∈ (0, 1). The set Λ depends on the swarm and is consequently random due to the random evolution of the particles. The probability symbol above, P, depends on the distribution of the particles in the swarm. Let us illustrate the use of confidence interval for a real-valued PSO. Typically the kind of results we expect is: g ∈ [m, M ] with a probability larger than, say .99 %. This does not aim at yielding a precise estimate for g but defines a "control area" for g, as well as a measure of the imprecision and variability of PSO.

Convergence of the swarm will not be an issue here. In fact we assume that the personal and global best converge : see assumptions A 2 and B 3 below. We are interested in the "step after" : localizing the limit of the particles with high probability, whatever their initialization and trajectories.

Formally, confidence set estimation forces us to inspect order two terms (i.e. the rate of convergence), typically convergence of the empirical variance. The word asymptotic just means that the sample size increases to infinity.

The outline of the paper is the following. In the next section the three main results are introduced. They are all related to weak convergence of random variables and vectors (see [START_REF] Billingsley | Convergence of probability measures[END_REF] for a classical monograph) and obtained under three different sets of assumptions.

The two first consider the trajectory of single particles. The sample consists in the path of a fixed particle. We show that two different regimes should be investigated depending on the limiting behavior of p n and g n . Briefly speaking : if the limits of p n and g n are distinct, the particles oscillate between them (which is a well-known characteristics of PSO), if the limits of p n and g n coincide, then particles converge at a fast, exponential, rate.

In the oscillating case a classical Central Limit Theorem is obtained relying essentially on martingale difference techniques. In the non-oscillating situation, the particle converges quickly and we have to use random matrices products to obtain a non-standard CLT. As by-products of these two subsections we will retrieve confidence sets of the form Λ (x s 1 . . . x s n ), depending on the n positions of each particle x s .

The third result states another classical CLT. The sample consists here in the whole swarm. This time the confidence set is of the form Λ x 1 n . . . x S n , depending on the S particles of the swarm when the iteration step is fixed to n. A numerical study and simulations are performed in a Python environment. A discussion follows. The derivations of the main theorems are collected in the last section.

Main results

The usual euclidean norm and associated inner product for vectors in R d are denoted respectively • and •, • . If X is a random vector with null expectation then E (X ⊗ X) = E XX t is the covariance matrix of X. The covariance matrix is crucial since it determines the limiting Gaussian distribution in the Central Limit Theorem. We will need two kinds of stochastic convergence in the sequel. convergence in probability of X n to X is denoted X n → P X. The arrow → stands for convergence in distribution (weak convergence).

Except in section 2.3 we consider a single particle in order to alleviate notations. We drop the particle index so that x s n = x n , p s n = p n and g s n = g n . In all the sequel we assume once and for all that :

c 1 = c 2 = c
Some authors took this choice for defining their working version of PSO, see for instance [START_REF] Clerc | Particle Swarm Optimization[END_REF] (see equation [3.2] p.39), [START_REF] Poli | Mean and variance of the sampling distribution of particle swarm optimizers during stagnation[END_REF].

Without this assumption, computations turn out to get intricate. Besides if c 1 = c 2 some constants become quite complicated and hard to interpret (see for instance condition A 2 and the two constants C and L involved in the asymptotic variance within Theorem 1 below).

At last we take for granted that particles are warm, reached an area of the domain were they fluctuate without exiting.

First case: oscillatory (p = g)

The following assumptions are required and discussed after the statement of Theorem 1. A 1 : For all n, x n ∈ Ω where Ω is a compact subset of R d . A 2 : The following holds :

1 √ N N n=1 (p n -Ep n ) → P 0, lim n→+∞ E p n -Ep n = 0, 1 √ N N n=1 (g n -Eg n ) → P 0, lim n→+∞ E g n -Eg n = 0.
A 3 : The inequality below connect c and ω :

0 < c < 12 1 -ω 2 7 -5ω .
Before stating the Theorem we need a last notation. Let δ = (δ 1 , ..., δ d ) ∈ R d . The notation diag (δ) stands for the diagonal d × d matrix with entries δ 1 , ..., δ d and δ 2 is the vector in R d defined by δ 2 = δ 2 1 , ..., δ 2 d .

Theorem 1 In addition to assumptions A 1-3 suppose that, when n goes to infinity Ep n → p and Eg n → g. Set:

L = 2c 1 -ω 1 + ω 1 + ω - c 2 - c 2 6 , C = c 12L 1 -ω 1 + ω 1 + ω - c 2 . Denote finally Γ = C • diag (p -g) 2 then: √ N 1 N N n=1 x n - 1 N N n=1 Ep n + Eg n 2 → N (0, Γ ) ,
where N (0, Γ ) denotes a Gaussian centered vector of R d with covariance matrix Γ .

Discussion of the Assumptions: We avoid here the assumption of stagnation: the personal and local best are not supposed to be constant but they oscillate around their expectation. The convergence occurs at a rate ensuring that neither g n nor p n are involved in the weak convergence of the particles x n . Condition A 2 is specific of what we intend by a convergent PSO. It ensures that p n or g n have no impact on the weak convergence behavior of the particles. With other words Assumption A 2 requires that the oscillations of p n and g n around their expectations are negligible. We tried here to model the stagnation phenomenon which consists in sequence of iterations during which g n (resp. p n ) remain constant hence g n = Eg n for n in N , N supposedly. Notice however that convergence of the expectation towards p and g is not mentioned at this step.

Note that assumption A 3 is exactly the condition found by Poli in Poli (2009) (see the last paragraph of section III) for defining order 2 stability. This condition may be extended to the case when c 1 = c 2 , see [START_REF] Cleghorn | Particle swarm stability: a theoretical extension using the non-stagnate distribution assumption[END_REF] and references therein. At last A 3 holds for the classical calibration appearing in [START_REF] Clerc | The particle swarm-explosion, stability, and convergence in a multidimensional complex space[END_REF] (constriction constraints) with c = 1.496172 and ω = 0.72984. The next Corollary simplifies Theorem 1 whenever the convergence of Ep n and Eg n to p and g holds with an appropriate rate. The centering term in the CLT is

then p+g 2 . Denote xN = 1 N N n=1 x n .
Corollary 1 If assumption A 1-3 hold and if :

1 √ N N n=1 (|p -Ep n | + |g -Eg n |) → 0, (2)
then Theorem 1 comes down to :

√ N xN - p + g 2 → N (0, Γ )
when N tends to infinity.

Remark 1 Condition 2 always holds under stagnation .

The Corollary above shows that the mean computed from the path of the particle converges to θ = (p + g) /2, losing the information on the location of p and g only. Note that the initial PSO version in 1 (with

c 1 = c 2 ) would lead to θ = (c 1 p + c 2 g) / (c 1 + c 2 ).
The next Corollary provides finally two kinds of by-products : asymptotic confidence intervals, in R, for the coordinates of θ and confidence regions for θ in R d . Let α ∈ [0, 1], q α be the 1 -α quantile of the standard Gaussian distribution and χ (d) 1-α be the 1 -α quantile of the Chi-square distribution with d degrees of freedom. Set also θ = (θ 1 , . . . , θ d ) and xN = (x N,1 , . . . , xN,d ).

Corollary 2 An asymptotic confidence interval at level 1 -α for θ is directly derived from Corollary 1 :

I 1-α (θ ) = [x N, -s (N, α) ; xN, + s (N, α)] with s (N, α) = |p -g | C
N q 1-α 2 . An asymptotic confidence region at level 1 -α for the vector θ = (p + g) /2 is :

Λ 1-α (θ) = t ∈ R d : N Γ -1/2 (t -xN ) 2 ≤ χ (d) 1-α = t = (t 1 , ..., t n ) : N C d =1 t -xN, p -g 2 ≤ χ (d) 1-α .
We note however that the vector θ may not be of crucial interest for the initial optimization problem conversely to g. This point will be discussed in the last section.

Remark 2 Assumption A2 deserves a last but technical remark. First, notice that 1 √ N N n=1 t n → P 0 does not imply lim n→+∞ E t n = 0 nor the converse. Take for instance t n = (-1) n t 0 with Et 0 = 0 then 1

√ N N n=1 t n → P 0 whereas E t n = E t 0 . Conversely take t n = t 0 / log (n) lim n→+∞ E t n = 0 but 1 √ N N n=1 t n = t 0 √ N N n=1 log -1 (n) cannot converge in probability to 0.
2.2 Second case: non-oscillatory and stagnant (p = g)

In this section we study again a single particle and suppose once and for all that x n ∈ R. We assume throughout this subsection that the particle is under stagnation that is p n = p for n sufficiently large (see assumption B 3 below). This assumption is strong but a more general framework leads to theoretical developments out of our scope. Starting from Equation (1), the PSO equation becomes this time:

x n+1 = (1 + ω) x n -ωx n-1 + c (r 1,n + r 2,n ) (p -x n ) .
Change the centering and consider x n -p = y n . The previous equation becomes :

y n+1 = (1 + ω -c + cε n ) y n -ωy n-1 , (3)
where ε n is the sum of two independent random variables with U [-1/2; 1/2] distribution.

Assuming that for all n y n = 0, we have then :

y n+1 y n = (1 + ω -c + cε n ) -ω y n-1 y n (4)
It is plain that y n+1 /y n defines a Markov chain (more precisely : a non-linear auto-regressive process) which will play a crucial role in the proofs and in the forthcoming results. Notice that N n=1 log

y n y n-1 = log |y N | -log |y 0 |.
We are ready to introduce a new set of assumptions. Denote π the stationary distribution of (y n /y n-1 ) n∈N and take (Z n ) n∈N a copy of (y n /y n-1 ) n∈N with Z 0 a realization of π. Then define :

µ x = E π log |Z 0 | , σ 2 x = Var π (log |Z 0 |) + 2 +∞ k=1 Cov π (log |Z 0 | , log |Z k |) . B 1 : x n -p x n-1 -p is a Harris recurrent Markov chain. B 2 : 1 + ω -c < ω/c < (1 + c)/4. B 3 : For sufficiently large n g n = p n = p = g is constant.
The definition of Harris recurrence needed above is in [START_REF] Meyn | Markov chains and stochastic stability[END_REF], beginning of Chapter 9. Theorem 2 Let ω ∈ (0, 1), c > 1, when B 1-3 hold, then:

1 √ n (log |x n -g n | -nµ x ) → N 0, σ 2 x .
Remark 3 The theorem above is not a Central Limit Theorem for x n . It is derived thanks to a CLT but it shows that the asymptotic distribution of |x n -g n | is asymptotically log-normal with approximate parameters nµ x and nσ 2 x .

Remark 4 The mean and variance µ x and σ 2 x are usually unknown but may be approximated numerically. We refer to the simulation section for more details.

Corollary 3 If p n = p an asymptotic confidence (non convex) region for g at level 1 -α denoted Λ 1-α below may be derived from the preceding Theorem:

Λ 1-α (g) = Λ + 1-α ∪ Λ - 1-α Λ + 1-α = x n + exp µ x + σ x √ n q α/2 , x n + exp µ x - σ x √ n q 1-α/2 Λ - 1-α = x n -exp µ x + σ x √ n q 1-α/2 , x n -exp µ x - σ x √ n q α/2
Remark 5 Under matrix form the equation 3 is purely linear but driven by a random matrix:

y n+1 y n = y n+1 = 1 + ω -c (1 + ε n ) -ω 1 0 y n y n-1
(5)

y n+1 =S n+1 y n , y n = S n S n-1 ...S 2 y 1 = T n y 1 , with T n = S n S n-1 ...S 2 .
It is plain here that a classical Central Limit Theorem cannot hold for the sequence (y n ) n∈N . We turn to asymptotic theory for the product of random matrices. We refer to the historical references: [START_REF] Furstenberg | Products of random matrices[END_REF] and [START_REF] Berger | Central limit theorem for products of random matrices[END_REF] who proved Central Limit Theorems for the regularity index of the product of i.i.d random matrices. Later [START_REF] Hennion | Limit theorems for products of positive random matrices[END_REF] generalized their results. But the assumptions of (almost surely) positive entries is common to all these papers. Other authors obtain similar results under different sets of assumptions (see Le Page, 1982, Benoist and[START_REF] Benoist | Central limit theorem for linear groups[END_REF], and references therein), typically revolving around characterization of the semi-group spanned by the distribution of S. These assumptions are uneasy to check here and we carried out to a direct approach with Markov chain fundamental tools.

The swarm at a fixed step

In this section we change our viewpoint. Instead of considering a single particle and sampling along its trajectory we will take advantage of the whole swarm but at a fixed and common iteration step. Our aim here is to localize the minimum of the cost function based on x 1 n , . . . , x S n . This time the particle index s varies up to S the swarm size, whereas the index n is fixed. In this subsection we assume that S ↑ +∞ and asymptotic is with respect to S. We do not drop n in order to see how the iteration steps influence the results. We still address only the case x s n ∈ R even if our results may be straigthforwardly generalized to x s n ∈ R d . We provide below a Central Limit Theorem suited to the case when the number of particles in the swarm becomes large. In order to clarify the method, we assume that for all particles x i n in the swarm p i n = g n = p. In other words, no local minimum stands in the domain Ω, which implies additional smoothness or convexity assumptions on the cost function f . This may be possible by a preliminary screening of the search space. Indeed a first (or several) run(s) of preliminary PSO(s) on the whole domain identifies an area where a single optimum lies. Then a new PSO is launched with initial values close to this optimum and with parameters ensuring that most of the particles will stay in the identified area.

So we are given x 1 n , ..., x S n where S is the sample size. Basically, the framework is the same as in the non oscillatory case studied above for a single particle. From (5) we get with y i n = x i n -p :

y i n = T i n y i 1 , T i n = Π n j=2 S j , S j = 1 + ω -c 1 + ε i j -ω 1 0 .
Assume that the domain Ω contains 0 and that for all s x i 0 , x i 1 i≤S are independent, identically distributed and centered then from the decomposition above, for all n and s, Ey i n = 0 and the y i n 1≤i≤S are i.i.d too. The assumptions we need to derive Theorem 3 below are : C 1 : The operational domain Ω contains 0 (and is ideally a symmetric set). C 2 : The couples x i 0 , x i 1 i≤S are i.i.d. and centered. C 3 : For all i in {1, ..., S} p i n = g n = p. When S is large the following Theorem may be of interest and is a simple consequence of the i.i.d. CLT.

Theorem 3 Under assumptions C 1-3 a CLT holds when S the number of particles in the swarm tends to +∞ :

1 √ S S i=1 x i n -g n → S→+∞ N 0, σ 2 n ,
where σ 2 n = E x 1 n -g n 2 is estimated consistently by :

σ 2 n = 1 S S i=1 x i n -g n 2 .
Remark 6 The convergence of σ 2 n to σ 2 n is a straightforward consequence of the weak and strong laws of large numbers. Denote x S = (1/S) S i=1 x i n . The Theorem above paves the way towards an asymptotic confidence interval.

Corollary 4 An asymptotic confidence interval at level 1 -α for g is :

Λ n (g) = x S - σ n √ S q 1-α/2 , x S + σ n √ S q 1-α/2 .

Simulation and numerical results

The Himmelblau's function is chosen as example for our experiments. It is a 2 dimensional function with four local optima in [-10, 10] 2 defined by : f (x, y) = x 2 + y -11 2 + x + y 2 -7 2 . Figure 3 illustrates the contour of this function. With the Himmelblau's function, we can observe the two different behaviors of the particles: oscillatory and non-oscillatory. The Himmelblau's function is positive and has four global minima in: (3, 2), (-2.81, 3.13), (-3.77, -3.28), (3.58, -1.84) where f (x, y) = 0. We use a ring topology (for a quick review of the different topologies of PSO see [START_REF] Lane | Particle swarm optimization with spatially meaningful neighbours[END_REF] for the algorithm in order to have both oscillating and non-oscillating particles. The latter converge quickly. The former go on running between two groups of particles converging to two distinct local optima.

Oscillatory case

We select particles oscillating between (3.58, -1.84) and (3, 2), these values could be both their personal best position or their neighborhood best position. In this case, the convergence of the p n and g n to (3.58, -1.84) or (3, 2) are satisfying the conditions of Corollary 1. We have to verify the Gaussian asymptotic behavior of

H s 1 (N ) = √ N 1 N N n=1 x s n -p+g 2
for each s oscillating particle. We launch PSO with a population of 200 particles and with 2000 iterations, ω = 0.72984 and c = 1.496172. A ring topology was used to ensure the presence of oscillating particles. A particle is said oscillating if between the 500th and the 2000th iteration, Assumptions A 1-3 holds.

A visual tool to verify the normality of H s 1 (N ) for a particle is a normal probability plot. Figures 4 and5 displays the normal probability plot of H s 1 (N ) respectively for the x axis and y axis. For each axis, the normality is confirmed: H s 1 (N ) fits well the theoretical quantiles. To check the formula of the covariance matrix Γ , the confidence ellipsoid is also a good indicator (see Figure 6). For a single particle, H s 1 (N ) is not necessarily always inside the confidence ellipsoid and does not respect the percentage of the defined confidence level. Figure 7 shows the trajectory of x s n and H s 1 (N ) on the y axis, H s 1 (N ) remains bounded. With 200 Monte-Carlo simulations of PSO (200 particles, and 2000 iterations), we select all the particles oscillating between (3.58, -1.84) and (3, 2), and for each of them we compute H s 1 (2000). Figure 8 displayed the density of H s 1 (2000) using 1150 oscillating particles. Almost 95% of the particles are inside the confidence ellipsoid of level 95% (represented in red).

Non-oscillatory case

We study now the behaviors of non-oscillating particles on the Himmelblau's function. We launch PSO with a population of 1000 particles and with 2000 iterations, ω = 0.72984 and c = 1.496172. A ring topology was used to ensure the presence of enough particle converging to each local optimum. We select particles converging to (3, 2), meaning that p n = g n = p for a sufficiently large n. For the weak convergence of the particle, we consider:

H s 2 (N ) = 1 √ N (log |x s N -g N | -N µ x ) .
First, it is easy to check the linear dependency of log |x s N -g N | with a single display of the trajectory. Figure 9 illustrates this phenomenon for a single particle. We observed numerical issues when we reach the machine precision, but a numerical approximation of µ x can be performed thanks to a linear regression. Fig. 9 |xn -gn| over 500 iterations in a logarithmic scale. After 300 iterations, we reach the computer precision. We take advantage of the linear behavior of log (|xn -gn|) on the 200 first iterations to perform a linear regression estimating µx.

Using many converging particles, a Monte Carlo approximation of µ x is done. For the approximation of σ x , a possibility is:

σ2 x = Var π (log |X 0 |) + 2 T k=1 Cov π (log |X 0 | , log |X k |) ,
where T = 20. With near 240 converging particles to (3, 2), we found that for the first coordinate: μx = -0.032, σx = 0.156.

We verify the asymptotic normality of H 2 (N ) with a normal probability plot using the approximation of µ x . Figure 10 displays the normal probability plot of H 2 (N ) on the first coordinate, the theoretical quantiles are well fitted by H 2 (N ). Fig. 11 Trajectories of H 2 (N ) on the first coordinate for five particles. The red dot represents the 95% confidence interval deduced from σx. Trajectories stop around the 400 iterations (after an heating phase) due to numerical precision.

Swarm at a fixed step

To check the Theorem 3, we study:

H n 3 (S) = 1 √ S σ n S i=1 x i n -g n .
In practice, we encountered some difficulties to verify Theorem 3 because of the convergence rate of the particles. Indeed, when p s n = g s n = p s , the particle s converges exponentially to g s n but the spread of the rate of convergence is large. As a consequence, at a fixed step of PSO, some particles could be considered as outliers because of a lower rate of convergence. Because of these particles qualified as belated, the asymptotic Gaussian behavior of H n 3 (S) is not verified. A solution is to filter the particles and remove the belated particles. Figures 12 and13 illustrate this phenomenon for the Himmelblau's function in 2D and with near 1500 converging particles to (3, 2) over 500 iterations. In Figure 12, we compute without any filtering H n 3 (S) and we notice that the Gaussian behavior is not verified and some jumps appeared. The presence of these "jumps" is due to belated particles which have a lower rate of convergence in comparison to the swarm. When we remove these particles with a classical outliers detection algorithm in Figure 13, Theorem 3 seems to be verified. (S) at the 200th iteration on the first coordinate, using all the particles. We observe a discontinuity of the probability plot due to belated particles.

Discussion

Our main theoretical contribution revolves around the three CLTs and the confidence regions derived from sampling either a particle's path or the whole swarm. Practically the confidence set Λ 1-α localizes, say g, with a high probability. High probability means larger than 1 -α. The simulations carried out in Python language tend to foster the theory.

This work was initiated in order to solve a practical issue in connection with oil industry and petrophysics. Yet in the previous section we confined ourselves to simulated data for several reasons. It appears that our method should be applied on real-data for a clear validation.

A second limitation of our work is the asymptotic set-up. The results are stated for samples large enough. This may not have an obvious meaning for people not familiar with statistics or probability. Practitioners usually claim that the Central Limit Theorem machinery works well for a sample size larger than thirty/forty whenever data are stable, say stationary (path with no jumps...). As a consequence the first iterations of the algorithm should be avoided to ensure a warming. When studying PSO the behavior of p n and g n turns out to be crucial too in order to ensure this stability, hence the validity of the theoretical results. However the control of p n and g n is a difficult matter, which explains our current assumptions on stagnation or "almost stagnation". However, in order to fix the question of the asymptotic versus non-asymptotic approach a work is on progress. We expect finite-sample concentration inequalities : the non-asymptotic counterparts of CLT.

In the oscillating framework (see Corollary 1) our results involve the parameter denoted θ, the center of the segment [p, g]. Here we miss the target g for instance.

We claim that our method may be adapted to tackle this issue. Briefly speaking : running two independent PSOs on the same domain with two distinct pairs (c 1 , c 2 ) and c 1 , c 2 will provide, under additional assumptions on the local and global minima of the cost function, two CLTs involving :

θ = c 1 g + c 2 p c 1 + c 2 and θ = c 1 g + c 2 p c 1 + c 2 .
A simple matrix inversion then recovers the CLT for e.g. g.

Derivations of the results

We start with some notations. First we recall that the sup-norm for square matrices of size d defined by:

M ∞ = sup x =0 M x x .
The tensor product notation is appropriate when dealing with special kind of matrices for instance covariance matrices. Let u and v be two vectors of R d then u ⊗ v = uv t (where v t is the transpose of vector v) stands for the rank one matrix defined for all vector x by (u ⊗ v) (x) = v, x u. Besides u ⊗ v ∞ = u v . The Hadamard product between vectors was mentioned earlier. Its matrix version may be defined a similar way. Let M and S be two matrices with same size then M S is the matrix whose (i, j) cell is (M S) i,j = m i,j s i,j . We recall without proof the following computation rule mixing Hadamard and tensor product. Let η, ε, u and v be four vectors in R d . Then:

(η u) ⊗ (ε v) = (η ⊗ ε) (u ⊗ v) , (6)
and the reader must be aware that the Hadamard product on the left-hand side operates between vectors whereas on the right-hand side it operates on matrices.

We will need F s n the filtration generated by the path of particle number s up to step n : {x s 0 , ..., x s n } and F S n the filtration generated by the swarm up to step n : {x s 0 , ..., x s n : s = 1, ..., S} . We will also denote later g n = E (g n )+ξ n and p n = E (p n )+ν n the expectationvariance decomposition of g n and p n where ξ n and ν n are centered random vectors and support all the variability of g n and p n respectively.

First case: oscillatory

The proof of Theorem 1 is developed in the next two subsections. The first subsection below provides some short preliminary results and the proof of the Theorem. The second subsection is devoted to the -somewhat long-derivation of the single Proposition 2. This Proposition 2 is needed within the proof of Theorem 1.

Preliminary results and proof of Theorem 1

Let us start with some Lemmas who will be invoked later.

Lemma 1 Let ε (i) n and η (j)
n be any coordinate of the random vectors ε n and η n . Clearly, ε (i) n and η (j) n are not independent but not correlated and follow the same type of distribution. Besides:

Eε (i) 2
= Eη (j) 2 = 1/6, Eε (i) η (j) = 0,

Eε (i) 3 η (j) = Eε (i) η (j) 3 = 0.
The proof is very simple hence omitted.

Lemma 2 When A 1-3 hold and when N goes to infinity: 7) we get in matrix form:

1 N N n=1 Ep n + Eg n 2 -Ex n 2 → 0. Proof: Denote q n = Ex n -(Ep n + Eg n ) /2. From (
q n+1 q n = 1 + ω -c -ω 1 0 q n q n-1 + K n 0
where:

K n+1 = (1 + ω) E (p n +g n ) -ωE (p n-1 +g n-1 ) -E (p n+1 +g n+1 ) .
The above equation may be rewritten q n+1 = Mq n + K n+1 with obvious but shorter notations. Notice that in

K n each E (p k +g k ), k ∈ {n -1, n, n + 1} may be replaced by E (p k +g k ) -(p + g).
Consequently, we know that K n tends to 0. Also notice that M ∞ < 1 because ω < 1.

Then we derive:

q n = n p=2 M n-p K p .
It suffices to notice now that if q n decays to 0, so does q n and Lemma 2 will be derived by a simple application of Cesaro's famous Lemma. Let us prove now that q n ↓ 0.

Take ε > 0 denote K ∞ = max ( K p ) and

C M = 1 -M ∞ -1 . First pick N * such that sup p≥N * +1 K p < ε/ (2C M ) q n ≤ N * p=2 M n-p ∞ K p + n p=N * +1 M n-p ∞ K p ≤ M n-N * ∞ N * p=2 M n-p ∞ K p + max p≥N * +1 K p n p=N * +1 M n-p ∞ ≤ C M K ∞ M n-N * ∞ + sup p≥N * +1 K p . Then let N † be such that M N † ∞ < ε/ (2C M K ∞ ).
It can be seen from the equations above that for n > N † + N * , q n < ε.

The following Proposition is crucial:

Proposition 1 Let S N = N n=1 z n then:

S N = N n=1 η n 2 (p n -g n ) -ε n z n + N n=1 rn + 1 c (z 1 -ωz 0 + ωz N -z N +1 ) . Under assumptions A 2 1 √ N N n=1
rn → P 0.

Besides when A 1 holds (the x n s are almost surely bounded) then S N / √ N converges weakly if and only if:

1 √ N N n=1 η n 2 (p n -g n ) -ε n z n converges weakly too.
Proof of Proposition 1: The first equation is obtained by summing from n = 1 to N both sides of Equation ( 8). Then we are going to prove that 1 √ N N n=1 rn → P 0 by Markov inequality which will be enough to get the rest of the Proposition. We start from rn = ε n

Ep n +Eg n 2 -Ex n + 1+ε n 2 (p n -Ep n + g n -Eg n ). By the first part of assumption A 2 1 √ N N n=1 (p n -Ep n ) → P 0 1 √ N N n=1 (g n -Eg n ) → P 0.
We turn to:

E N n=1 ε n Ep n + Eg n 2 -Ex n 2 = Eε 2 1 N n=1 Ep n + Eg n 2 -Ex n 2 because ε n is an i.i.d. centered sequence.
The proposition now follows from Lemma 2. This completes the proof of the Proposition.

Proposition 1 shows that the limiting distribution of S N is completely determined by the limiting distribution of

N n=1 η n 2 (p n -g n ) -ε n z n . If a
Central Limit Theorem holds for the previous series, then the limiting distribution will depend on the covariance matrix of η n 2 (p n -g n ) -ε n z n . The latter will be decomposed in several terms.

Proof of Theorem 1: As a preliminary step, let us rewrite (1) on a single line :

x n+1 = (1 + ω) x n -ωx n-1 +c (r 1,n + r 2,n ) p n + g n 2 -x n +c (r 1,n -r 2,n ) p n -g n 2 .
Taking expectation:

Ex n+1 = (1 + ω) Ex n -ωEx n-1 + c Ep n + Eg n 2 -Ex n . (7)
By subtracting, denoting the centered process z n = x n -Ex n and e the vector of R d defined by e = (1, 1, ..., 1) the PSO equation becomes :

z n+1 = (1 + ω -c) z n -ωz n-1 -cε n z n + cr n + c 2 η n (p n -g n ) , (8)
with

ε n = r 1,n + r 2,n -e η n = r 1,n -r 2,n , rn = ε n Ep n + Eg n 2 -Ex n + e + ε n 2 (p n -Ep n + g n -Eg n ) .
Equation ( 8) will be a starting point especially for studying single particle trajectories.

We turn to the Theorem. From:

1 N N n=1 x n - 1 N N n=1 Ep n + Eg n 2 = 1 N N n=1 z n + 1 N N n=1 Ex n - Ep n + Eg n 2 , we see from the proof of Lemma 2 that 1/ √ N N n=1 [Ex n -(Ep n + Eg n ) /2
] → 0 as N tends to infinity. Following Proposition 1, the Theorem 1 comes down to proving:

1

√ N N n=1 η n 2 (p n -g n ) -ε n z n → N (0, Γ ) .
We can first remark that u n = η n 2 (p n -g n )-ε n z n is a vector-valued martingale difference sequence with respect to the filtration F s n = F n . We confine ourselves to proving a Levy-Lindeberg version of the CLT for the series of u n in two steps (Theorem 2.1.9 p. 46 and its corollary 2.1.10 in [START_REF] Duflo | Random Iterative Models[END_REF]: first ensuring convergence of the conditional covariance structure of 1/ √ N N n=1 u n , then checking the Lyapunov condition holds (hence the Lindeberg's uniform integrability that ensures uniform tightness of the sequence).

First step:

The conditional covariance sequence of u n is the sequence of matrices defined by 1

N N n=1 E (u n ⊗ u n |F n ).
We show in this first step that this sequence converges in probability to Γ . We start with elementary calculations:

1 N N n=1 E (u n ⊗ u n |F n ) = 1 N N n=1 E η n 2 (p n -g n ) ⊗ η n 2 (p n -g n ) |F n + 1 N N n=1 E ((ε n z n ) ⊗ (ε n z n ) |F n ) - 1 2N N n=1 E ([η n (p n -g n )] ⊗ [ε n z n ] |F n ) - 1 2N N n=1 E ([ε n z n ] ⊗ [η n (p n -g n )] |F n ) = Eε 2 1 N I N n=1 1 4 [(p n -g n ) ⊗ (p n -g n )] + z n ⊗ z n = 1 6 I 1 4 Γ N,p,g + Γ N,z ,
where we used Lemma 1 and denoted Γ N,p,g = 1

N N n=1 [(p n -g n ) ⊗ (p n -g n )] and : Γ N,z = 1 N N n=1 z n ⊗ z n ,
By Proposition 2 stated and proved in the next subsection, Γ N,z -H N Γ N,p,g tends in probability to zero. Hence whenever the limit in the r.h.s below exists:

lim N →+∞ 1 N N n=1 E (u n ⊗ u n |F n ) = 1 24 + H 6 lim N →+∞ I Γ N,p,g .
Under assumption A 1-2 lim N →+∞ Γ N,p,g -Γ N,p,g = 0 in matrix norm, where

Γ N,p,g = 1 N N n=1 E [(p n -g n ) ⊗ E (p n -g n )]
.

To prove this it suffices to write:

Γ N,p,g -Γ N,p,g = 1 N N n=1 (ν n -ξ n ) ⊗ (p n -g n ) + 1 N N n=1 E (p n -g n ) ⊗ (ν n -ξ n ) ,
and apply the bounds derived from the proof of Proposition 2, for instance:

1 N N n=1 (ν n -ξ n ) ⊗ (p n -g n ) ∞ ≤ 1 N N n=1 (ν n -ξ n ) ⊗ (p n -g n ) ∞ ≤ 1 N N n=1 ν n -ξ n p n -g n ≤ sup n p n -g n 1 N N n=1 ( ν n + ξ n ) → P 0.
Finally, the convergence of Γ N,p,g to Γ is a consequence of the adaptation of Cesaro's Lemma for tensor products.

Second step: A Lyapunov condition holds almost trivially here. Namely we are going to prove that:

1 N 2 N n=1 E u n 4 → N →+∞ 0.
Simple calculations provide:

u n 2 ≤ 2 η n 2 (p n -g n ) 2 + 2 ε n z n 2 ≤ 1 2 sup n,i {(p n,i -g n,i ) 2 } η n 2 + 2 sup n,i {z 2 n,i } ε n 2 ≤ 2 |Ω| 2 ∞ η n 2 + ε n 2 , with |Ω| ∞ = sup s∈Ω |s| hence E u n 4 ≤ C |Ω| 4 ∞ ,
for some constant C then 1/N 2 n≤N E u n 4 tends to zero when N tends to infinity. This completes the proof of Theorem 1.

Proof of Proposition 2

We are ready to state and prove this important intermediate Proposition 2 Under assumption A 3 denote:

H = c 2 24 2c 1 -ω 1 + ω 1 + ω - c 2 - c 2 6 -1 , then when N tends to +∞ Γ N,z - H N N n=1 diag δ 2 n → P 0, with δ n = p n -g n .
The proof of the Proposition will immediately follow from this Lemma.

Lemma 3 Let E n be a sequence of i.i.d centered random matrices with finite moment of order 4, let u n and v n two sequence of random vectors almost surely bounded and such that (u n , v n ) is for all n independent from E n then for the

• ∞ norm: 1 N N n=1 E n (u n ⊗ v n ) → P 0.
The proof of the above Lemma is a straightforward application of Cauchy-Schwartz inequality.

Proof of the Proposition 2: We take advantage of Equation ( 8) and note that

(ε n z n ) ⊗ (ε n z n ) = (ε n ⊗ ε n ) (z n ⊗ z n ) = (ε n ⊗ ε n ) -σ 2 ε I (z n ⊗ z n ) + σ 2 ε I (z n ⊗ z n ) ,
where σ 2 ε = 1/6 by Lemma 1 and I stands for identity matrix. After some tedious calculations we obtain:

z n+1 ⊗ z n+1 = (1 + ω -c) 2 (z n ⊗ z n ) + ω 2 z n-1 ⊗ z n-1 (9) + c 2 σ 2 ε [I (z n ⊗ z n )] + c 2 4 [η n (p n -g n )] ⊗ [η n (p n -g n )] -ω (1 + ω -c) [z n ⊗ z n-1 + z n-1 ⊗ z n ] + C n .
For the sake of completeness, we list now all the eleven terms contained in C n . In order to simplify notations let [[u : v]

] = u ⊗ v + v ⊗ u. Notice for further use: [[u : v]] ∞ ≤ 2 u v .
Then:

C n = c (1 + ω -c) [[z n : rn ]] -[[z n : ε n z n ]] + 1 2 [[z n : η n (p n -g n )]] (10) + ωc [[z n-1 : ε n z n ]] -[[z n-1 : rn ]] - 1 2 [[z n-1 : η n (p n -g n )]] -c 2 [[cε n z n : rn ]] + 1 2 [[ε n z n : η n (p n -g n )]] - 1 2 [[r n : η n (p n -g n )]] + c 2 rn ⊗ rn + c 2 (ε n ⊗ ε n ) -σ 2 ε I (z n ⊗ z n ) .
From 9 we sum over all indices n in {1, ..., N } we get:

1

-(1 + ω -c) 2 -ω 2 N n=1 z n ⊗ z n = c 2 σ 2 ε N n=1 [I (z n ⊗ z n )] + c 2 4 N n=1 [η n (p n -g n )] ⊗ [η n (p n -g n )] -ω (1 + ω -c) N n=1 [z n ⊗ z n-1 + z n-1 ⊗ z n ] + N n=1 C n + z 1 ⊗ z 1 -z N +1 ⊗ z N +1 + ω 2 (z 0 ⊗ z 0 -z N ⊗ z N ) .
Now we need to go slightly further in the computations and to find a recurrent equation for z n ⊗ z n-1 . Let us start again from (8):

z n+1 ⊗ z n = (1 + ω -c) z n ⊗ z n -ωz n-1 ⊗ z n -c (ε n z n ) ⊗ z n + cr n ⊗ z n + c 2 [η n (p n -g n )] ⊗ z n , z n ⊗ z n+1 = (1 + ω -c) z n ⊗ z n -ωz n ⊗ z n-1 -cz n ⊗ (ε n z n ) + cz n ⊗ rn + c 2 z n ⊗ [η n (p n -g n )]
.

Summing over all indices n in {1, ..., N } above we come to

(1 + ω) N n=1 [z n ⊗ z n-1 + z n-1 ⊗ z n ] = z 1 ⊗ z 0 + z 0 ⊗ z 1 -z N ⊗ z N +1 -z N +1 ⊗ z N + 2 (1 + ω -c) N n=1 [z n ⊗ z n ] + c N n=1 [[r n : z n ]] -[[(ε n z n ) : z n ]] + 1 2 [[η n (p n -g n ) : z n ]] .
Plugging the last equation in line ( 9) we get finally:

1

-(1 + ω -c) 2 -ω 2 N n=1 z n ⊗ z n = c 2 σ 2 ε N n=1 [I (z n ⊗ z n )] + c 2 4 N n=1 [η n (p n -g n )] ⊗ [η n (p n -g n )] -2ω (1 + ω -c) 2 1 + ω N n=1 [z n ⊗ z n ] + ωc (1 + ω -c) 1 + ω N n=1 [[(ε n z n ) : z n ]] -ωc (1 + ω -c) 1 + ω N n=1 [[r n : z n ]] -ωc (1 + ω -c) (1 + ω) N n=1 1 2 [[η n (p n -g n ) : z n ]] + N n=1 C n . Denoting κ = 1 -(1 + ω -c) 2 -ω 2 + 2ω (1+ω-c) 2 1+ω = 2c 1-ω 1+ω 1 + ω -c 2 : κ N n=1 z n ⊗z n = c 2 6 N n=1 [I (z n ⊗ z n )]+ c 2 4 N n=1 [η n (p n -g n )]⊗[η n (p n -g n )]+R N , with R N = -ω (1 + ω -c) 1 + ω ([[z 1 : z 0 ]] -[[z N : z N +1 ]]) + z 1 ⊗ z 1 -z N +1 ⊗ z N +1 (11) + ω 2 (z 0 ⊗ z 0 -z N ⊗ z N ) + N n=1 C n + ωc (1 + ω -c) 1 + ω N n=1 [[(ε n z n ) : z n ]] -[[r n : z n ]] - 1 2 [[η n (p n -g n ) : z n ]] .
It is proved in Lemma 4 that R N /N → P 0. First, let us unravel the matrix:

T N = κ N N n=1 z n ⊗ z n - c 2 6 1 N N n=1 [I (z n ⊗ z n )] .
Denote [T N ] ij the (i, j) cell of matrix T n . It is simple to see that:

[T N ] ii = κ - c 2 6 1 N N n=1 (z n ⊗ z n ) ii , [T N ] ij = κ 1 N N n=1 z n ⊗ z n i,j (12)
for i = j. Now denote:

Λ N = 1 N N n=1 [η n (p n -g n )] ⊗ [η n (p n -g n )] = 1 N N n=1 (η n ⊗ η n ) [(p n -g n ) ⊗ (p n -g n )]
.

Taking conditional expectation w.r.t F n we get -all non diagonal term vanish:

E (Λ N |F n ) = 1 N N n=1 [E (η n ⊗ η n )] [(p n -g n ) ⊗ (p n -g n )] = Eη 2 1 N N n=1 I [(p n -g n ) ⊗ (p n -g n )] ,
with I the identity matrix. Noting that the difference:

Λ N -E (Λ N |F n ) = 1 N N n=1 (η n ⊗ η n ) -Eη 2 1 I [(p n -g n ) ⊗ (p n -g n )]
vanishes when N tends to infinity by applying Lemma 3 we get with Landau notation in probability:

T N = c 2 24N N n=1 I [(p n -g n ) ⊗ (p n -g n )] + o P 1 N .
From ( 12) we obtain simultaneously that [T N ] ij → P 0 for i = j and

[T N ] ii - c 2 24 κ -c 2 6 1 N N n=1 I [(p n -g n ) ⊗ (p n -g n )] → P 0
which is precisely the statement of Proposition 2.

Lemma 4 R N /N → P 0 in matrix norm.

Proof: Let us start from (11). First by assumption A 1 :

1

N ([[z 1 : z 0 ]] -[[z N : z N +1 ]])+z 1 ⊗z 1 -z N +1 ⊗z N +1 +ω 2 (z 0 ⊗ z 0 -z N ⊗ z N ) → P 0.
It must be also noticed that:

ωc (1 + ω -c) 1 + ω N n=1 [[(ε n z n ) : z n ]] -[[r n : z n ]] - 1 2 [[η n (p n -g n ) : z n ]]
involves three terms between double brackets that already appear in C n , up to constants. Proving that each term of (1/N ) N n=1 C n tends in probability to 0 will be sufficient to complete the proof of the Lemma.

Our aim is to perform successive applications of Lemma 3. For the sake of completeness we remind now the eleven terms contained in C n and mentioned earlier.

C n = c (1 + ω -c) [[z n : rn ]] -[[z n : ε n z n ]] + 1 2 [[z n : η n (p n -g n )]] + ωc [[z n-1 : ε n z n ]] -[[z n-1 : rn ]] - 1 2 [[z n-1 : η n (p n -g n )]] -c 2 [[cε n z n : rn ]] + 1 2 [[ε n z n : η n (p n -g n )]] - 1 2 [[r n : η n (p n -g n )]] + c 2 rn ⊗ rn + c 2 (ε n ⊗ ε n ) -σ 2 ε I (z n ⊗ z n ) .
In the list above terms numbered 2, 3, 4, 6, 11 vanish by applying successively (6) and Lemma 3. Take for instance term 1, get rid of the constants and focus on the first tensor product in the double bracket namely:

z n ⊗ (ε n z n ) = (e z n ) ⊗ (ε n z n ) = (e ⊗ ε n ) (z n ⊗ z n ) .
Recall that e is the vector with all components valued at 1. Then we can apply Lemma 3 which shows that (1/N ) N n=1 (e ⊗ ε n ) (z n ⊗ z n ) vanishes in probability. It is not hard to see that terms 3, 4, 6 and 11 may be treated the same way.

Term number 8 namely depends on both ε n and η n . An application of (6) leads to considering, up to a constant and commutation to:

(η n ⊗ ε n ) (z n ⊗ (p n -g n )) .
The reader can check that the matrix E n = (η n ⊗ ε n ) is stochastically independent from z n ⊗ (p n -g n ), that the sequence of matrices E n are independent and centered. Lemma 3 applies here.

The terms numbered 1, 5, 7, 9, 10 depend on rn . Now consider

rn = ε n Ep n + Eg n 2 -Ex n + e + ε n 2 (ξ n + ν n ) = r n1 + r n2 ,
and the five terms containing them in C n . We split rn and deal separately with r n1 and r n2 to derive the needed bounds in probability.

The strictly positive K

1 (ω, c, τ ) exists because √ x log 2 |(ω/x) -λ + c| is bounded on [0, A τ ]. The first condition reads: 1 2c (ω/x)-λ+c (ω/x)-λ-c g (s) ds ≤ ρ 1 C 1 / |x|, 0 ≤ x ≤ A τ whenever C 2 K 1 (ω, c, τ ) ≤ ρ 1 C 1 (15)
and ρ 1 will be fixed after the second step.

Second step: Now we turn to A τ ≤ x ≤ M τ . We still have g (x) = C 1 / |x| but we need to focus on the bounds of the integral.

This time the lower bound of the integral

(ω/x) -λ -c ∈ [-λ -τ, M τ ] and the upper bound (ω/x) -λ + c ∈ [2c -λ -τ, 2c + M τ ].
We are going to require that (ω/x) -λ -c ≥ -M τ it suffices to take λ + τ ≤ M τ and this comes down to the following set of constraint on τ : {τ ≥ c -ω} ∪ {τ ≤ c -1} . We keep the second and assume once and for all that:

τ ≤ c -1. ( 16
) Then for x ∈ [A τ , M τ ], 1 2c (ω/x)-λ+c (ω/x)-λ-c g (s) ds = 1 2c M τ (ω/x)-λ-c g (s) ds + 1 2c (ω/x)-λ+c M τ g (s) ds (17) ≡ I 1 + I 2 .
We want to make sure that the upper bound (ω/x) -λ + c is larger than M τ . This will hold if

(ω/M τ ) -λ + c ≥ M τ hence if 2c -τ -λ ≥ M τ . We imposed previously that λ + τ ≤ M τ . So: τ < c - ω c ⇒ M τ < c ⇒ 2c -τ -λ ≥ M τ
but the constraint τ < c -ω/c is weaker than ( 16) consequently ( 17) holds. Focus on the first term I 1 in (17) and consider:

I 1 = 1 2c M τ (ω/x)-λ-c g (s) ds = 1 2c M τ (ω/x)-λ-c C 1 |s| ds.
Consider the (only) two situations on the sign of (

ω/x) -λ -c = (ω/x) -(1 + ω). If x < ω/ (1 + ω) then (ω/x) -(1 + ω) > 0 and I 1 ≤ C 1 c √ M τ
Notice by the way and for further purpose that:

sup x∈[A τ ,ω/(1+ω)] |x|I 1 ≤ C 1 c ω 1 + ω M τ . If x ≥ ω/ (1 + ω) then (ω/x) -(1 + ω) ≤ 0 and: I 1 = 1 2c M τ (ω/x)-λ-c g (s) ds = 1 2c 0 (ω/x)-(1+ω) C 1 |s| ds + 1 2c M τ 0 C 1 |s| ds = C 1 c |(ω/x) -(1 + ω)| + √ M τ .
Again:

|x|I 1 ≤ C 1 c |x (1 + ω) -ω| + √ xM τ .
From the bounds above we see that:

sup x∈[A τ ,M τ ] |x|I 1 ≤ C 1 c |M τ (1 + ω) -ω| + M τ .
The reader will soon understand why we need to make sure that the right had side in equation above is strictly under C 1 . It is not hard to see that the function τ -→ |M τ (1 + ω) -ω| + M τ is increasing and continuous on [0, c -1]. If we prove that for some δ ∈ ]0, 1[:

1 c |M 0 (1 + ω) -ω| + M 0 = 1 -3δ < 1,
then the existence of some τ + > 0 such that:

1 c |M τ + (1 + ω) -ω| + M τ + = 1 -2δ < 1 (18)
will be granted. But

1 c ω c λ + ω c < 1 2 1 + 1 c < 1.
We turn to I 2 in (17):

I 2 = 1 2c (ω/x)-λ+c M τ g (
C 2 K 2 (ω, c, τ ) ≤ C 1 δ. ( 19
)
Finally comparing ( 15) and ( 19), we see that both conditions cannot be incompatible. Accurate choices of the couple C + 1 , C + 2 are given by the summary bound:

C + 2 ≤ C + 1 min δ K 2 , 1 -δ K 1 . ( 20
)
It is now basic to see that the quadruple C + 1 , C + 2 , τ + , ρ + 1 yields the drift condition (13) for x / ∈ C.

Third step: The remaining step is to check the inequality for some ρ 2 : R g (t) P (t, x) dt ≤ ρ + 1 g (x) + ρ 2 , for any x in C -that is any |x| > M τ (rather x > M τ here as explained above since x > 0). We see that:

0 ≤ ω x ≤ ω M τ ,
and:

1 2c with non-null I 1 and I 2 . Focus on:

I 1 = C 1 2c M τ (ω/x)+λ-c 1 √ s ds = C 1 c √ M τ -(ω/x) + λ -c sup x∈[B τ ,M τ ] √ xI 1 ≤ C 1 c M τ .
At last we see that for Set finally ρ - 1 = (1 + 1/c) /2 < 1. From all that was done above we get: R g (t) P (t, x) dt ≤ C 1 |x|

1 c + C 2 |x| K - 2 ≤ ρ - 1 C 1 |x| ,
whenever:

C 2 K 2 ω, c, τ -≤ 1 -1/c 2 C 1 .

Fig. 1

 1 Fig. 1 Display of constraint A 3 in the plane (ω, c).

Fig. 2

 2 Fig. 2 Display of constraint B 2 in the plane (c, ω).

Fig. 3

 3 Fig. 3 Contour of the Himmelblau's function in the space [-6, 6] 2 .

Fig. 4

 4 Fig. 4 Normal probability plot of H s 1 (N ) on the first coordinate.

Fig. 5

 5 Fig. 5 Normal probability plot of H s 1 (N ) on the second coordinate.

Fig. 6 Fig. 7

 67 Fig.6Trajectory of H s 1 (N ) for an oscillating particle in [-10, 10] 2 . The confidence ellipsoid at a level of 85% is displayed in red. Around 99% of the trajectory of H s 1 (N ) is inside the ellipse.

Fig. 8

 8 Fig. 8 Density of H 1 (2000) with 1150 particles issued from Monte-Carlo simulations. The red ellipse is the 95% confidence ellipsoid.

Fig. 10

 10 Fig. 10 Normal probability plot of H 2 (N ) on the first coordinate.

Figure 11

 11 Figure 11 illustrates different trajectories of H 2 (N ) on the first coordinate which are bounded by the 95% confidence interval deduced from σx .

Fig. 12

 12 Fig. 12 Normal probability plot of H n=200 3

  Fig. 13 Normal probability plot of H n=200 3 (N ) at the 200th iteration on the first coordinate, computed without outlying or belated particles.

1 c

 1 |M 0 (1 + ω) -ω| + M 0 = 1 c ω c λ + ω c .If we assume that λ < ω/c < (1 + c) /4 (assumption B 2 ) then since c > 1:

  The values of the constants C 1 and C 2 were fixed above. Then denote: g (s) ds > 0, then clearly for any x in C: R g (t) P (t, x) dt ≤ ρ + 2 , so that (13) holds.Part 2 (x ≤ 0) We go on with x < 0 and λ > 0, set y = -x ≥ 0, Since g is even and in view of the proposed C we just have to prove exactly the following drift condition with x > 0:1 2c(ω/x)+λ+c (ω/x)+λ-c g (s) ds ≤ ρ 1 g (x) + ρ 2 1 x∈C . First step: Take x / ∈ C. We split [0, M τ ] in two subsets, [0, M τ ] = [0, B τ ] ∪ [B τ , M τ ] with B τ = ω/ (M τ -λ + c) is chosen such that 0 ≤ x ≤ B τ implies the following inequality on the lower bound of the integral:(ω/x)+λ-c > M τ . Clearly B τ ≤ M τ for all τ . Then: ds ≤ C 2 log 2 |(ω/x) + λ + c| .Let: sup0≤x≤B τ |x| log 2 |(ω/x) + λ + c| = K 1 (ω, c, τ ) < +∞. The strictly positive K 1 (ω, c, τ ) exists because |x| log 2 |(ω/x) + λ + c| is bounded on [0, B τ ]. The initial condition reads: ds ≤ ρ 1 C 1 / |x|, 0 ≤ x ≤ B τ whenever C 2 K 1 (ω, c, τ ) ≤ ρ 1 C 1and ρ 1 will be fixed later.Second step: Now we turn to B τ ≤ x ≤ M τ . We still have g (x) = C 1 / |x| but we need to focus on the bounds of the integral.This time the lower bound of the integral (ω/x) + λ -c ∈ [λ -τ, M τ ] and the upper bound(ω/x) + λ + c ∈ [λ + 2c -τ, M τ + 2c]. If we assume that τ ≤ λ then (ω/x) + λ -c ≥ 0 ≥ -M τ . Besides in order that λ + 2c -τ ≥ M τ we just have need that 2c ≥ M τ or τ ≤ c -ω/ (2c) . As a consequence the assumption: τ ≤ min λ, c -ω 2c (21) allows to write for x ∈ [B τ , M τ ], ds ≡ I 1 + I 2 ,

  M τ ≤ 1 i.e. τ ≤ c -ω, sup x∈[B τ ,M τ ] √ xI 1 ≤ C 1 /c. This condition combined with (21) let us set in the sequel: τ ≤ τ -= min (λ, c -ω) .

First, let us focus on r n1 = ε n Ep n +Eg n 2

-Ex n only. Once again Lemma 3 may be invoked for the specific halves of terms 1, 5 and 9 by using the same methods as above. The half part of terms 7 and 10 of (10) containing ε n may be bounded the following way (denote

Lemma 2 ensures that both right hand sides above tend to zero (since sup n z n is almost surely bounded) and Cesaro's Lemma terminates this part.

We should now inspect the terms numbered 1, 5, 7, 9, 10 in (10) with respect to r n2 = e+ε n 2 (ξ n + ν n ). Terms 7 and 9 may be controlled by Lemma 2 and Lemma 3 respectively. Let us deal with term 1:

Take for instance N n=1 ν n . Applying Markov inequality:

then Assumption A 2 together with Cesaro's Lemma again ensure that 1

. Term 5 vanishes in probability with the same technique at hand. Now term 10 gives:

which also tends to zero in probability under A 1-2 . By the way we mention that the cross-product [[r 1n : r 2n ]] in term 10 vanishes due to Lemma 2. Our task is almost done. Let us deal now with the remaining terms of R N . Both:

tends to 0 in probability by Lemma 3 and (6). At last, the remaining [[r n : z n ]] is basically the same as term 1 in (10) and may be addressed the same way. This terminates the proof of the Lemma.

Second case: non-oscillatory and stagnant

We start from 4 and set X n = y n /y n-1 . Then we get :

where ε n = r 1,n + r 2,n -1 has a "witch hat" distribution (convolution of two uniform distributions) with support [-1, +1]. We focus now on the above homogeneous Markov chain X n and we aim at proving that a CLT holds for h (X n ) = log |X n | namely that for some µ and σ 2 :

which will yield:

We aim at applying Theorem 1 p.302 in [START_REF] Jones | On the markov chain central limit theorem[END_REF]. We need to check two points: existence of a small set C and of a function g with a drift condition (see [START_REF] Meyn | Markov chains and stochastic stability[END_REF].

Denote P (t, x) the transition kernel of X n . It is plain that P (t, x) coincides with the density of the uniform distribution on the set

The Theorem 2 is a consequence of the two Lemmas below coupled with the abovementioned Theorem 1 p.302 in [START_REF] Jones | On the markov chain central limit theorem[END_REF].

) is a small set for the transition kernel of X n .

Proof:

We have to show that for all x ∈ C and Borel set A in R:

where ε > 0 and Q is a probability distribution. The main problem here comes from the compact support of P (t, x). Take x such that |x| ≥ M then:

where ε n has compact support [-1, +1]. It is simple to see that with M = M τ = ω/ (c -τ ) the above bound becomes:

The intersection of the supports of 1+ω-2c+τ +cε n and 1+ω-2c+τ +cε n is the set [1 + ω -c -τ, 1 + ω -c + τ ] whatever the value of x in C. The probability measure Q mentioned above may be chosen with support

Now we turn to the drift condition. Our task consists in constructing a function g : R → [1, +∞[ such that for all x: R g (t)

where 0 < ρ 1 < 1 and ρ 2 ≥ 0. Besides, in order to get a CLT on log |X n | we must further ensure that for all x:

Note however that, if (13) holds for g but ( 14) fails, then both conditions will hold for updated function g * = ηg with constant η > 1 and ρ 2 = ηρ 2 such that ( 14) holds.

Lemma 6 Take for g the even function defined by g (x) = C 1 / |x| for |x| ≤ M τ and g (x) = C 2 (log |x|) 2 for |x| > M τ . Assume that:

Then it is always possible to choose three constants τ, C 1 and C 2 such that ( 13) holds for a specific choice of ρ 1 and ρ 2 .

Proof:

The proof of the Lemma just consists in an explicit construction of the abovementioned τ , C 1 , and C 2 . This construction is detailed for the sake of completeness.

At this point and in order to simplify the computations below we will assume that the distribution of ε n is uniform on [-1, +1] instead of the convolution of two U [-1/2,1/2] distributions.

Set λ = 1 + ω -c, assume that λ > 0 (the case λ < 0 follows the same lines) and notice that:

the last inequality stemming from parity of g. We should consider two cases x > 0 and x ≤ 0. The proof takes 2 parts (x > 0 and x < 0 respectively). Both are given again for completeness and because the problem is not symmetric. Each part is split in three steps: the two first steps deal with x / ∈ C, the third with

is chosen such that 0 ≤ x ≤ A τ implies the following inequality on the lower bound of the integral:

Then:

1 2c

This will be combined with the constraint of the first step

. The new condition:

ensures that 1 2c

Third step: The remaining step is to check the inequality:

for any x in C -that is here any x > M τ . Adapting the method given above is straightforward and leads to the desired result with a given ρ - 2 .

We are ready to conclude. Take

Conditions ( 20) and ( 22) hold for the couple (C * 1 , C * 2 ) . For such a couple we have: R g (t) P (t, x) dt ≤ ρ + 1 g (x) + ρ + 2 , x > 0, R g (t) P (t, x) dt ≤ ρ - 1 g (x) + ρ - 2 , x ≤ 0, and for all x: R g (t) P (t, x) dt ≤ max ρ + 1 , ρ -

This finishes the proof of the Lemma.


Dr Michela Zuffo

Dr Xiao Xie

Dr Anton Granzhan
email: anton.granzhan@curie.fr[

Strength in Numbers: Development of a Fluorescence Sensor Array for Secondary Structures of DNA

Keywords: fluorescence sensor array, DNA secondary structures, G-quadruplex, multivariate analysis, highthroughput conformational analysis

High-throughput assessment of secondary structures adopted by DNA oligonucleotides is currently hampered by the lack of suitable biophysical methods. Fluorescent sensors hold great potential in this respect; however, the moderate selectivity that they display for one DNA conformation over the others constitutes a major drawback to the development of accurate assays. Moreover, the use of single sensors impedes a comprehensive classification of the tested sequences. Herein, we propose to overcome these limitations through the development of a fluorescence sensor array constituted by easily accessible, commercial dyes. Multivariate analysis of the emission data matrix produced by the array allows to explore the conformational preferences of DNA sequences of interest, either by calculating the probability of group membership in the six predefined structural categories (three G-quadruplex groups, double-stranded, and two groups of single-stranded forms), or by revealing their particular structural features. The assay enables rapid screening of synthetic DNA oligonucleotides in a 96-well plate format.

Introduction

DNA is highly polymorphic, being capable of forming singlestranded (random coil), double-stranded (A-, B-and Z-DNA), triple-stranded, [1] three-and four-way (Holliday) junctions, [2][3] i-motifs, [4] G-quadruplex (G4) [5] as well as tetrahelical non-G4 structures. [6][7] Among the plethora of secondary structures, G4-DNA are particularly intriguing due to their plasticity and ability to form various topologically diverse isoforms, defined by the number (intra-vs. intermolecular) and mutual orientation of the strands (parallel, hybrid or anti-parallel topologies), as well as the length and topology of the loops. [8][9][10] In view of this conformational diversity, experimental probing of secondary structures is not straightforward, even with short oligonucleotides. High-resolution methods such as NMR or X-ray crystallography can provide detailed structural information, but are time-consuming, require specialized equipment, and often fail if a sample presents a mixture of isoforms. Other biophysical techniques, such as thermal denaturation, [11] thermal difference spectroscopy, [12] hydrodynamic, [START_REF] Garbett | Chaires in G-Quadruplex DNA: Methods and Protocols[END_REF] mass-spectrometric, [START_REF] Marchand | [END_REF][15] chromatographic [16] and gel-electrophoretic methods [START_REF] Moon | Jarstfer in G-Quadruplex DNA: Methods and Protocols[END_REF] may be employed in order to assess the secondary structure of DNA oligonucleotides; however, they provide only limited structural information and are poorly suited for high-throughput analysis. Vibrational circular dichroism (VCD) [START_REF] Andrushchenko | [END_REF][19] and electronic circular dichroism (CD) spectroscopy can furnish a more detailed picture of DNA conformation (such as quantitative analysis of structural elements in G4-DNA), especially when coupled with chemometric analysis. [20][21] Nonetheless, major drawbacks of CD spectroscopy are the requirement of significant amount of material, long acquisition times and lack of high-throughput equipment.

Optical sensors represent an attractive alternative to many biophysical methods and have found countless applications in diverse fields. Among optical methods, fluorescence is most widely utilized, due to simplicity, sensitivity and rapidity of experimental protocols. Moreover, suitable instruments (fluorimeters and plate readers) are available in most laboratories, allowing widespread implementation of high-throughput assays. In recent years, fluorescent sensing has found large application in probing secondary structures of DNA, notably G4-DNA. [22] Thus, numerous fluorescent probes for detection of G4-DNA were discovered or rationally designed by many groups. [START_REF] Largy | Teulade-Fichou in Quadruplex Nucleic Acids[END_REF][START_REF] Vummidi | [END_REF][25] In the best cases, such probes demonstrate selective response towards one or another G4 topology. [26][27][28][29][30][31][32][33] However, most probes do not provide an exclusive selectivity. Moreover, even if they allow a DNA sample to be ascribed to a certain topology group, currently available fluorescent probes do not allow a comprehensive conformational analysis of DNA. Finally, individual probes may provide false outputs with DNA sequences that adopt unusual or non-G4 conformations. [34][35] In this context, application of two dyes for conformational analysis of G4-DNA structures represents a promising step forward. [36][37] Thus, a fluorimetric method that would allow a full conformational analysis of a given DNA sample would be highly desirable.

Sensor arrays (i.e., combinations of several probes) represent useful tools in domains where exquisite specificity is problematic to achieve. Rather than relying on specific interactions with a single analyte, they yield a differential response pattern for each target, enabling its identification. Optical sensor arrays have been developed for numerous analytes and applications (e.g., metal cations, nutrients or contaminants in food, bacteria in patient samples, pollutants in air, soil and water). [38][39] However, most of them only aim at distinguishing individual analytes (e.g., different small molecules of interest), rather than attempting more far-reaching categorizations. Moreover, the applications of this method in the field of nucleic acids sensing are scarce and mainly limited to the recognition of specific sequences, in line with usual applications. [40][41][42] Herein, we describe the first fluorescence sensor array designed for secondary structure analysis of DNA samples. Our work aims at providing a proof-of-concept of the applicability of such an approach to the categorization of oligonucleotides based on their secondary structure instead of their specific identity.

Results and Discussion

For the implementation of an initial sensor array, we selected 11 commercially available dyes with documented affinity for nucleic acids (Figure 1 and Supporting Information, Table S1). The assay is performed in a 96-well plate format, suitable for testing of small amounts of DNA (4-5.5 nmol of oligonucleotide) in an automatized, high-throughput fashion. Following data pre-treatment, DNA conformation is assigned after multivariate analysis of emission dataset for all DNA-dye pairs, recorded at the wavelengths suitable for each dye (Figure 2 and Supporting Information, Figure S1). . Excitation (dashed) and emission (solid lines) spectra of A) ethidium bromide and B) NMM (c = 2.5 µM in K-100 buffer) in the presence of two molar equiv. of G4 (c-myc, blue) or duplex DNA (ds26, red), recorded after 1 h of incubation at room temperature. Emission spectra were also recorded for the dye alone (black lines). Excitation (green) and emission (blue) intervals represent the passbands of the filters used in the plate reader measurements.

We began our study with a set of 34 DNA samples (33 oligonucleotides of 15 to 34-nt length and highly polymerized calf thymus DNA, used at roughly equivalent nucleotide concentration) with well-established conformations (Supporting Information, Table S2). These samples were chosen to uniformly represent the five groups of interest, namely parallel, anti-parallel and hybrid G4s, single-stranded and double-stranded structures.

The CD spectra of all samples were in agreement with literature reports, confirming their expected conformations (Supporting Information, Figure S2). The raw fluorescence intensity data for each dye (c = 2.5 µM) in the presence of all DNA samples (2 molar equiv.) demonstrated an extreme variability of fluorimetric responses (Supporting Information, Figure S3). As such, these results were not suitable for classification of samples through univariate analysis, as demonstrated by Tukey plots obtained by grouping the data according to their expected conformation (Supporting Information, Figure S4). This multi-dimensional dataset (12 emission channels) was thus subjected to intra-channel normalization, followed by Principal Component Analysis (PCA) as described in the experimental section. PCA is a multivariate analysis method that describes the variance in a sample set through a reduced number of principal components (PCs), greatly simplifying the interpretation of multivariate data. [43] Of note, PCA is an unsupervised method that clusters data points based on their similarity, rather than utilizing a predefined classification to orient their disposition: this feature is particularly useful to test the discriminatory power of our sensor array, without considering the conformational premises of samples.

The variance of the full dataset (three replicates for each of the 34 DNA samples) could be described at 95% with six PCs out of 12 (Supporting Information, Table S3), in contrast to most arrays, in which the first two components are sufficient. However, a scree plot of dye eigenvalues (Supporting Information, Figure S5) suggested that only the first four PCs were statistically significant: altogether, they represented 89% of total variance. Analysis of eigenvalues according to the Kaiser rule yielded similar results (Supporting Information, Table S3). A score plot of PC1 vs. PC2 would thus represent alone 65% of sample variance. Visual inspection of the plot allowed us to identify six clusters, highlighted by different colors (Figure 3 A-B). Remarkably, the clustering mirrors, to a large extent, the expected structural classification of samples, with the group of single-stranded samples further dividing into two clusters. These were found to correspond to the sequences containing (ss2, ss3, ss4, ss6, ss7) or not containing (ss1, ss5, dT26) guanosine residues (light blue and orange datasets in Figure 3 A-B, respectively). Depending on the specific dye, such behavior might be ascribed either to the prompt photoinduced oxidation of dG (EdG•+/dG = 1.47 V vs. NHE) [44] or to the conservation of some structured domains, allowing partial intercalation of dyes in the guanine-containing strands. 2F8U) [45][46] and HIV-PRO-1. [47] Notably, some sequences were found to fall relatively far from the centers of the expected conformation clusters (highlighted in Figure 3 A-B). This is the case of Bcl2Mid, HIV-PRO-1 and, to some extent, G4CT. Interestingly, this positioning correlates with the presence of structural peculiarities, which distinguish them from other members of their groups and could impact the response of our array. Specifically, hybrid G4 Bcl2Mid contains an additional A:T base pair, partially hindering one of the external G tetrads. [45][46] Similarly, HIV-PRO-1 contains an additional Watson-Crick G:C base pair of top of a bi-layered anti-parallel G4 structure (Figure 4). [47] Finally, G4CT was shown to adopt two possible conformations, parallel and anti-parallel G4, at high and low potassium concentrations, respectively. [48] Except for these cases, all clusters appeared well separated, considering both the PC1 vs. PC2 (Figure 3, A) and the PC1 vs. PC3 (Figure 3, B) plots. This is particularly remarkable, taking into account the high intragroup polymorphism characterizing the tested sequences.

Next, we attempted to reduce the size of the array, in order to implement a simpler but equally efficacious assay. First, we examined the coefficients of each dye for the first four PCs (Supporting Information, Table S4) and kept only the dyes providing maximal contribution to each of them (PC1: ThT, PC2: Proflavine, PC3: NMM, PC4: Hoechst 33258). However, the separation achieved with such four-dye array did not prove satisfactory (Supporting Information, Figure S6). We thus attempted a less restrictive approach, allowing up to 8 dyes to be part of the array. As evidenced by loading values (grey arrows in Figure 3 A), some pairs of dyes provide similar information (e.g., TO and TO-PRO-3), allowing the exclusion of one of the pair components. We thus retained Berberine, TO-PRO-3 and TMPyP4 (720 nm reading), which gave higher loadings with respect to their companions (i.e., Crystal Violet, TO and TMPyP4-670 nm reading) in the first two, and most significant, PCs. Contextually, we excluded EtBr in favor of Hoechst 33258, which had previously been selected for the four-sensor array. The resulting eight-sensor array proved efficacious and preserved the pattern observed with the full array (Figure 3 C). In fact, only small overlaps were observed in the PC1 vs. PC2 plot between the zones of G-containing single strands, duplexes and anti-parallel G4s. We then attempted to further reduce the array size by eliminating DODCI, according to the same criteria applied for EtBr. However, this resulted in a significant increase of the overlap regions (Figure S7), prompting us to proceed to the following steps with eight sensors. As a next step, we analyzed the dataset originating from this optimized array by Linear Discriminant Analysis (LDA). Similar to PCA, this method defines linear discriminant functions that best separate data points from a set. However, LDA is a supervised method, which utilizes the information on each analyte's group membership to optimize the separation of groups. [43] The practical interest of such analysis resides in the possibility to assign an unknown analyte to one of the predefined groups. We thus defined the relevant functions (up to 5 canonical variables) using the 34 DNA samples as a training set. Validation of the method by the leave-one-out test proved satisfactory, with only one misclassified DNA sample out of 34 (2.9% error, Supporting Information, Table S5), confirming the potential of the sensor array approach even when applied to a task as complex as ours. The classification is determined upon calculation of the probability for each sample's membership in each of the six predefined groups, taking into account the whole set of canonical variables. A visual representation of the outcome is provided by the 2D LDA plot (Canonical Variables 1 vs. 2, Figure 5). As previously observed in PCA, the six groups form as many independent clusters, with minimal overlaps. A similar result is observed in the 3D plot obtained upon addition of the third canonical variable (Supporting Information, Figure S8).

Further, we tested several practical applications of our assay. First, we examined multiple DNA sequences whose CD spectra were problematic to interpret. In some cases, CD spectra suggest a simplified structure with respect to the conformation established by high-resolution methods. For example, CD signatures of the sequences Pu24T and c-kit87up (cf. Supporting Information, Table S2) suggest a parallel G4 conformation (Supporting Information, Figure S9). However, their solution structures are much more complex, as demonstrated by NMR analysis (PDB: 2A5P and 2O3M, respectively). [49][50] Interestingly, PCA analysis highlighted these peculiarities, positioning them in between the parallel and hybrid groups (Supporting Information, Figure S10). This was also confirmed by visual inspection of the LDA plot and the derived classification, which ascribed both Pu24T and c-kit87up to the hybrid G4 group with high probability (Figure 5 and Table 1). Likewise, solution structure of 22GGG as established by NMR (PDB: 2KF8) features a bi-tetrad, antiparallel (basket-type) conformation, [51] whereas its CD spectrum suggests a hybrid conformation (Supporting Information, Figure S2). Instead, PCA analysis (Supporting Information, Figure S11 A) indicated the similarity of 22GGG to other members of the antiparallel G4 group, as further confirmed by LDA assignment (Supporting Information, Figure S11 B and Table S6). In some other cases, CD spectra allow dual interpretation, as in the case of UpsB-Q3 [52] and EBR1 [53] (cf. Supporting Information, Table S2). CD spectra of both sequences display a negative maximum at 240 nm and a positive maximum at 263 nm, with a shoulder between 283 and 310 nm (Figure S12). The presence of these two positive bands might arise either from a single hybrid conformation, or from the co-existence of two distinct structures. PCA (Supporting Information, Figure S13) and LDA assessments (Figure 5 and Table 1) prompted us to confirm the assignment of UpsB-Q3 to the hybrid G4 group. For EBR1, instead, the data points are located in the middle of the PC1 vs. PC2 score plot, that is, between the parallel and anti-parallel G4 groups. In good agreement, EBR1 was not classified as hybrid by LDA (Figure 5 and Table 1). Therefore, these results support the interpretation of its CD spectrum in favor of a mixture of secondary structures.

Next, we examined the proficiency of our sensor array at dealing with secondary structures different from the ones used in the training set. This is the case of several sequences displaying CD signatures non-typical for single strands (Supporting Information, Figure S14), despite the lack of obvious hydrogen bonding patterns allowing the formation of double-stranded or G4 structures, namely scr26, ss8, 24non096, dA26 and (TGA)8T, as well as VK2 forming a peculiar non-G4 tetraplex structure (cf. Supporting Information, Table S2). [7] The results of PCA analysis revealed a positioning of data points that was ambiguous, but clearly falling outside the expected single-strand groups (Supporting Information, Figure S15 A). A similar result was observed using LDA (Supporting Information, Figure S15 B). Although our approach does not allow to precisely identify the conformation of such oligonucleotides, it still provides valuable information on the presence of unexpected secondary structures.

Finally, we evaluated the capability of our assay to analyze mixtures of analytes. For this purpose, we tested samples containing varying molar fractions of ds26 (self-complementary duplex) and 24TTA (hybrid G4), as well as c-myc and Bm-U16 (parallel and anti-parallel G4s, respectively). In both cases, CD spectra of mixtures represented arithmetical sum of the ones of individual components, resulting in complex shapes unsuitable for simple qualitative analysis (Supporting Information, Figure S16). A PC1 vs. PC2 score plot for the ds26 / 24TTA mixtures revealed a progressive shift from the group of duplexes to the one of hybrid G4s (Figure 6 A), with only pure samples clearly falling into the corresponding clusters. Interpretation of the data for mixtures c-myc / Bm-U16 was more complex, since all samples containing more than 20 mol% of c-myc spotted in the parallel G4 group zone (Figure 6 B). This is likely due to higher affinity of some of the sensors (e.g., NMM) for the parallel G4 structure. In all cases, though, data sets for all mixtures were clearly distinguishable from pure samples.

Conclusions

In this work, we presented a novel fluorescence-based method, based on a sensor array composed of easily available dyes (which allows an easy implementation by other users) and allowing the assessment of secondary structures adopted by DNA oligonucleotides in a high-throughput fashion. Our first-generation array is able to satisfactorily discriminate among the three topological classes of G4-DNA, single-stranded and duplex oligonucleotides, correctly assigning 97.1% of tested sequences through the use of multivariate analysis tools (PCA and LDA), available in commonly used and/or freeware analysis software (such as Origin and R). Moreover, we demonstrated that the cases appearing rather far from the barycenter of expected groups (HIV-PRO-1, Bcl2Mid, G4CT) can be interpreted as an indication of particular structural features of a given sample, distinguishing them from typical members of the expected group. Thus, this method may be exploited not only for categorization of samples but also for the discovery of particular motifs.

We explored some practical applications of this assay and demonstrated that it provides useful insights on a variety of sequences whose conformation cannot be assessed by simple interpretation of CD spectra. The assay also proved efficacious at distinguishing mixtures of analytes, folded into different conformations, with respect to pure samples. Of note, the assay design allows plenty of room for improvement through, e.g., optimized selection of dyes and wavelength settings, and may be easily adapted to other DNA conformations, such as i-motif, triplehelical or junction structures, through inclusion of corresponding sequences into the training set. It may also constitute the basis for a high-throughput method to experimentally assess the G4folding potential and the preferentially adopted topology of novel DNA sequences. Of course, this approach, similar to most in vitro biophysical techniques (NMR, ESI-MS, CD spectroscopy) is limited to relatively short synthetic oligonucleotides. However, the strong advantage of sensor arrays, along with its versatility, is the potential of high-throughput analysis (up to thousands of sequences), which is currently unreachable with other techniques.

Experimental Section

Buffer and dye solutions: All experiments were performed using a 'K-100' buffer containing 0. ) was diluted with K-100 buffer to c ≈3 mM in base pairs, as calculated by absorption measurement at 260 nm based on εnucleotide = 6650 cm -1 M -1 , [54] and then further diluted so as to obtain a working solution with a comparable nucleotide concentration as in oligonucleotide samples (i.e., 123 µM, considering 22 as an average length of oligonucleotides, cf. Supporting Information, Table S2). Working solutions were subsequently annealed (5 min at 95 °C) and then let equilibrate overnight. Annealed solutions were stored at 4 °C. For DNA mixtures, the components were annealed separately at a concentration of 5.6 µM, according to the same procedure and, after equilibration, mixed in relevant quantities (DNA1 : DNA2 = 100:0, 80:20, 60:40, 40:60, 20:80, 0:100, v/v).

CD spectra: CD spectra were recorded with a Jasco J-1500 spectropolarimeter. Spectra were recorded using working solutions of the pure DNA samples or of their relevant mixtures in K-100 buffer, in quartz cuvettes with a path length of 0.5 cm. For all samples, total strand concentration was 5.6 µM, except for 46AG (2.8 µM) and ct DNA (123 µM nucleotides). In the case of DNA mixtures, the total oligonucleotide concentration was also 5.6 µM. Parameters for spectra acquisition: wavelength range, 210 to 330 nm; scan speed, 50 nm min -1 ; number of scans, 3; data pitch, 0.5 nm; bandwidth, 2 nm; integration time, 1 s; temperature, 22 °C. Spectra were subsequently corrected for the blank. Finally, spectra were converted to molar dichroic absorption Δε [M -1 cm -1 ] = θ/(32980 × c × ℓ), where θ is the CD ellipticity in millidegrees (mdeg), c is DNA concentration in M, and ℓ is the pathlength in cm.

Fluorescence excitation and emission spectra: Fluorescence excitation and emission spectra were recorded with a HORIBA Jobin-Yvon Fluoromax-3 spectrofluorimeter, in quartz cuvettes (path length of 1 × 0.5 cm for emission and excitation beams, respectively). Sample solutions contained 2.5 µM of dye and 5 µM of the relevant DNA (G4 c-myc or duplex ds26), unless otherwise stated, in K-100 buffer. After mixing, the solutions were left to equilibrate for 1 h in the dark, before performing the measurements. Emission and excitation wavelengths are listed in Supporting Information, Table S1 and correspond to those used in the following high-throughput assay.

Assay experimental procedure: Working solutions of DNA samples (5.6 µM or equivalent in K-100 buffer, 1.1 mL) or buffer alone were filled into 1.4-mL round-bottom tubes (Micronic). Working solutions of the dyes prepared as described above (25 µM, 5 mL) or buffer alone were filled into a 12-column reservoir (Porvoir Sciences). The solutions were transferred into 96-well black, flat-bottom microplates (Corning model 3650) using an automated pipetting system (Gilson PIPETMAX-268). Specifically, 90 µL of DNA solution (or buffer) were transferred into each of the 12 columns of the plate. Afterwards, 10 µL of dye solution (or buffer) were transferred into the appropriate columns. Final concentrations were 2.5 µM of the dye and 5 µM of DNA (or equivalent, cf. Supporting Information, Table S2), in a total volume of 100 µL per well. The final solvent mixture contained 1% v/v of DMSO. The plates were stirred for 1 min at 400 rpm and then let equilibrate in the dark for 1 hour at ambient temperature. Fluorescence emission was recorded using a microplate reader (BMG FluoStar Omega), exciting each dye at the appropriate wavelengths, with the aid of appropriate filters (see Supporting Information, Table S1 for excitation and emission wavelengths, and Figure S1 for filter passbands). The instrument gain was set for each channel and kept constant throughout all analyses. All samples were analyzed by three independent preparations and measurements.

Multivariate analysis:

The acquired emission intensity data were normalized to a 0 to 1 scale separately for each dye. The emission intensity of the dye alone was included in the pre-treatment and then excluded from multivariate analysis, since it does not belong to any conformational group. Normalization is recommended for PCA to avoid biases due to different orders of magnitude for the various variables. Other pre-treatment methods (sample/free dye fluorescence intensity ratio, [55] autoscaling [56] ) were also tested, but did not produce significantly different results. Multivariate analysis (PCA and LDA) was performed using Origin Pro 2018b (OriginLab, Northampton, MA). Data for PCA are presented as score plots of PC1 vs. PC2, unless otherwise stated. Data for LDA are presented as Canonical Variables 1 vs. 2 plots, unless otherwise stated. Leave-one-out test was used for internal validation of the LDA method.

FULL PAPER

Entry for the

Figure 1 .

 1 Figure 1. Principle of the proposed assay. A) Selection of fluorescent dyes; B) preparation of dye-DNA sample combinations in a 96-well plate format and fluorescence emission measurement; C) data pre-treatment (intra-channel normalization); D) multivariate analysis.

Figure 2

 2 Figure2. Excitation (dashed) and emission (solid lines) spectra of A) ethidium bromide and B) NMM (c = 2.5 µM in K-100 buffer) in the presence of two molar equiv. of G4 (c-myc, blue) or duplex DNA (ds26, red), recorded after 1 h of incubation at room temperature. Emission spectra were also recorded for the dye alone (black lines). Excitation (green) and emission (blue) intervals represent the passbands of the filters used in the plate reader measurements.

Figure 3 .

 3 Figure 3. A) PC1 vs. PC2 and B) PC1 vs. PC3 score plots derived from PCA analysis of the training set (34 DNA samples × 3 replicates, plotted separately as colored dots) using the data from the full sensor array (11 dyes, 12 variables). C) PC1 vs. PC2 score plot representing the PCA analysis of the data from the reduced sensor array (8 dyes). Grey arrows: loading vectors for each variable. Ellipses represent 90% confidence intervals for each group.

Figure 4 .

 4 Figure 4. Schematic representation of G4 structures adopted by Bcl2Mid (PDB: 2F8U)[45][46] and HIV-PRO-1.[47]

Figure 6 .

 6 Figure 6. PC1 vs. PC2 score plot derived from PCA analysis of the training set (empty dots) and A) ds26 / 24TTA mixtures or B) Bm-U16 / c-myc mixtures (filled dots).

1 :A 1 : 2 : 6 : 4 : 3 : 40 B

 11264340 ds26 / 24TTA 100:0 2: ds26 / 24TTA 80:20 3: ds26 / 24TTA 60:40 4: ds26 / 24TTA 40:60 5: ds26 / 24TTA 20:80 6: ds26 / 24TTA 0:100 Bm-U16 / c-myc 100:0 Bm-U16 / c-myc 80:20 Bm-U16 / c-myc 0:100 Bm-U16 / c-myc 40:60 5: Bm-U16 / c-myc 20:80 Bm-U16 / c-myc 60:

Table 1 .

 1 Assignment of Pu24T, c-kit87up, UpsB-Q3 and EBR1 samples to the relevant conformation groups by LDA.

					Assignment probability [a]
	Sequence							Predicted group
		Hybrid G4	Anti-parallel G4	Parallel G4	Duplex	Single strands (with dG)	Single strands (without dG)
	Pu24T	0.99	3.88 × 10 -6	4.15 × 10 -11	8.67 × 10 -15	4.18 × 10 -6	2.42 × 10 -42	Hybrid
	c-kit87up	0.98	1.66 × 10 -2	3.02 × 10 -11	8.83 × 10 -13	1.04 × 10 -4	1.50 × 10 -27	Hybrid
	UpsB-Q3	1	2.25 × 10 -5	1.88 × 10 -28	1.27 × 10 -38	4.34 × 10 -19	2.62 × 10 -29	Hybrid
	EBR1	3.77 x 10 -6	0.99	6.90 × 10 -18	8.78 × 10 -20	6.18 × 10 -5	1.08 × 10 -11	Anti-parallel
	[a] Analysis performed using average values calculated from three independent replicates for each sample.

Figure 5. LDA analysis of the training set (34 DNA samples × 3 replicates, plotted separately as colored dots). Ellipses represent 85% confidence intervals for each group. Samples colored in magenta represent the test samples (Pu24T, c-kit87up, EBR1 and UpsB-Q3, 3 replicates per sample), analyzed according to the established canonical variables.

  and ds 6 + 7) were prepared by mixing equal volumes of the corresponding single strands. The final concentration of heteroduplexes was 2.8 µM, to account for the doubled number of nucleotides. Calf thymus DNA (ct DNA, Invitrogen, 10 mg mL-1

	1 M KCl and 0.01 M lithium cacodylate in MilliQ
	water, pH 7.2. Dyes were purchased from Sigma-Aldrich (Thiazole Orange,
	Crystal Violet, DODCI, TMPyP4), TCI (Berberine hemisulfate, Thioflavin
	T), Fluka (Hoechst 33258, Proflavine, Ethidium Bromide), Invitrogen (TO-
	PRO-3) or Frontier Scientific (NMM), and used as supplied without further
	purification. Stock solutions of dyes were prepared in DMSO at a
	concentration of 4 mM and stored at -20 °C. Working solutions (25 µM)
	were prepared by diluting the stock solutions with a mixture of milliQ water
	and DMSO (9:1, v/v) containing 0.1 M KCl and 0.01 M lithium cacodylate
	(pH 7.2), and stored at 4 °C. The absence of precipitation in such
	conditions was verified for all dyes by centrifugation and visual inspection
	of the resulting solutions.
	DNA samples: Oligonucleotides were purchased from Eurogentec (RP-
	HPLC purification grade) and used without further purification. Stock
	solutions of oligonucleotides (strand concentration: 100 µM, except for
	46AG: 50 µM) were prepared in K-100 buffer and stored at 4 °C. Working
	solutions (c = 5.6 µM, except for 46AG: c = 2.8 µM) were prepared by
	dilution of stock solutions in the same buffer. Heteroduplexes (ds 3 + 4, ds
	5 + 8

Table of Contents

 of FULL PAPER Text for Table of ContentsMichela Zuffo, Xiao Xie, Anton Granzhan*

Page No. -Page No. Strength in Numbers: Development of a Fluorescence Sensor Array for Secondary Structures of DNA

Acknowledgements

The authors thank Dr. M.-P. Teulade-Fichou for helpful and stimulating discussions. Financial support by Institut Curie (postdoctoral fellowship to M.Z.) and National Research Agency (ANR-17-CE07-0004-01, to A.G.) is gratefully acknowledged.


Estelle Couallier
email: estelle.couallier@univ-nantes.fr

Shuli Liu

Erika Clavijo

Liliana Villafaña Lopéz

Matthieu Frappart

Recovery of lipids from microalgae extracts by membrane processes: Comparison of cross-flow and shear-enhanced filtration performances

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

Recovery of lipids from microalgae extracts by membrane processes: Comparison of cross-flow and shear-enhanced filtration performances

Introduction

The biorefinery of renewable resources like microalgae offers great opportunities to substitute biomolecules to traditional raw materials in various industry sectors. Such strategies necessitate innovative choices of soft and energy-efficient processes to guarantee the integrity of fragile molecules and develop eco-friendly production. For large-scale production (food, energy or green chemistry), a wet processing of biomass has been proposed, that avoids expensive drying steps and reduces solvent use. However the energetically efficient extraction of biomolecules at low cost and industrial scale is not yet mature. Biomass wet treatment includes 1-the harvesting, 2the cell disruption step to release the valuable biochemical compounds in the aqueous phase, 3the fractionation step (extraction, concentration and purification). The integration of membrane processes into the microalgae downstream processing concerns the harvesting and the concentration of microalgae, but membrane filtration is also a promising clean separation process for the fractionation step.

In this work we focused on the recovery of lipids from Parachlorella kessleri, cultivated in starving conditions to enhance their lipid production. Clavijo et al. [1] demonstrated that the supernatant after bead milling and centrifugation of Parachlorella kessleri contains an important part of emulsified lipids in aqueous phase (until 50 % of total lipids).According to literature, the lipid recovery from microalgae is mostly performed with supercritical CO 2 on dried matter or solvent extraction. In this study the recovery of lipids from aqueous extracts was performed by membrane processes.

Methods

A model solution was formulated, based on the analysis of ground Parachlorella kessleri [1]. Then the performances (retention, flux, fouling, cleanability) of polyacrylonitrile (PAN), polyethersulfone (PES) and polyvinylidene fluoride (PVDF) membranes to concentrate lipids from Parachlorella kessleri aqueous extracts were evaluated, using the model solution. The most appropriate material and conditions (TMP) were then selected and validated on real microalgae fractions. Two filtration modes were compared allowing : the crossflow filtration in a plane module (CF) and the shear-enhanced filtration in a rotating disk module (RD).

Results and discussion

The PAN 500kDa membrane presented the best performances (flux, lipid retention and cleanability) in crossflow and rotating disk filtrations. It was selected to test the filtration of real aqueous extracts from Parachlorella kessleri. The performances of the PAN membrane in CF module with the real products were similar to the ones with the model solution despite the composition differences. The lipids were totally retained whereas some of the hydrophilic compounds (polysaccharides and proteins depending on the conditions) could permeate. The RD filtration allowed reducing the membrane fouling and using a higher transmembrane pressure, thus it led to higher fluxes. In the case of real aqueous extracts RD filtration, the water permeability was enhanced by the accumulation of hydrophilic compounds whereas the shear rate limited the membrane clogging by the lipids. The RD filtration using a PAN membrane presented the best performances.

Conclusions

The membrane filtration is a relevant separation process to concentrate the lipids, reduce the volume of water to treat for further purification steps and it offers interesting opportunities for the fractionation of biomolecules from microalgae. An optimization of the coupling with the upstream operations (culture, cells disruption, clarification) would enhance the separation.

Figure 1 .

 1 Figure 1. Crossflow and rotating-disk filtrations of an emulsion (model solution of concentrated microalgae aqueous extract) and two real aqueous extracts (SUP1 and SUP2) with a PAN membrane.