Eric Carlen
email: carlen@math.rutgers.edu

Pierre Degond
email: pierre.degond@math.univ-toulouse.fr

Bernt Wennberg
email: wennberg@chalmers.se

Kinetic limits for pair-interaction driven master equations and biological swarm models

Keywords: Master equation, kinetic equations, binary interactions, propagation of chaos, Kac's master equations, swarms, correlation AMS Subject classification: 35Q20, 35Q70, 35Q82, 35Q92, 60J75, 60K35, 82C21, 82C22, 82C31, 92D50

We consider a class of stochastic processes modeling binary interactions in an N -particle system. Examples of such systems can be found in the modeling of biological swarms. They lead to the definition of a class of master equations that we call pair interaction driven master equations. We prove a propagation of chaos result for this class of master equations which generalizes Mark Kac's well know result for the Kac model in kinetic theory. We use this result to study kinetic limits for two biological swarm models. We show that propagation of chaos may be lost at large times and we exhibit an example where the invariant density is not chaotic. hospitality during their visits there.

alignment tendency: at each collision, the two particles change their velocity to their average velocity up to a certain noise. One of the goals of the present paper is to provide a rigorous justification of this kinetic model, at least in the space-homogeneous case.

Here, we also propose a different and to our knowledge original binary collision process which mimics the Vicsek alignment dynamics. In this process called 'Choose the Leader (CL)', one of the two colliding particles (the follower) decides to take the velocity of the other one (the leader) up to some noise. The choice of the leader and the follower is random with equal probabilities. We propose a kinetic formulation of the CL process and rigorously establish it in the space homogeneous case. One of the advantages of the CL process, from the mathematical viewpoint, is that it leads to a closed hierarchy of marginal equations (or BBGKY hierarchy). We will make use of this opportunity to provide explicit computations of the correlations, i.e. of the distance to statistical independence.

The BGD and CL processes are special examples of a general class of pair-interaction processes. The paper will study these processes in the space-homogeneous case. Having in mind the special examples of the BDG and CL processes, we assume that the particle velocities are two-dimensional vectors of constant norm. However, this assumption could be easily waived. The main theorem is that the chaos propagation property is true for these pair-interaction processes. The derivation of the BDG and CL kinetic equations follow from this theorem. In the case of the BDG operator, we recover the collision operator or [START_REF] Bertin | Boltzmann and hydrodynamic description for selfpropelled particles[END_REF]. The proof of the theorem generalizes some of the combinatorial arguments of Mark Kac [START_REF] Kac | Foundations of kinetic theory[END_REF].

This result can be seen as paradoxical at first sight. Indeed, the BDG and CL processes build-up correlations, in the sense that particles tend to eventually become close to eachother (in velocity space). This correlation build-up must be present in any model that is to display "swarming" or "flocking" behavior, and it should lead to a breakdown of the statistical independence of the particles. The resolution of this apparent paradox is seen through an inspection of the time scales. Indeed, the theorem is only valid on any finite time interval at the kinetic scale. At this scale, the average number of collisions undergone by each particle on a finite time interval is finite and uniformly bounded. If N is large and one selects a typical pair of particles and a finite time t, none of the particles that have collided with the first particles by that time have also collided with any particle that has collided with the other particle at time t. If one fixes N, and lets t become large, pair correlations do develop. However, if one fixes t and lets N become large, the amount of pair correlation built up by time t goes to zero as N increases.

At equilibrium, the solution of the master equation converges to the so-called invariant density. Due to the special properties of the CL process, it is possible to provide closed expressions for the marginals of the invariant density. When the noise distribution is appropriately scaled as N → ∞, we show that the invariant density cannot be chaotic. Indeed, while the single particle marginal density is uniform, the two-particle marginal density is not. Therefore, the two-particle marginal density is not a tensor product of two copies of the single-particle marginal density, as it should if it would be chaotic. Numerical simulations, reported in a forthcoming work [START_REF] Carlen | Kinetic hierarchy and chaos propagation in biological swarm models[END_REF], confirm this result. Again, it may seem paradoxical that the invariant measure has strong pair correlations and yet chaos is propagated, and again, it is consideration of the time scales that resolves the paradox: The time for the N particle system to approach equilibrium grows so rapidly with N that the stationary measure in not relevant to time evolution on the kinetic scale.

The message of the paper is that the chaos property may be true even for processes that seemingly build-up correlations. However, in this case, the correlation build-up capacities of the processes under considerations only manifest themselves at scales which are large compared to the kinetic scale. To describe these systems at these large scales, kinetic theory is not valid anymore, and alternate theories must be devised. So far, the subject is widely open in the literature and constitute a fascinating area of research.

We conclude this section by a few more bibliographical remarks. In [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], an alternate kinetic model for the Vicsek system has been proposed. It consists of a nonlinear Fokker-Planck equation. It has been derived from a mean-field limit of the Vicsek system in [START_REF] Bolley | Mean-field limit for the stochastic Vicsek model[END_REF]. In biological swarm modeling, most of the authors make use of particle (aka 'Individual-Based') models (see e.g. [START_REF] Aoki | A simulation study on the schooling mechanism in fish[END_REF][START_REF] Chuang | State transitions and the continuum limit for a 2D interacting, self-propelled particle system[END_REF][START_REF] Couzin | Collective Memory and Spatial Sorting in Animal Groups[END_REF][START_REF] Cucker | Emergent Behavior in Flocks[END_REF][START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF]) and sometimes, fluid-like hydrodynamic models (see e.g. [START_REF] Chuang | State transitions and the continuum limit for a 2D interacting, self-propelled particle system[END_REF][START_REF] Mogilner | A non-local model for a swarm[END_REF][START_REF] Topaz | A nonlocal continuum model for biological aggregation[END_REF]). The use of kinetic models is more rare. There is a kinetic version of the Cucker-Smale model [START_REF] Carrillo | Asymptotic flocking dynamics for the kinetic Cucker-Smale model[END_REF]. This kinetic Cucker-Smale model takes the form of a nonlinear non-local Fokker-Planck equation which has been rigorously derived from the mean-field limit of the discrete Cucker-Smale model in [START_REF] Bolley | Stochastic mean-field limit: non-Lipschitz forces & swarming[END_REF]. Kinetic models have also been proposed in the context of fish schools [START_REF] Degond | A macroscopic model for a system of swarming agents using curvature control[END_REF], bacteria and cell motion [START_REF] Hillen | The diffusion limit of transport equations derived from velocity-jump processes[END_REF][START_REF] Painter | Modelling cell migration strategies in the extracellular matrix[END_REF] and ant-trail formation [START_REF] Boissard | Trail formation based on directed pheromone deposition[END_REF]. In most cases, their justification is purely formal. The present paper is the first step towards a justification of kinetic models in biological swarm modeling.

The paper is organized as follows. In section 2, we derive the master equations for the two examples of biological swarm models that we will consider, the BDG and the CL processes. We also provide the definition of the most general pair interaction driven master equation. Section 3 is devoted to the proof of the chaos propagation theorem for pair interaction driven master equations and its application to the BDG and CL dynamics. This proof is modeled in part on the orginal approach of Kac, but with some differences. Section 4 investigates the invariant measure for the CL process and shows that, under some suitable scaling of the noise distribution, it violates the chaos property. A conclusion is drawn in section 5.

2 Motivation: biological swarm models

Two examples

We consider a population of N agents, which to be concrete, we take to be fish in a shallow pond, each swimming at unit speed. In this work, we are only concerned with the evolution of the velocities, and we neglect the precise spatial location of the fish; that is, we assume that all N under consideration are sufficiently close to interact with one another.

The shallow pond is essentially a planar domain, and so the individual velocity vectors belong to the unit circle S 1 . The state space of the system is therefore the torus

T N = [S 1 ] ×N .
The state of the swarm, or school, is then specified by giving a vector

v = (v 1 , . . . , v N ) ∈ T N .
All N agents, or fish, are considered as part of a local cluster, and all are interacting with one another. The evolution of v = (v 1 , . . . , v N ) will be modeled in various ways, all based on the following general scheme: there is a steady Poisson stream of jump times, at which a pair (i, j) is selected at random from {1, . . . , N}, and then these two fish adjust their velocities in some way

(v i , v j ) -→ (v ′ i , v ′ j ) . (2.1)
To complete the specification of the dynamics, we need to give the precise rule for updating the velocities in (2.1). Here are the rules we shall consider.

(1) BDG dynamics: This rule is designed to lead to a kinetic model first investigated by Bertin, Droz and Gregoire [START_REF] Bertin | Boltzmann and hydrodynamic description for selfpropelled particles[END_REF]. The idea is that the pair of agents adjust their velocities cooperatively to achieve the same direction of motion, apart from some noise in their adjustments. More precisely, define

v i,j = v i + v j |v i + v j | . (2.2)
Now think of v i,j ∈ S 1 as unit complex number. Let W i and W j be two more unit complex numbers, chosen independently at random from a probability distribution g(w)dw on S 1 , and define

v ′ i = W i v i,j and v ′ j = W j v i,j . (2.3)
Regarding S 1 as the unit circle in the complex plane, W j v i,j simply means the product in the complex plane of the random variables W i and v i,j , and likewise for the other term. Thus, if we write

W i = e iΘ and v i,j = e iθ , then W i v i,j = e i(Θ+θ)
so that the noise is additive in the angles. In the case of no noise, W i and W j are the constant random variable W i = W j = 1, and in the case of small noise, they are random variables that are strongly peaked around 1. We suppose that g(w) is symmetric; i.e., g(w) = g(w * ), (2.4) where w * denote the complex conjugate of w, and that g(w) is somewhat peaked near w = 1.

Remark 2.1 This rule could be referred to as the "Maxwellian BDG" dynamics, in reference to the fact that the selection of the pair (i, j) is independent of their relative velocity, like the Maxwellian molecular interaction in rarefied gas dynamics. A more general setting would make the collision probability of the pair (i, j) depend on their relative velocity v i v * j , but this will be discarded here, for reasons developed below (see remark 3.2).

(2) "Choose the leader" (CL) dynamics: In this variant, one of the two agents in the pair decides to adopt the other agent's velocity, though it does not get this velocity exactly right: The new velocity it adopts is the velocity of the other agent up to some noise term. The process is written as follows: if the pair selected is (1, 2), and agent 1 decides to adopt the velocity of agent 2, then the velocity of agent 1 is updated as follows:

v 1 → v ′ 1 := Zv 2 ,
where Z is an independent random variable with values in S 1 and probability g. As before, we regard S 1 as the unit circle in the complex plane, and Zv 1 simply means the product in the complex plane of the random variables Z and v 1 . Here again, we assume that g is symmetric and satisfies (2.4) for the sake of simplicity. We also have in mind that g is peaked around 1.

We use a fair coin toss, modeled by a Bernoulli variable B, to decide which agent adopts the velocity of the other, and then we have the following description of the jump when pair (i, j) The selected pair of velocities is updated according to

v ′ i = Bv i + (1 -B)Zv j v ′ j = BZv i + (1 -B)v j (2.5)
and all other velocities are unchanged.

There are many other variants on these basic examples, but for now, let us focus on these two and seek a passage from this description of the interactions of individual agents to an evolution equation for the statistical distribution of the velocities in the system. For this we shall employ methods of kinetic theory that have been developed for a similar problem concerning colliding molecules in a gas. We shall use a probabilistic approach of Marc Kac [START_REF] Kac | Foundations of kinetic theory[END_REF], using a so-called Master equation.

Master Equations

General framework

We now derive master equations describing the evolution of the probability density for the state of the system as it undergoes our stochastic processes. An advantage with starting at the microscopic level, i.e., the level of individual agents, is that the modeling is much clearer before any large N limits are taken.

We note that both the BDG or the CL dynamics are Markovian. In general, we consider a Markov process on T N and we denote by V k ∈ T N its state just after the kth jump. We define its Markov transition operator Q as usual by

Qϕ( v) = E{ϕ( V k+1 ) | V k = v } ,
for any continuous test function ϕ on T N . Now let F k ( v) denote the probability density of V k (with respect to the uniform measure on T N ). Then by definition, one has

E(ϕ( V k+1 )) = T N ϕ( v)F k+1 ( v)d N v .
On the other hand, by standard properties of the conditional expectation,

E(ϕ( V k+1 )) = E(E{ϕ( V k+1 ) | V k }) = T N Qϕ( v)F k ( v)d N v .
That is,

F k+1 = Q * F k , (2.1)
where

Q * is the adjoint of Q in L 2 (T N , d N v).
The next step is to construct a time continuous process which will lead to a timecontinuous master equation. The state of the process is now a function of time v(t) ∈ T N and the probability density F ( v, t) is a function of the continuous time parameter t instead of the discrete jump index k. We assume that F ( v, t + dt) only depends on F ( v, t) and not on the past values F ( v, s) for s < t. In this way, we can construct a time-continuous Markov process. Thus, we assume that between time t and t + dt the probability that a collision occurs is

λQ * F ( v, t) dt + o(dt),
where λ is a constant. We assume that the probability that multiple collisions occur in the time interval [t, t + dt] is negligible. Thus, since there are N particles colliding independently, we have, for dt:

F ( v, t + dt) = NλQ * F ( v, t) dt + (1 -NλF ( v, t)) dt + o(dt).
Equivalently, we can write

F ( v, t + dt) -F ( v, t) = Nλ[Q * F ( v, t) dt -F ( v, t)] dt + o(dt),
which leads to the Master Equation in the limit dt → 0:

d dt F ( v, t) = λN[Q * -I]F ( v, t),
where I is the identity operator. This equation must be complemented by the initial condition

F ( v, 0) = F 0 ( v)
where F 0 is the initial probability distribution. In the remainder, we scale time in such a way that λ = 1 and we define

L = N[Q -I] and L * = N[Q * -I]. (2.2)
We summarize the previous discussion in the following definition:

Definition 2.1
The time-continuous master equation associated to a discrete Markov process of transition operator Q is written

d dt F ( v, t) = L * F ( v, t), (2.3)
F ( v, 0) = F 0 ( v). (2.4)
where

L * = N[Q * -I] with Q * the adjoint of Q and F 0 is the initial probability distribution.
In the next sections, we determine the master equations of the BDG and CL processes successively.

The BDG dynamics

We state the following:

Proposition 2.

The master equation of the BDG dynamics is written (2.3) with

L * = N N 2 -1 i<j (Q * (i,j) -I) , (2.5)
and the binary interaction operator Q * (i,j) given by:

Q * (i,j) F ( v) = T 2 F (v 1 , . . . , y i , . . . , y j , . . . , v N ) g v i y * i,j g v j y * i,j dy i dy j . (2.6)
Proof: ¿From the definition (2.3) of the BDG dynamics, we get

Qϕ( v) = E{ϕ( V k+1 ) | V k = v } = Eϕ(v 1 , . . . , W i v i,j , . . . , W j v i,j , . . . , v N ) = 2 N(N -1) i<j T 2 ϕ(v 1 , . . . , w i v i,j , . . . , w j v i,j , . . . , v N ) g(w i ) g(w j ) dw i dw j .
To compute the adjoint of Q, we note that for any probability density F on T N ,

T N F (v 1 , . . . , v N ) × T 2 ϕ(v 1 , . . . , w i v i,j , . . . , w j v i,j , . . . , v N ) g(w i ) g(w j ) dw i dw j dv 1 . . . dv N = = T N T 2 F (v 1 , . . . , v i , . . . , v j , . . . , v N ) g y i v * i,j g y j v * i,j dv i dv j × ϕ(v 1 , .
. . , y i , . . . , y j , . . . , v N ) dv 1 . . . dy i , . . . dy j . . . , dv N ,

where we have introduced the variables

y i = w i v i,j and
y j = w j v i,j .
Now changing the names of variables, we finally have

Q * F (v 1 , . . . , v N ) = 2 N(N -1) i<j T 2 F (v 1 , . . . , y i , . . . , y j , . . . , v N )g v i y * i,j g v j y * i,j dy i dy j , (2.7)
where

y i,j = y i + y j |y i + y j | .
Therefore, we can write

Q * = N 2 -1 i<j Q * (i,j) ,
with Q * (i,j) defined by (2.6). Then, using (2.2), we get eq. (2.5).

Remark 2.2 It is useful to note that the adjoint L of L * which corresponds to the Markov transition operator, is defined by

L = N N 2 -1 i<j (Q (i,j) -I) , (2.8)
with

Q (i,j) ϕ( v) = T 2 ϕ(v 1 , . . . , w i v i,j , . . . , w j v i,j , . . . , v N ) g(w i ) g(w j ) dw i dw j .

The CL dynamics

We now derive the master equation for the CL model. We introduce the notation (v 1 , . . . , v i , . . . , v N ) for the n -1 tuple formed by removing v i from v and

[F ] i (v 1 , . . . , v i , . . . , v N ) := S 1 F (v 1 , . . . , v N )dv i
for the marginal of F obtained by integrating v i out. We show the Proposition 2.3 The master equation of the CL dynamics is written (2.3) with L * given by (2.5) and the binary interaction operator Q * (i,j) by:

Q * (i,j) F ( v) = 1 2 [F ] i (v 1 , . . . , v i , . . . , v N ) + [F ] j (v 1 , . . . , v j , . . . , v N ) g(v * i v j ) .
(2.9)

Proof: From the definition (2.5) of the CL process, the Markov transition operator Q is given by:

Qϕ( v) = 1 N(N -1) i<j S 1 [ϕ(v 1 , . . . , zv j , . . . , v j , . . . , v N ) +ϕ(v 1 , . . . , v i , . . . , zv i , . . . , v N )] g(z) dz . (2.10)
To compute the adjoint of Q, we note that for any probability density F on T N ,

T N F (v 1 , . . . , v N ) S 1 ϕ(v 1 , . . . , zv j , . . . , v j , . . . , v N )g(z)dz dv 1 . . . dv N = = T N-1 [F ] i (v 1 , . . . , v i , . . . , v N ) S 1
ϕ(v 1 , . . . , zv j , . . . , v j , . . . , v N ) g(z)

dz dv 1 . . . dv i , . . . , dv N (2.11)

We next introduce a new variable y i = zv j (or equivalently z = v * j y i ). Evidently dz = dy i . Additionally, since v i has been integrated out, it has disappeared from (2.11). Therefore, we can change the name y i into v i without any confusion. We can then rewrite (2.11) as

T N F (v 1 , . . . , v N ) S 1 ϕ(v 1 , . . . , zv j , . . . , v j , . . . , v N )g(z)dz dv 1 . . . dv N = = T N [F ] i (v 1 , . . . , v i , . . . , v N ) g(v * j v i ) ϕ(v 1 , . . . , v i , . . . , v j , . . . , v N ) dv 1 . . . dv i , . . . , dv N . (2.12)
Using this formula and its analog for i and j exchanged, we see that the master equation for the CL dynamics is given by (2.2), (2.5) with Q (i,j) given by (2.9).

Remark 2.3 Again, we note that the adjoint L of L * is defined by (2.8) with Q (i,j) , the adjoint of Q * (i,j) , given by:

Q (i,j) ϕ( v) = S 1
ϕ(v 1 , . . . , zv j , . . . , v j , . . . , v N )g(z)dz .

Extension: Pair-Interaction driven Master Equation

The master equations of the BDG and CL dynamics are two examples of a class of master equations which we will call 'Pair Interaction driven Master Equations', defined below.

Definition 2.4 (Pair Interaction Driven Master Equation) A pair interaction driven

Master equation is an equation of the form

∂ ∂t F ( v, t) = L * F ( v, t),
describing the evolution of probability densities on some product space

X N with elements v = (v 1 , . . . , v N ) where L * = N i<j p i,j ( v) (Q * (i,j) -I) .
The operators Q (i,j) are Markov operators on functions on X N such that Q (i,j) ϕ = ϕ whenever ϕ does not depend on either v i or v j . The pair selection probabilities p i,j ( v) are such that p i,j ( v) ≥ 0 and

i<j p i,j ( v) = 1 .
We have given two examples already: the CL and BDG master equations. The Kac Master equation [START_REF] Kac | Foundations of kinetic theory[END_REF] is another example and is described below.

Example 2.1 The Kac Master equation. In this example, X N is the sphere in R N of radius √ N ,

Q (i,j) ϕ( v) = π -π ρ(θ)ϕ(R i,j,θ v)dθ , R i,j,θ v = (v 1 , v 2 , . . . , cos θ v i + sin θ v j , . . . , -sin θ v i + cos θ v j , . . . , v n ), p i,j = 2 N(N -1)
and ρ is a probability density on S 1 . The operators Q (i,j) in the Kac model are self adjoint with respect to the uniform probability measure on the sphere S N -1 , which is therefore the invariant measure for this process. In other words the Kac process is reversible meaning that is satisfies detailed balance: if you saw a movie of the process running backwards, there would be no clue that it was running backwards.

By contrast to the Kac model, the BDG and CL models do not have detailed balance property and time reversibility. If you ran the movie backwards, you would see pairs of fish with similar velocities changing them to differ in a random way. For these non-reversible processes, it is not so easy to determine the invariant measure, though it will exist and be unique for each N for our processes under mild assumptions on the noise distribution g. In section 4, it will be possible to determine the marginals of the invariant density of the CL dynamics in closed form. However, this simplification is not possible for the BDG dynamics.

In the next section, we show that pair interaction driven master equations with uniform selection probabilities p i,j = 2/(N(N -1)) do have the propagation of chaos property, and therefore, satisfy a kinetic equation at the kinetic time scale. However, in section 4, we show that the equilibrium density of the CL dynamics cannot satisfy the propagation of chaos property, meaning that this property may break down at larger time scales.

Propagation of Chaos 3.1 Definition

To pass to a kinetic description, and then on to a hydrodynamic description, the key step is a propagation of chaos result. That may seem unlikely in the cases of the CL and BDG dynamics which are expected to build pair correlations. However, the time scales at which pair correlations built up may be longer than the kinetic time scale at which a kinetic model is expected to be valid. In the present section, we shall see that chaos is propagated in both the BDG and CL models at the kinetic time scale. In section 4, we prove that the invariant measure of the CL dynamics is not chaotic; it exhibits pair correlation. These two observations are not self-contradictory since propagation of chaos holds only on a finite time scale while the invariant measure is reached as time tends to infinity. Definition 3.1 (Chaos) Let X N be the N-fold cartesian product of a polish space X equipped with some reference measure µ. Let f be a given probability density on X. For each N ∈ N, let F (N ) be a probability density on X N with respect to µ ⊗N . The sequence

{F (N ) } N ∈N of probability densities on X N is f -chaotic in case (1.) Each F (N ) is a symmetric function of {v 1 , v 2 , . . . , v N } (2.

) For each fixed k, and any bounded measurable function

φ on R k , lim N →∞ X N φ(v 1 , v 2 , . . . , v k )F (N ) ( v)dσ = X k φ(v 1 , v 2 , . . . , v k ) k j=1 f (v j )d k v .
Kac [START_REF] Kac | Foundations of kinetic theory[END_REF] proved that the semigroup e tL * associated to Kac's master equation propagates chaos. More precisely, Kac's Theorem is stated as follows:

Theorem 3.2 (Propagation of chaos) Let {F (N ) } N ∈N be f -chaotic. Then the family {e tN (Q-I) F (N ) } N ∈N is f (•, t)-chaotic where f (v, t) is the solution of ∂f ∂t (v, t) = Q(f, f )(v, t) with f (v, 0) = f (v) , (3.1
)

with Q(f, f ) = 2 R π -π [f (v ′ )f (w ′ ) -f (v)f (w)] ρ(θ) dθdw ,
and v ′ = cos θ v + sin θ w, w ′ = -sin θ v + cos θ w.

Eq. (3.1) is called the Kac-Boltzmann equation. In this section we prove a propagation of chaos result valid in the general class of pair interaction driven Master equations. We shall use this result to discuss the kinetic limits of the BDG and CL dynamics.

Propagation of chaos for pair interaction driven master equations

Consider a general Master equation

∂ ∂t F = L * F, (3.2)
for a probability density F on T N of the form

L * F = N(Q * -I)F = 2 N -1 i<j (Q * (i,j) -I)F, (3.3)
where Q i,j is a Markovian operator acting on F through v i and v j alone. The goal of this section is to prove the following:

Theorem 3.3 Let {F (N )
0,N } N ∈N be f -chaotic. Then for each t > 0, the family of marginals {e tL * F (N ) 0,N } N ∈N associated to eq. (3.2), where L * is a pair-interaction operator of the form (3.3), is f (•, t)-chaotic where f (•, t) satisfies the following Boltzmann equation:

∂ ∂t f (v, t) = 2 S 1 Q * (1,2) f ⊗2 (v, w) dw -f (v, t) , (3.4
)

associated to the initial condition f (v, 0) = f .
Before proving Theorem 3.3, we make some preliminary comments. Let initial data F 0,N be given, and let us compute the evolution of F

t , the single particle marginal at time t. For any test function ϕ(v) of the single coordinate v ∈ S 1 , we have

S 1 ϕ(v) F (1) t (v) dv = T N ϕ(v 1 ) e tL * F 0,N ( v) dv 1 . . . dv N = T N e tL ϕ(v 1 ) F 0,N ( v) dv 1 . . . dv N .
A similar relation holds for the two particle marginal and so on. So, to study the evolution of low dimensional marginals, it is helpful to understand the behavior of expressions of the form e tL ϕ( v) when ϕ( v) depends on only finitely many coordinates in v.

It is clear that in general, for a bounded continuous function ϕ on T N ,

Lϕ ∞ ≤ 2N ϕ ∞ .
However, if ϕ depends only on v 1 , . . . , v k , tighter bounds are valid. This is because

i, j > k ⇒ Q (i,j) ϕ = ϕ ,
and so

Lϕ = 2 N -1 i<j (Q (i,j) -I)ϕ = 2 N -1 k i=1 N j=i+1 (Q (i,j) -I)ϕ , (3.5)
and thus, as soon as k ≥ 1:

Lϕ ∞ ≤ 2 N -1 k(N - k + 1 2 ) 2 ϕ ∞ ≤ 4k ϕ ∞ . (3.6)
We can now state the following fundamental lemma:

Lemma 3.4 Let ϕ be a function depending only on v 1 , . . . , v k . We can regard ϕ as a function on X N for each N ∈ N, N ≥ k. Then, the power series

e tL ϕ = ∞ k=0 t k k! L k ϕ, (3.7)
converges absolutely in L ∞ , uniformly in N ∈ N * and t ∈ [0, T ] for any T < 1/4.

Proof: Consider first the case in which ϕ depends only on one variable. Without loss of generality, owing to the permutation symmetry of the problem, we can set this variable to v 1 . Then from (3.5), Lϕ is an average of functions depending on only two velocities. Likewise, L 2 ϕ is a combination of terms depending only on three velocities and so on. By what we have noted above, we can expect the following formula:

L k ϕ ≤ 4 k k! ϕ ∞ . (3.8)
To show (3.8) we prove that L k ϕ is of the form

L k ϕ = 2 N -1 k s∈S k ψ (k) s , (3.9)
where the set S k is a set of multi-indices s = (1, s 1 , . . . , s k ), such that

Card S k ≤ k! k j=1 (N - j + 1 2 ). (3.10)
The function ψ (k) s depends only on the k + 1 variables (v 1 , v s 1 , . . . , v s k ) and satisfies

ψ (k) s ∞ ≤ 2 k ϕ ∞ . (3.11)
Of course, (3.8) results from (3.9), (3.10) and (3.11) and the remark that

1 (N -1) k k j=1 (N - j + 1 2 ) ≤ 1.
We note that, if k + 1 > N, some of the variables (v 1 , v s 1 , . . . , v s k ) may be identical, but the arguments below are still valid in this case. The proof of (3.9) is by induction. For k = 1, using (3.5), we have

Lϕ = 2 N -1 N j=2 (Q (1,j) -I)ϕ.
Therefore, letting S 1 = {2, . . . , N} and ψ

L k+1 ϕ = 2 N -1 k s∈S k Lψ (k) s = 2 N -1 k+1 s∈S k m∈s,ℓ>m (Q (m,ℓ) -I)ψ (k)
s .

The expression m ∈ s means that m takes any of the indices {1, s 1 , . . . , s k } present in the multi-index s. We know from the computation in formula (3.6) that there are (k + 1)(N -k+2 2 ) such pairs (m, ℓ) for a given s. Defining S k+1 as the set of the soconstructed multi-indices {s ′ = (s, ℓ)}, we can write

L k+1 ϕ = = 2 N -1 k+1 s ′ ∈S k+1 ψ (k+1) s ′ , (3.12)
with ψ

(k+1) s ′ = (Q (m,ℓ) -I)ψ (k) s .
Clearly,

ψ (k+1) s ′ is a function of the k + 2 variables (v 1 , v s 1 , . . . , v s k , v m ) and we have Card S k+1 = (k + 1)(N - k + 2 2 ) Card S k = (k + 1)(N - k + 2 2 ) k! k j=1 (N - j + 1 2 ) = (k + 1)! k+1 j=1 (N - j + 1 2 ). (3.13)
Finally, we have, using the Markov property

ψ (k+1) s ′ ∞ ≤ 2 ψ (k) s ∞ ≤ 2 k+1 ϕ ∞ . (3.14)
Now, collecting (3.12), (3.13), (3.14) shows that the induction hypothesis is valid at rank k + 1. This proves (3.8). Using (3.8), we deduce that

t k k! L k ϕ ∞ ≤ (4t) k ϕ ∞ ,
uniformly in N, and so for t < 1/4, we have the uniform absolute convergence of the series (3.7).

If now ϕ is a function p variables, p ≥ 2, formula (3.8) is changed into

L k ϕ ≤ 4 k k! p + k -1 k ϕ ∞ (3.15) ≤ 4 k k! (k + p -1) p-1 (p -1)! ϕ ∞ ≤ C p 4 k k! (k + 1) p-1 ϕ ∞ . (3.16)
The proof follows the same lines as above. The only thing to note is that Card S k is now changed into

Card S k = k j=1 (p + j -1)(N - p + j 2 ),
the other expressions remaining identical. Then, we get

t k k! L k ϕ ∞ ≤ C p (4t) k (k + 1) p-1 ϕ ∞ , (3.17)
uniformly in N. For any given p ≥ 2, the right-hand side of (3.17) is still the general term of a convergent series for t < 1/4. Therefore, the series (3.7) is still absolutely uniformly convergent for t < 1/4, which ends the proof of Lemma 3.4.

Remark 3.1 ¿From the last proof, we note that the series (3.7) is not uniformly convergent with respect to p.

We now have the following Lemma 3.5 Suppose that F 0,N is a symmetric probability density on T N , and suppose that ϕ (k) depends only on v 1 , . . . , v k and is L ∞ . Define ϕ (k+1) by

ϕ (k+1) (v 1 , . . . , v k+1 ) = 2 k i=1 (Q (i,k+1) -I)ϕ (k) (v 1 , . . . , v k ) . (3.18)
Then, if N ≥ k + 1, we have

T N F 0,N Lϕ (k) dv 1 . . . dv N = T N F 0,N (ϕ (k+1) + φ(k+1) ) dv 1 . . . dv N , (3.19)
where φ(k+1) only depends on v 1 , . . . , v k+1 and is such that

φ(k+1) ∞ ≤ 6 k (k -1) N -1 ϕ (k) ∞ . ( 3

.20)

Proof: We compute, using (3.5) and that N ≥ k + 1:

T N F 0,N Lϕ (k) dv 1 . . . dv N = = 2 N -1 k i=1 N j=i+1 T N F 0,N (Q (i,j) -I)ϕ (k) dv 1 . . . dv N = 2 N -k N -1 k i=1 T N F 0,N (Q (i,k+1) -I)ϕ (k) dv 1 . . . dv N + + 2 N -1 i<j≤k T N F 0,N (Q (i,j) -I)ϕ (k) dv 1 . . . dv N = k i=1 T N F 0,N (ϕ (k+1) + φ(k+1) ) dv 1 . . . dv N ,
where we have used the symmetry in the second equality, and where

φ(k+1) (v 1 , . . . , v k+1 ) = 2 k -1 N -1 k i=1 (Q (i,k+1) -I)ϕ (k) + + 2 N -1 i<j≤k (Q (i,j) -I)ϕ (k) . (3.21)
This shows (3.19). Now, from the Markov property of Q (i,j) and from (3.21), we get (3.20), which ends the proof of the Lemma. Now consider any ϕ (m) depending only on v 1 , . . . , v m . Since if F 0,N is symmetric, so is each (L * ) k F 0,N . Therefore, we can repeatedly apply the previous lemma and so on. Using (3.18), we inductively define ϕ (m+1) , ϕ (m+2) , . . . ϕ (m+k) . We note that, if ϕ (m) does not depend on N, neither does ϕ (m+k) . Therefore, the functions ϕ (m+k) are good candidates to express what happens in the limit N → ∞. However, for a given N, formula (3. [START_REF] Graham | Convergence rate on path space for stochastic particle approximations to the Boltzmann equation[END_REF]) is only valid until m + k ≤ N. A special treatment is required for indices k such that m + k > N. The contribution of the remainder φ(k+1) needs also to be estimated. These are the goals of the following Lemma.

Lemma 3.6 We assume that {F 0,N } is f -chaotic (see definition 3.1). Then, for t < 1/4, we have:

lim N →∞ T N F 0,N e tL ϕ (m) dv 1 . . . dv N = = ∞ k=0 t k k! T k+m k+m j=1 f (v j ) ϕ (m+k) (v 1 , . . . , v k+m ) dv 1 . . . dv k+m . (3.

22)

Proof: ¿From Lemma 3.4, we have lim

N →∞ T N F 0,N e tL ϕ (m) dv 1 . . . dv N = lim H→∞ lim N →∞ S N H (t) , with S N H (t) = H k=0 t k k! T N F 0,N L k ϕ (m) dv 1 . . . dv N , (3.23)
Indeed, according to Lemma 3.4, the series (3.23) converges absolutely, uniformly with respect to N and we can interchange the H → ∞ and N → ∞ limits. Now, for k+m ≤ N, according to the inductive definition of ϕ (m+k) and to (3.19), we can write:

T N F 0,N L k ϕ (m) dv 1 . . . dv N = T N F 0,N ϕ (m+k) + k j=1 L k-j φ(m+j) dv 1 . . . dv N .
Taking an index H such that H + m ≤ N, we have:

S N H (t) = H k=0 t k k! T N F 0,N ϕ (m+k) dv 1 . . . dv N + T N F 0,N R N H ϕ (m) dv 1 . . . dv N , (3.24) with R N H ϕ (m) = H k=0 t k k! k j=1 L k-j φ(m+j) .
Now, taking successively N → ∞ then H → ∞, we show that the first term of (3.24) tends to the right-hand side of (3.22) and the second one tends to zero. We start with the first term. By the chaos assumption, we have, as N → ∞:

lim N →∞ H k=0 t k k! T N F 0,N ϕ (m+k) dv 1 . . . dv N = = H k=0 t k k! T m+k k+m j=1 f (v j ) ϕ (m+k) dv 1 . . . dv N . (3.25)
We show that the series at the right-hand side of (3.25) is absolutely convergent, uniformly with respect to t in any interval [0, T ] with T < 1/4. The proof follows the same lines as that of Lemma 3.4 and we only sketch it. Using (3.18) and the Markov property of Q (i,k+1) , we have

ϕ (m+k) ∞ ≤ 4(m + k -1) ϕ (m+k-1) ∞ ≤ 4 k k! m + k -1 k ϕ (m) ∞ (3.26) ≤ C m 4 k k! (k + 1) m-1 ϕ (m) ∞ .
Therefore, since k+m j=1 f (v j ) is a probability:

t k k! T m+k k+m j=1 f (v j ) ϕ (m+k) dv 1 . . . dv N ≤ C m (4t) k (k + 1) m-1 ϕ (m) ∞ .
This is the general term of a convergent series which converges uniformly as stated above.

We deduce that lim

H→∞ lim N →∞ H k=0 t k k! T N F 0,N ϕ (m+k) dv 1 . . . dv N = = ∞ k=0 t k k! T m+k k+m j=1 f (v j ) ϕ (m+k) dv 1 . . . dv N ,
which is the right-hand side of (3.22). We now consider the second term of (3.24). First, using (3.15) with the pair (k, p) replaced by (k -j, m + j), we get

L k-j φ(m+j) ∞ ≤ 4 k-j (m + k -1)! (m + j -1)! φ(m+j) ∞ . (3.27)
Then, (3.20) with k replaced by m + j -1 yields

φ(m+j) ∞ ≤ 6 (m + j -1) 2 N -1 ϕ (m+j-1) ∞ . (3.28)
Finally, using (3.26) with k = j -1, we obtain

ϕ (m+j-1) ∞ ≤ 4 j-1 (m + j -2)! (m -1)! ϕ (m) ∞ . (3.29)
Collecting (3.27), (3.28) and (3.29) leads to

L k-j φ(m+j) ∞ ≤ 3 2 4 k (m + k -1)! (m -1)! m + j -1 N -1 ϕ (m) ∞ ,
and:

k j=1 L k-j φ(m+j) ∞ ≤ 3 2 4 k (m + k -1)! (m -1)! (m + k) 2 2(N -1) ϕ (m) ∞ ≤ C m 4 k k! (k + 1) m+1 N ϕ (m) ∞ ,
which finally gives:

t k k! k j=1 L k-j φ(m+j) ∞ ≤ C m (4t) k (k + 1) m+1 N ϕ (m) ∞ ,
Therefore, we have

T N F 0,N R N H ϕ (m) dv 1 . . . dv N ≤ 1 N C m ϕ (m) ∞ ∞ k=0 (4t) k (k + 1) m+1 ,
and deduce that lim

H→∞ lim N →∞ T N F 0,N R N H ϕ (m) dv 1 . . . dv N = 0,
which ends the proof.

We now have all the elements to prove Theorem 3.3.

Proof of Theorem 3.3: We need to show the existence of a function of a single velocity variable f (v, t) such that lim

N →∞ T N F 0,N e tL ϕ (m) dv 1 . . . dv N = Tm m j=1 f (v j , t) ϕ (m) dv 1 . . . dv N , (3.30)
for all functions ϕ (m) of m velocity variables (v 1 , . . . , v m ) and all m ≤ 1.

First applying this definition with m = 1 defines f (v, t) (if it exists), by duality:

T 1 f (v, t) ϕ(v) dv = lim N →∞ T N F 0,N e tL ϕ dv 1 . . . dv N . (3.31) Now, applying (3.30
) with m = 2, ϕ (2) (v 1 , v 2 ) = η (1) (v 1 )ξ (1) (v 2
), and using (3.31) leads to lim

N →∞ T N F 0,N e tL ϕ (2) dv 1 . . . dv N = lim N →∞ T N F 0,N e tL η (1) dv 1 . . . dv N lim N →∞ T N F 0,N e tL ξ (1) dv 1 . . . dv N . (3.32)
Reciprocally, (3.32) implies (3.30) for m = 2 and general functions ϕ (2) (v 1 , v 2 ) by the density of linear combinations of tensor products. Finally, the same argument applied to arbitrary m would require to check the property for ϕ (m) equal to the tensor product of m one-dimensional functions. This proof follows closely the proof for the case m = 2 and will be omitted. So, we now focus to the case m = 2. We first consider the case of small t < 1/4. Then, using (3.22), eq. (3.32) is equivalently written:

∞ K=0 t K K! T K+2 K+2 j=1 f (v j ) ϕ (K+2) dv 1 . . . dv K+2 = ∞ k=0 t k k! T k+1 k+1 j=1 f (v j ) η (k+1) dv 1 . . . dv k+1 × × ∞ ℓ=0 t ℓ ℓ! T ℓ+1 ℓ+1 j=1 f (v j ) η (ℓ+1) dv 1 . . . dv ℓ+1 .
We start with the right-hand side of this formula, denoted by R. By distributing the various terms in the product, we can write

R = ∞ K=0 t K K! T K+2 K+2 j=1 f (v j ) K m=0 K m η (m+1) (v 1 , v 3 , . . . , v m+2 )ξ (K-m)+1 (v 2 , v m+3 , . . . , v K+2 ) dv 1 . . . dv K+2 .
Therefore, the result is proved if we show that for any symmetric function of F (v 1 , v 2 , . . . , v K+2 ), we have

T K+2 K m=0 K m η (m+1) (v 1 , v 3 , . . . , v m+2 )ξ (K-m)+1 (v 2 , v m+3 , . . . , v K+2 ) F (v 1 , v 2 , . . . , v K+2 ) dv 1 . . . dv K+2 = = T K+2 ϕ (K+2) (v 1 , v 2 , . . . , v K+2 ) F (v 1 , v 2 , . . . , v K+2 ) dv 1 . . . dv K+2 . (3.33)
Indeed, at the initial step K = 1, by direct computation we get:

T 3 ϕ (3) (v 1 , v 2 , v 3 ) F (v 1 , v 2 , v 3 ) dv 1 . . . dv 3 = = T 3 (η (2) (v 1 , v 3 )ξ (1) (v 2 ) + η (1) (v 1 )ξ (1) (v 2 , v 3 )) F (v 1 , v 2 , v 3 ) dv 1 . . . dv 3 .
Then, (3.33) is easily proved by induction, using elementary properties of the binomial coefficients. This shows that e tL * F 0,N is f (•, t)-chaotic for small t < 1/4. It remains to show that f (v, t) is a solution of (3.4). Again, assuming small t < 1/4 and Using (3.22) with m = 1 and (3.31), we get

T 1 f (v, t) ϕ (1) (v) dv = ∞ k=0 t k k! T k+1 k+1 j=1 f (v j ) ϕ (k+1) (v 1 , . . . , v k+1 ) dv 1 . . . dv k+1 .
The convergence of the series at the right-hand side is uniform for t < 1/4. So, we can differentiate this formula with respect to t and obtain:

T 1 ∂f ∂t (v, t) ϕ (1) (v) dv = = ∞ k=0 t k k! T k+2 k+2 j=1 f (v j ) ϕ (k+2) (v 1 , . . . , v k+2 ) dv 1 . . . dv k+2 . (3.34)
Applying (3.22) and (3.30) with m = 2, we can re-write the right-hand side of (3.34) and get

T 1 ∂f ∂t (v, t) ϕ (1) (v) dv = T 2 2 j=1 f (v j , t) ϕ (2) dv 1 dv 2 = 2 T 2 f (v 1 , t)f (v 2 , t) (Q (1,2) -I)ϕ (1) dv 1 dv 2 = 2 T 2 (Q * (1,2) -I)(f ⊗ f ) ϕ (1) dv 1 dv 2 ,
which is the weak form of (3.4). This shows that f (v, t) is a weak solution of (3.4) for small t < 1/4. In the last step of the proof, we need to remove the restriction on t < 1/4. However, since this bound is independent of the initial data, we can partition any interval [0, T ] by intervals [t k , t k+1 ] of length t k+1 -t k < 1/4 and apply the result on each of these intervals with initial data F N (t k ) and f (v, t l ). This shows that e tL * F 0,N is f (•, t)-chaotic on any finite-size interval [0, T ] where f (v, t) is the solution of (3.4). This ends the proof. Remark 3.2 In Theorem 3.3 the assumption that the pair selection probabilities p i,j ( v) are uniform, i.e.

p i,j = 2 N(N -1) , (3.35)
is crucial. Indeed, we need at least that (i) there exists a uniform constant C (independent of N) such that

p i,j ≤ C N(N -1)
.

This is required to get the fundamental property (3.6), with the constant 4k replaced by 4Ck.

(ii) p i,j = p i,j (v i , v j ). In this way, the key fact in the proof of Lemma 3.4 that Lϕ, for ϕ depending on only one velocity, is an average of terms only depending on two velocities and so on for L k ϕ remains true.

(iii) p i,j (v, w) = p(v, w) is independent of (i, j) to preserve the permutation symmetry of the problem.

We see that these three properties together imply that (3.35). Indeed, maintaining that

i<j p i,j = 1, (3.36)
for all velocity configurations necessitates that p is a constant, and the normalization constraint (3.36) leads to (3.35).

Remark 3.3 Theorem 3.3 can be extended to interaction processes involving multiple interactions as soon as the number of particles involved in an elementary interaction is finite and bounded independently of N. For instance, it will hold with a ternary interaction process, in which any interaction involves triples of particles. More generally, it will hold with a p-fold interaction process where any interaction involves exactly p particles. As long as the interactions involve a finite number p of interactions, with p constant or bounded by a constant P independent of N, the combinatorial arguments which have been developed above can be extended. Again, the master equation must combine the elementary interaction operators by means of uniform selection probabilities for the same arguments as those developed in Remark 3.2. Boltzmann operators with multiple interactions have been previously considered in [START_REF] Bobylev | On the self-similar asymptotics for generalized non-linear kinetic Maxwell models[END_REF].

Remark 3.4 As a by-product of Theorem 3.3, we have the existence of a weak solution of the nonlinear problem (3.4). Additionally, the proof is constructive. However, the Theorem does not prove uniqueness, as another solution could exist not being related to a solution of the corresponding master equation.

3.3 Application to the BDG and CL dynamics

The BDG dynamics

In the BDG case, thanks to (2.6), we have

Q * (1,2) f ⊗2 (v, w) = T 2 f (y 1 ) f (y 2 ) g(vy * i,2 ) g(wy * 1,2 ) dy 1 dy 2 ,
and thus

S 1 Q * (1,2) f ⊗2 (v, w) dw = T 2 f (y 1 ) f (y 2 ) g(vy * i,2
) dy 1 dy 2 .

Thus defining

Q + (f, f )(v) = T 2 f (y 1 ) f (y 2 ) g(vy * i,2 ) dy 1 dy 2 , the one particle marginal f (•, t) such that {e tL * F (N ) 0,N } is f (•, t)-chaotic satisfies the kinetic- type equation ∂f ∂t (v, t) = T 2 f (y 1 ) f (y 2 ) g(vy * i,2 ) dy 1 dy 2 -f (v, t) , (3.37)
with f (v, 0) = f (v).

The CL dynamics

In the CL case, thanks to (2.9) we have:

Q * (1,2) f ⊗2 (v, w) = 1 2 (f (v) + f (w)) g(v * w) ,
and thus

1 S Q * (1,2) f ⊗2 (v, w) dw = 1 2 f (v) + S 1 f (w) g(v * w) dw = 1 2 [f (v) + f ⋆ g(v)] ,
where ⋆ denote the convolution. Thus defining

Q + (f, f )(v) = 1 2 [f (v) + f ⋆ g(v)] , the one particle marginal f (•, t) such that {e tL * F 0,N } is f (•, t)-chaotic satisfies the kinetic- type equation ∂f ∂t (v, t)(v, t) = 1 2 [g ⋆ f (v, t) -f (v, t)], (3.38)
with f (v, 0) = f (v).

In this treatment, we have assumed that g is independent of N. But if we let the variance of g N go to zero with N, we find lim

N →∞ g N ⋆ f (v) = f (v) ,
and then we have

∂f ∂t (v, t) = 0 .
That is, chaos is propagated, but nothing at all happens on the kinetic time scale. On a much longer time scale, correlations develop and a new approach is needed to describe the bulk limit. This is what is shown in the next section by investigating the invariant densities, i.e. the equilibria of the master equation.

4

The Invariant densities F ∞,N for the CL dynamics

Preliminaries

Both the BDG and CL processes are clearly ergodic as long as g is continuous, say, and so for each there will be a unique invariant density F ∞ , i.e., a unique density

F ∞ with Q * F ∞ = F ∞ .
Since the process is symmetric under permutations of the variables, it is clear that F ∞ will be symmetric. It is not easy to write F ∞ down in closed form. However, in the case of the CL dynamics, a very special property is true, namely the hierarchy of equations for the marginals (or BBGKY hierarchy) is closed at any order. This special feature will provide more information on the marginals of F ∞ . In particular, we show that under some specific scaling of the noise with respect to N, the invariant measure is not chaotic. This is not in contradiction to Theorem 3.3, since it is valid for fixed noise and on finite time intervals. The invariant density for the CL dynamics is the function F ∞ which cancels the righthand side of (2.5) with Q * (i,j) given by (2.9). Therefore, it satisfies

F ∞ ( v) = 1 N(N -1) i<j [[F ∞ ] i (v 1 , . . . , v i , . . . , v N )+ +[F ∞ ] j (v 1 , . . . , v j , . . . , v N ) g(v * i v j ) . (4.1)
While it is not easy to write F ∞ down in closed form, we can at least say what F ∞ is not:

In general F ( v) = 1 does not solve Q * F = F , i.e.
, F ∞ is not the uniform density. Indeed, if we replace F ∞ by 1 on the right hand side of (4.1), we find

i<j 2 N(N -1) g(v * j v i ) .
This will equal 1 for all v if and only if g(z) = 1 for all z. However, it is easy to see that for fixed smooth g, and large N,

i<j 2 N(N -1) g(v * j v i ) ≈ 1 ,
with high probability if v is selected at random, uniformly on T N , and so one can expect that for fixed g, the invariant density F ∞ becomes more and more uniform as N increases.

However, a non-uniform invariant density can be found in the large N limit if one includes some N dependence in the noise density g in such a way that it more and more closely approximates a δ function. This is perhaps justified in the context of biological modeling: if the population is small, one fish may be less interested in carefully mimicking his neighbor than when the population is large. One can imagine that the larger the group, the more important it is to follow behavioral rules closely. However, biological data are needed to support this claim.

While it does not seem easy to write F down in closed form, it is possible to obtain analytical expressions of its marginals. This is the aim of the next section.

Marginals

For any symmetric density F and any m = 1, 2, N -1, define the m-variable marginal density F (m) on T m by

F (m) (v 1 , . . . , v m ) = T N-m F (v 1 , . . . , v m , v m+1 , . . . , v N ) dv m+1 . . . dv N .
For typical binary collision process master equations, the evolution of the m-variable marginal density depends on the (m + 1)-variable marginal density, and one gets a hierarchy of evolution equations, the so-called BBGKY hierarchy. To break the hierarchy at some manageable level (e.g. m = 1 in the usual kinetic theory case) in the large N limit, one generally needs special assumptions on the initial data, such as the f -chaotic property, for a suitable f . One typically runs into the same hierarchy problem when trying to compute the marginals of the invariant density, except that in some cases, one has Q = Q * , and then the invariant density is simply uniform.

For the CL dynamics, the hierarchy breaks itself. Before stating the result, we recall that the Fourier series of a given function φ(w), w ∈ S 1 is defined by φ(k) = S 1 w -k φ(w) dw, k ∈ Z.

We have:

Proposition 4.1 Let F (1)
∞ and F

(2) ∞ be the one and two-variable marginal invariant densities of the CL process. We have:

(i) F (1)
∞ is uniform, i.e.

F (1) ∞ (v 1 ) = 1, ∀v 1 ∈ S 1 .

(ii) F

∞ is given by:

F (2) ∞ (v 1 , v 2 ) = F (v * 1 v 2 ), ∀(v 1 , v 2 ) ∈ T 2
, with F (w), w ∈ S 1 given by its Fourier series

F (k) = 1 N -1 g(k) 1 - N -2 N -1 g(k) -1 , (4.2)

Noise scaling

The reason for considering a scaling of the noise intensity is the following. For fixed g ∈ L 1 (S 1 ), lim ℓ→∞ g * ℓ = 1, and since for large N, most of the weight in the average is on large values of ℓ, F will be nearly uniform for large values of N, and the correlations are washed out. Therefore, we recover here that the invariant density is nearly uniform for large values of N, a fact which has already been noticed (see section 4.1). But if g is taken to depend on N itself, this need not be the case. As a typical example, we can consider g(z * w) to be the kernel of e ∆/N on S 1 ; i.e., the heat kernel on S exists and is non trivial (i.e. not equal to the Kronecker δ k,0 ). Then, the corresponding correlation F N associated to g N through (4.2) or (4.3) satisfies:

lim N →∞ F N (k) := F ∞ (k) = 1 1 -γ(k)
.

Proof: Using g N in place of g in (4.2), we find

F N (k) = 1 N -1 g N (k) 1 - N -2 N -1 (1 + ( g N (k) -1)) -1 = 1 N -1 g N (k) 1 N -1 - N -2 N -1 ( g N (k) -1) -1 = g N (k) [1 -(N -2)( g N (k) -1)] -1 . (4.7)
The result follows from inserting (4.6).

Example 4.1 The heat kernel (4.5). In this case, we have

γ(k) = -k 2 , F ∞ (k) = 1 1 + k 2 .
Hence, the correlation function is a Lorentzian in Fourier space.

With this noise scaling, the two-variable marginal invariant density of the N-particle CL process F

(2)

∞,N (v 1 , v 2 ) is such that F (2) ∞,N (v 1 , v 2 ) → F ∞ (v * 1 v 2 ), as N → ∞
, where F ∞ is not the uniform distribution. Therefore, non-trivial correlations remain in the large N limit and in particular, the invariant density {F ∞,N } is not chaotic. This result is in marked contrast to the case studied by Kac, in which the invariant density is the uniform density on the sphere S N -1 ( √ N ). This family is well known to be G-chaotic where G(v) denotes the centered unit Gaussian on R. The lack of chaos in the invariant density might seem to be a strong obstacle to propagation of chaos. But as we have seen in Theorem 3.3, this is not the case. As already noticed, there is no contradiction between this two seemingly paradoxical results. Theorem 3.3 is valid under fixed noise and on a finite time interval. By contrast, the lack of chaos property of the invariant density for the CL model is shown under N-dependent noise intensity and in the infinite time limit. But these two results show that the chaos property can be valid for a finite time interval at the kinetic time scale and be lost at larger times.

Conclusion

We have considered a class of pair-interaction stochastic processes in an N-particle system and their associated pair interaction driven master equations. We have proved a chaos propagation theorem for this class of master equations and have used this result to study the kinetic limits of two biological swarm models, the BDG and CL processes. By investigating the invariant density of the CL process, we have shown that the chaos property may be lost at large times. This work shows that the chaos property may be true even for processes that seemingly build-up correlations but may not be uniformly valid in time. Correlation build-up manifests itself at large time scales. In order to restore the validity of kinetic theory at these large scales, new theories must must be developed. This is a fascinating and widely open area of research.

  (1,j) -I)ϕ, we can write Lϕ according to formula (3.9). Clearly, Card S 1 = N -1 in accordance to(3.10). Finally, by the Markov property, we have ψ (1) j ∞ ≤ 2 ϕ ∞ , which is consistent with (3.11). Therefore, (3.8) is proved for k = 1. Now, we assume that (3.8) is true for k and try to deduce it for k +1. By the induction hypothesis, ψ (k) s depends only on k + 1 variables so (3.5) applies and we have

1 Proposition 4 . 2

 142 at time 1/N. Then g N (k) := e -k 2 /N . Suppose the scaled noise intensity g N is such thatlim N →∞ (N -2)( g N (k) -1) := γ(k) (4.6)

Acknowledgements:

The first author gratefully acknowledges support from the Région Midi-Pyrénées in the frame of a Chaires Pierre-de-Fermat and support under U.S. National Science Foundation grant DMS 0901632. The second author acknowledges support from the ANR under contract 'CBDif-Fr' (ANR-08-BLAN-0333-01). The first and third authors would both like to thank the Institute of Mathematics, Toulouse, for its kind

The third author acknowledges support from the Swedish Research Council.

or equivalently

(4.3)

Eq. ( 4.3) defines F as an average of convolution powers of g.

Proof:

We begin with F

∞ . First of all, it is easy to see that for all i > 1 and all j , then

Therefore, by integrating (4.1) with respect to (v 2 , . . . , v N ) and using (2.4), we have:

Then, using the permutation symmetry of [F ∞ ] 1 , one finds for each j ≥ 2,

where the * denotes convolution on S 1 . Substituting it into (4.4), we find

) . The only solution of this equation for any L 1 density g is the uniform density F (1)

and

Next, for i = 1 and j = 2, using the permutation symmetry of [F ∞ ] 1 and (2.4):

and

Finally, for i = 1 and j > 2, we have

and

and for i = 2 and j > 2,

and

We can now compute the two variable marginals of both sides of (4.1), and we find

where

This simplifies to

Fourier transforming both sides, we have

)) and so

where δ i,j is the usual Kronecker symbol. It follows that F

∞ (k 1 , k 2 ) has the form F (k 1 )δ k 1 ,-k 2 , i.e. that F

with F defined by its Fourier transform according to (4.2) (owing to the evenness of ĝ, a consequence of (2.4)). This ends the proof of proposition 4.1.


Hong Liang

Aude Pauly

Luisa Riancho

Christophe Baudouin
email: frbaudouin@aol.com

Françoise Brignole-Baudouin

Toxicological evaluation of preservative-containing and preservative-free topical prostaglandin analogs on a 3D-reconstituted corneal epithelium system

Keywords: Cornea, Immunology, Inflammation, Drugs

Aims: Using an established three-dimensional (3D) toxicological model based on reconstituted human corneal epithelium (HCE), this study investigated the tolerability of four topical intraocular pressure (IOP)-lowering agents: the commercial solutions of benzalkonium chloride (BAC)-containing 0.005% latanoprost, 0.004% travoprost, 0.03% bimatoprost

containing 0.02%, 0.015%, and 0.005% BAC, respectively and the preservative-free (PF) tafluprost. Solutions of 0.01% and 0.02% BAC alone were also evaluated for comparison.

Methods: The 3D-HCEs were treated with either solutions for 24 hours followed or not by a 24-hour recovery period. We used a modified MTT procedure to assess cell viability in the HCE. Frozen sections of HCE were analyzed using fluorescence microscopy for the evaluation of apoptosis (TUNEL), inflammation (ICAM-1), and proliferation (Ki67). Corneal epithelial tight junctions (occludin and ZO-1) were also assessed by en-face confocal microscopy in response to the different eye drops.

Results:

The MTT test revealed that the cytotoxicity of antiglaucoma eye drops was primarily related to the concentration of their common BAC preservative (0.02%BAC-latanoprost > 0.015%BAC-travoprost > 0.005%BAC-bimatoprost). PF-tafluprost did not induce obvious cytotoxicity and showed the least expression of inflammatory or apoptotic markers and revealed preservation of membrane immunostaining of tight junction proteins in comparison with BAC-containing solutions.

Conclusion:

The toxicological model of 3D reconstructed corneal epithelia model confirmed the ocular surface cytotoxicity of BAC-containing antiglaucomatous solutions. Compared to the formulations containing the toxic preservative BAC, PF-tafluprost was well tolerated without inducing significant corneal epithelium deterioration.

Introduction

Tafluprost is a newly synthesized PGF2alpha-agonist without any toxic preservative [START_REF] Nakajima | New fluoroprostaglandin F(2alpha) derivatives with prostanoid FP-receptor agonistic activity as potent ocular-hypotensive agents[END_REF] .It demonstrated no cytotoxicity in human conjunctival epithelial cell lines [START_REF] Brasnu | In vitro effects of preservative-free tafluprost and preserved latanoprost, travoprost, bimatoprost in a conjunctival epithelial cell line[END_REF] and was well tolerated following short-and repeated instillations in rabbits [START_REF] Liang | Conjunctival and corneal reactions in rabbits following short-and repeated exposure to preservative-free tafluprost, commercially available latanoprost and 0.02% benzalkonium chloride[END_REF] . The problems raised by potentially toxic glaucoma treatments should really require more attention to be paid due to the increased number of patients in the world, life-span duration of treatments that require not only efficacy, but also safety, good tolerability, and optimal compliance for patients.

Supplied by SkinEthic® Laboratories, the reconstructed three-dimensional (3D) model of human corneal cells (HCE) is an appropriate alternative toxicological method to the classical Draize test in rabbits. [START_REF] Doucet | Reconstituted human corneal epithelium: a new alternative to the Draize eye test for the assessment of the eye irritation potential of chemicals and cosmetic products[END_REF] It was proven to have the same characteristic as the corneal epithelium of the human eye in morphology and thickness, and is suitable for investigating the undesired effects of ophthalmic drugs. [START_REF] Pauly | Multiple endpoint analysis of the 3D-reconstituted corneal epithelium (SkinEthicTM) after treatment with benzalkonium chloride: modified MTT procedure and new markers for the early detection of toxic damage[END_REF] The objective of this study was to investigate in this 3D-HCE system the tissue changes after the treatment of the currently available commercial prostaglandin (PG)-derived eye drops, namely the commercial solutions of 0.005% latanoprost, 0.004% travoprost, 0.03% bimatoprost (containing 0.020%, 0.015%, and 0.005% of benzalkonium chloride (BAC), respectively), and preservative-free (PF) tafluprost. We analyzed a panel of five selected biomarkers in order to assess the phenomena of apoptosis, inflammation, proliferation and tight junction impairment after a contact with the different antiglaucoma eye drops, which would contribute to a better understanding of human corneal cell reactions following topical antiglaucoma treatments.

Materials and methods

Tissue model and antiglaucoma solution treatments

The 3D-HCE model was supplied by SkinEthic® Laboratories (Nice, France) and consists of immortalized HCE cells grown vertically on a 0.5cm 2 insert permeable polycarbonate filter.

All the experimentations were conducted as previously published 5 : 30µl of sterile phosphatebuffered saline (PBS) used as negative control solution, BAC solutions at 0.020% and 0.010% used as positive controls and the commercial solutions of 0.020% BAC-containing latanoprost (Xalatan ® ; Pfizer, NY, USA), 0.015% BAC-containing travoprost (Travatan ® , Alcon, TX, USA), 0.005% BAC-containing bimatoprost (Lumigan ® , Allergan, CA, USA) or PFtafluprost (Taflotan ® , Santen Oy, Tampere, Finland) were applied on the apical surface of 3D-HCE for 24 hours (h) and 24 h followed by an additional 24h-recovery period (24h+24hrecovery). The recovery period (24 hours) was chosen in order to observe if toxic effects on HCE were reversible or not.

Six series of 3D-HCE were used for each solution: two series for cell viability MTT testing, two series for immunohistological analyses on cryosections, and two series for immunofluorescent labeling on the most superficial layers of 3D-HCE by en-face confocal microscopic analyses.

Modified MTT test

The experiments were conducted in triplicate. The 3D-HCEs were transferred to 24-well plates containing 300μl of the MTT solution at 0.5g/ml in culture medium with the application onto the apical surface. After the 3h-incubation at 37°C, the 3D-HCEs were incubated with 750 μl isopropanol. After the agitation for 2h, the results were analyzed by the optical density (OD, absorbance) at 570 nm versus OD690 nm, and they were expressed as a percentage of cell viability compared to the negative control, PBS. Analyses were performed using Safire technology (Tecan, Lyon, France).

Immunofluorescence analyses

After incubations with different solutions, the 3D-HCE samples were transferred into Petri dishes containing PBS to be separated into two pieces using a surgical scalpel. One piece of tissue was embedded in OCT® embedding medium (Tissue-Tek, Miles Inc., Elkhart, IN, USA), and frozen at -80°C for future 10μm vertical cryosections. The other piece was fixed in 4% paraformaldehyde (PFA) for 20 min before immunofluorescent labeling of the tight junction proteins occludin and ZO-1.

1/ Detection of apoptosis (TUNEL assay), inflammation (CD54) and proliferation (Ki67) on 3D-HCE cryosections

-Apoptosis, TUNEL assay: A terminal deoxynucleotidyl transferase-mediated dUTP-nick end labeling (TUNEL) kit containing TUNEL enzyme and TUNEL label (Roche Diagnostics, Meylan, France) was used to detect apoptosis in the tissue layers. The nuclei were labeled with DAPI and the cryosections were mounted in an anti-fade medium (Vectashield; Vector Laboratories, Burlingame, CA, USA).

-CD54 (ICAM-1) and Ki67 immunostainings: The samples were firstly fixed with 4% PFA for 10 min and then permeabilized with 0.01%-diluted Triton X100® (Sigma Chemical Company, Saint Louis, MO, USA) for 5 min. The following primary antibodies were incubated: the mouse anti-human ICAM-1 (IgG1, BD Biosciences, Pharmingen, San Diego, CA, USA; 1:100), the mouse anti-human Ki67 (Immunotech, Marseilles, France; 1:25), and negative isotypic control mouse IgG1 (BD Biosciences). Alexa 488 conjugated-goat antimouse IgG at a 1:500 dilution (Invitrogen-Molecular Probes, Eugene, OR, USA) was used as second antibody. The nuclei were labeled with propidium iodide (PI) and cryosections were mounted. Samples were analyzed under a laser confocal microscope equipped with digital camera (E800, PCM 2000, Nikon, Champigny-sur-Marne, France). We counted the immunopositive cells under x 20 objective of the microscope in three different areas, and the results were calculated as the average of counts, and finally presented as cells/field (x20 objective) after each treatment.

2/ Confocal immunofluorescence on entire epithelia for tight junction staining

The rabbit anti-human occludin (IgG1, Dako, Glostrup, Denmark; 1:100 dilution) and rabbit anti-human ZO-1 (IgG1, Santa Cruz Biotechnology, Santa Cruz, CA) were used for tight junction staining. Alexa 488-conjugated goat anti-mouse was used as second antibody, and finally nuclei were labeled with DAPI. Samples were then analyzed under a laser confocal microscope (E800, PCM 2000, Nikon) for detecting occludin and ZO-1 expressions.

Quantification and Statistical Analysis

TUNEL-, ICAM-1-, and Ki67-positive cells were quantified manually, using a microscopic grid on images under ×400 magnification. The results were expressed as mean cell number per field (cells/field). Standard deviations were indicated.

Control and treatment groups for MTT data and immunopositive cell counts were compared using two-way analysis of variance (ANOVA) followed by the Fisher's adjustment (Statview V; SAS Institute Inc., Cary, NC, USA).

Results

Cell viability: MTT test

PBS that was used as negative control did not affect the cell viability neither at 24h nor 24h+24h-recovery (Fig. 1). PF-tafluprost showed the same level of cell viability as PBS at 24h (90.2%) and at 24h+24h-recovery (100%) without any statistically significant differences compared to PBS. The preservative solutions of BAC at 0.010% and 0.020% induced a cell viability decrease at 24h (65.9% and 42.5% respectively, P<0.01 compared to PBS) and at 24h+24h-recovery (62.9% and 43.4% respectively, P<0.01 compared to PBS).

The preservative-containing antiglaucoma solutions induced a BAC-concentration-dependent decrease of the cell viability: at 24h, 85.5% for 0.005%BAC-bimatoprost, 75.5% for 0.015%BAC-travoprost (P<0.01 compared to PBS group), and 40.1% for 0.020%BAClatanoprost (P<0.05 compared to PF-tafluprost, 0.005%BAC-bimatoprost and 0.015%BACtravoprost groups, P<0.01 compared to PBS and P<0.05 compared to 0.010%BAC group), at 24h+24h-recovery, 83.0% for 0.005%BAC-bimatoprost, 65.5% for 0.015%BAC-travoprost (P<0.01 compared to PBS group, P<0.02 compared to PF-tafluprost group), and 53.2% for 0.020%BAC-latanoprost (P<0.05 compared to PF-tafluprost and 0.005%BAC-bimatoprost).

Immunofluorescence analyses and quantification of apoptosis (TUNEL)

Few apoptotic cells were observed after PBS (Fig. 2A) or PF-tafluprost (Fig. 2B) incubations (Fig. 2H): 17.0 cells/field for PBS and 13.3 cells/field for PF-tafluprost at 24h; 14.5 cells/field for PBS and 15.4 cells/field for PF-tafluprost at 24h+24h-recovery. In accordance with our previous study, BAC at 0.010% and 0.020% significantly increased the number of TUNELpositive cells at 24 h, and at 24 h+24 h-recovery compared to PBS group (P<0.02 for 0.010%BAC, P<0.01 for 0.020%BAC).

The number of apoptotic cells increased after the BAC-containing antiglaucoma eye drops in a BAC-concentration dependent manner at 24h. In addition, the corneal epithelial cells could not recover after the recovery period and more cells tended to undergo apoptosis with an increase of TUNEL-positive cells at 24h+24h-recovery. Bimatoprost with 0.005% BAC (Fig. 3C) induced moderate expression of apoptosis (Fig. 2H, 23.4 cells/field at 24h and 24.8 cells/field at 24h+24h-recovery), and 0.015%BAC-travoprost (Fig. 2E) with 37.0 cells/field at 24h, and 42.6 cells/field at 24h+24h-recovery (P<0.03 compared to PF-tafluprost and 0.005%BAC-bimatoprost groups at 24h and 24h+24h-recovery). The highest levels of apoptotic cells were found homogenously disseminated after the application of 0.020%BAClatanoprost (Fig. 2F) and of 0.020% BAC (Fig. 2G): 43.4 cells/field and 44.3 cells/field respectively at 24h; 68.3cells/field and 61.3 cells/field respectively at 24h+24h-recovery.

These TUNEL-positive cells were observed especially in the apical cell layers, but also in the middle epithelial layers. Latanoprost with 0.020%BAC induced significantly higher amount of apoptotic cells than did PF-tafluprost, 0.005%BAC-bimatoprost at 24h (P<0.003 for both) and 24h+24h-recovery (P<0.003 for both), as well as 0.010%BAC and 0.015%BACtravoprost at 24h+24h-recovery (P<0.003).

Immunofluorescence analyses and quantification of the inflammation marker: ICAM-1 (CD54)

PBS-treated 3D-HCE cultures (Fig. 3A, 3H) were found to express CD54 at levels of approximately 113.2cells/field at 24h, and 74.0 cells/field at 24h+24h-recovery. PF-tafluprost induced 121.5cells/field at 24h and 79.0 cells/field at 24h+24h-recovery without any statistically significant differences with PBS-treated samples. BAC 0.005%-bimatoprost induced an increase of CD54-positive cells: 155.5 cells/field at 24h and 134.0 cells/field at 24h+24h-recovery (P<0.001 compared to PBS, P<0.002 compared to 0.02%BAC-latanoprost, P<0.008 compared to 0.015%BAC-travoprost at 24h and 24h+24-recovery). This level of increased CD54 activation was also found with 0.01%BAC: 158.8 cells/field at 24h, and 128.5 cells/field at 24h+24h-recovery (P<0.001 compared to PBS, P<0.002 compared to 0.02%BAC-latanoprost, P<0.008 compared to 0.015%BAC-travoprost at 24h and 24+24recovery). For 0.005%BAC-bimatoprost and 0.010%BAC, the CD54 positive cells were located throughout all the epithelial layers (Figs. 3C and3D). However, with more higher concentrations of BAC, this number decreased, probably due to the erosion of corneal epithelium with toxic solutions: 74.0 cells/field after 0.015%BAC-travoprost treatment, 73.3 cells/field after 0.020%BAC-latanoprost treatment, and only 64.3 cells/field after 0.020% BAC treatment at 24h (for the three groups: P<0.001 compared to PBS group, P <0.005 compared to 0.010% BAC).

Immunofluorescence analyses of cell proliferation: Ki67

After PBS treatment (Fig. 4A), some proliferating cells were observed, not restricted to the basal layer at 24h or 24h+24h-recovery. PF-tafluprost induced similar aspect as did PBStreated samples. However, more numerous proliferating cells, with a greater number located in the basal layer were found after 0.005%BAC-bimatoprost (Fig. 4C), 0.01%BAC (Fig. 4D), 0.015%BAC-travoprost (Fig. 4E), 0.02%BAC-latanoprost (Fig. 4F) and 0.02% BAC (Fig. 4G). At 24h+24h-recovery, we could still find some proliferative cells after the treatments with PBS or PF-tafluprost. In this time, following treatment with the solutions containing BAC (0.010%BAC, 0.020%BAC, and the three BAC containing antiglaucoma solutions), no or very rare proliferative cells were observed, again most likely due to the deep impairment of corneal cells after toxic challenge.

En-face confocal microscopic analysis of tight junction: occludin and ZO-1

The 3D-HCEs disclosed a thin occludin immunostaining in the most superficial cells, forming a persistent green ring around, leaving a diffuse cytoplasmic staining, after treatments with PBS (Fig. 5A) or PF-tafluprost (Fig. 5B). We could still find a slight occludin expression in the apical cells with 0.005%BAC-bimatoprost (Fig. 5C). This kind of occludin expression clearly disappeared after the treatments with 0.010%BAC (Fig. 5D), 0.015%BAC-travoprost (Fig. 5E), 0.020%BAC-latanoprost (Fig. 5F) and 0.020% BAC (Fig. 5G) at 24h, leaving only dense green patches. After the 24h recovery period, no significant change in occludin protein distribution was found compared to treatment without the recovery period. The same tendency was found using ZO-1: a thin green ring was found just after the treatment with PBS (Fig. 6H) or PF-tafluprost (Fig. 6I).

Discussion

The toxicological model of 3D-reconstructed cornea epithelial model confirmed the cytotoxicity of BAC-containing solutions with a better approach of what exists in vivo or in vitro. Using a selection of pertinent biomarkers on this 3D-HCE, we could detect and analyze the presence of cell apoptosis, activation/inflammation, proliferation/turnover and cellular tight junctions after application of different test solutions. This model, while respecting the ethical guidelines of animal experimentation, specially the 3R rule (refining, reducing and replacing the use of animals), [START_REF] Manciocco | The application of Russell and Burch 3R principle in rodent models of neurodegenerative disease: the case of Parkinson's disease[END_REF] contributes to a better visualization of the human epithelial structure than the traditional monolayer cell cultures and represents a highly valuable tool between animal and cellular models for toxicity detection. The fluorescence techniques conjugated with confocal microscopy on 3D-reconstructed corneal epithelia were proven to be suitable for the investigation of toxicological markers and gave relevant results compared with the known human data. [START_REF] Pauly | Multiple endpoint analysis of the 3D-reconstituted corneal epithelium (SkinEthicTM) after treatment with benzalkonium chloride: modified MTT procedure and new markers for the early detection of toxic damage[END_REF] As on monolayer cell cultures, the 24 hours recovery experiments on 3D-HCE model allow the apoptotic process to further develop and to induce the inhibition of proliferation. The expression of the proliferation marker, KI-67 nuclear antigen decreased after these 24 hours recovery while the apoptosis-and inflammation-related markers increased. These results confirmed that BAC, depending on the time and concentration used, acts as a pro-inflammatory and pro-apoptotic agent able to irreversibly impair the normal epithelium turn-over even after BAC withdrawal. Several in vitro studies have already pointed out these deleterious effects of BAC and BAC-containing antiglaucoma drugs. MTT assay is a simple quantitative method usually used to assess the cell viability in these 3D-HCE models. [START_REF] Cotovio | In vitro assessment of eye irritancy using the Reconstructed Human Corneal Epithelial SkinEthic HCE model: application to 435 substances from consumer products industry[END_REF] It seems more reliable and robust than immunofluorescence and confocal microscopy that are known to provide visual and convincing results for the analysis of cell damage; however, these techniques are deemed to give adjunct results for the purpose of comparison of tested solutions. But, even if based on a smaller number of cells and a quantitative reading that depends on the experimenter, we proved in this study, that these techniques also allow quantification and provide a good understanding of all toxic phenomena and their impact on various biological processes. Consistent with similar analyses in monolayer cell cultures and in rabbit model, in the present study using this new 3D model, we demonstrated that PF-tafluprost presented almost the same aspects as did the negative control with persistent epithelial cell tight junctions and appropriate proliferative level, with no induction of cellular apoptosis or inflammation when compared to the formulations containing toxic preservative BAC. BAC-induced ocular toxicity has been well-documented in cell culture, animal models and in patients. [START_REF] Baudouin | Preservatives in Eyedrops: the Good, the Bad and the Ugly. Progress in retinal and eye research[END_REF] The reconstituted 3D-HCE model, resembling the corneal epithelium of the human eye in morphology and thickness, [START_REF] Nguyen | Three-Dimensional Construct of the Human Corneal Epithelium for In Vitro Toxicology[END_REF] was already used as a new alternative to the classic Draize eye test for the assessment of the eye irritation potential of chemicals and cosmetic products. [START_REF] Doucet | Reconstituted human corneal epithelium: a new alternative to the Draize eye test for the assessment of the eye irritation potential of chemicals and cosmetic products[END_REF][START_REF] Doucet | A new in vitro human epithelial model for assessing the eye irritating potential of formulated cosmetic products[END_REF] It confirmed the cytotoxicity of BAC, and this model was sensitive enough to distinguish between different concentrations of BAC, from 0.001% to 0.5%. [START_REF] Pauly | Multiple endpoint analysis of the 3D-reconstituted corneal epithelium (SkinEthicTM) after treatment with benzalkonium chloride: modified MTT procedure and new markers for the early detection of toxic damage[END_REF] The versatility and sensitivity of the 3D-HCE models were also investigated and this model was proposed to replace some uses of animal testing in pre-clinical studies. In addition, the occludin gene expression was proposed as an early in vitro sign for mild eye irritation assessment. [START_REF] Meloni | Occludin gene expression as an early in vitro sign for mild eye irritation assessment[END_REF] Even though multiple studies consistently showed BAC toxicity and potential advantages of BAC-free compounds, [START_REF] Baudouin | Preservatives in Eyedrops: the Good, the Bad and the Ugly. Progress in retinal and eye research[END_REF] there still exist some debate on its cytotoxicity in ocular tissues and potential harmfulness for glaucoma patients. In contrast with our results and all literature data on BAC, [START_REF] Baudouin | Preservatives in Eyedrops: the Good, the Bad and the Ugly. Progress in retinal and eye research[END_REF][START_REF] Baudouin | Expression of inflammatory membrane markers by conjunctival cells in chronically treated patients with glaucoma[END_REF][START_REF] Baudouin | Conjunctival epithelial cell expression of interleukins and inflammatory markers in glaucoma patients treated over the long term[END_REF][START_REF] Saint | Effects of benzalkonium chloride on growth and survival of Chang conjunctival cells[END_REF][START_REF] Pisella | Conjunctival proinflammatory and proapoptotic effects of latanoprost and preserved and unpreserved timolol: an ex vivo and in vitro study[END_REF][START_REF] Asbell | Effects of topical antiglaucoma medications on the ocular surface[END_REF][START_REF] Leung | Prevalence of ocular surface disease in glaucoma patients[END_REF][START_REF] Hong | Allergy to ophthalmic preservatives[END_REF][START_REF] Guenoun | In vitro study of inflammatory potential and toxicity profile of latanoprost, travoprost, and bimatoprost in conjunctiva-derived epithelial cells[END_REF] one recent study did not find any loss of viability when testing BAC alone or BAC-containing latanoprost after a short contact duration from 10 to 60 minutes by using the same model of 3D-corneal epithelium as used in this study. [START_REF] Khoh-Reiter | Evaluation of the cytotoxic effects of ophthalmic solutions containing benzalkonium chloride on corneal epithelium using an organotypic 3-D model[END_REF] These conflicting results were principally due to the differences in experimental conditions, such as the choice of model or appropriateness of duration of contact. In patients, the age, the quality of tear film, the individual sensitivity of the ocular surface, a preexisting sensitivity or allergy, and the overall number of BAC-containing eye drops and duration of administration could all be considered to influence the impact of toxic preservative on the ocular tissues and may explain wide individual variations. In experimental conditions, the contact duration of BAC is essential to consider and may explain some apparently conflicting results: Khoh-Reiter and Jessen stated in their study that BAC is in contact with the eye structures only for few seconds after instillation. Another study proved that a rapid dilution in the tear film reduced the concentration of BAC to near zero in minutes. [START_REF] Friedlaender | The dilution of benzalkonium chloride (BAC) in the tear film[END_REF] However, Champeau et al. detected the presence of BAC in the conjunctival epithelium nine days after a single administration, which suggests prolonged impregnation of BAC in ocular tissue instead of simple dilution in tear film. [START_REF] Champeau | Effect of ophthalmic preservatives on the ocular surface: conjunctival and corneal uptake and distribution of benzalkonium chloride and chlorhexidine digluconate[END_REF] In rabbits, radiolabelled BAC was found in the palpebral and bulbar conjunctiva and all corneal layers, and single-drop administration resulted in high tissue levels in the anterior ocular tissues retaining for up to 120 hours. [START_REF] Green | Benzalkonium Chloride Kinetics in Young and Adult Albino and Pigmented Rabbit Eyes[END_REF] In our model of 3D-culture of corneal epithelium, a 24-hour incubation was therefore considered to better mimic the accumulation of BAC over the long term and allowed to better investigate the cell alterations in the deeper layers of the corneal epithelium. As BAC is usually considered a potent and useful enhancer of drug penetration into deep ocular structures, all experimental models should be considered in light of a significant contact time, depending not on tear concentrations of BAC but instead on much more prolonged contact times, as most likely occurs with a compound that impairs the tear film, disrupts cell membranes, increases corneal epithelium permeability, and is repeatedly used over extremely long periods of time. [START_REF] Baudouin | Preservatives in Eyedrops: the Good, the Bad and the Ugly. Progress in retinal and eye research[END_REF] Complementarily, the 24h recovery appeared justified as for the monolayer cultures for evaluating late effects, either in the way of possible compensation and recovery or as an amplification of toxic reactions.

In conclusion, our results using the new 3D-HCE model for evaluating BAC-containing prostaglandin analogs and PF-tafluprost were in accordance with previously published works investigating BAC toxicity. The present study provides a valuable method to evaluate the impacts of eye drops in vitro, being more accurate than single monolayers and easier to perform than animal models. The development of 3D-HCEs thus enriches the existing models to study the bioactivity of eye drops in the cornea by investigating and quantifying a set of cellular functions like cellular viability, inflammation, proliferation, apoptosis or cell junctions.

Figures and legends

Fig. 1: Cell viability MTT test

Cell viability of PBS-, PF-tafluprost-, 0.005%BAC-bimatoprost-, 0.010%BAC-, 0.015%BAC-travoprost-, 0.020%BAC-latanoprost-and 0.020% BAC-treated 3D-HCE after 24h-incubation without or with a 24-h recovery period. Antiglaucoma treatments induced a BAC-dose dependent decrease in cell viability as assessed by the MTT test. PF-tafluprost presented the higher level of viability than all BAC-containing eye drops. * P <0.01 compared to PBS at the same time point; # P <0.03 compared to 0.010% BAC at the same time point; $ P <0.002 or P <0.03 ($$) compared to 0.020% BAC at the same time point; ♦ P <0.05 compared to 0.020%BAC-latanoprost at the same time point; • P <0.02 compared to 0.015%BAC-travoprost at the same time point.

Fig. 2: Apoptosis analysis

Localization of TUNEL positive cells (green) on 3D-HCE samples after PBS (A), PFtafluprost (B), 0.005%BAC-bimatoprost (C), 0.010%BAC (D), 0.015%BAC-travoprost (E), 0.020%BAC-latanoprost (F) and 0.020% BAC (G) incubation after 24h without (left column) or with the 24h-recovery period (right column). Quantification of TUNEL-positive cells (H) was performed after 24 h of treatment with all solutions followed or not by a 24h-recovery period.

At 24h or 24h+24h-recovery, no or very rare apoptotic cells were observed after PBS (A) or PF-tafluprost (B) treatments. Bimatoprost with 0.005% BAC (C) induced moderate expression of apoptosis, and 0.010%BAC (D), 0.015%BAC-travoprost (E), 0.020%BAClatanoprost (F) 0.020% BAC (G) all induced a great number of TUNEL-positive cells principally in the apical cell layers, and also in the middle epithelial layers. The nuclei were stained with DAPI (blue).

The quantification of TUNEL-positive cells shows the increase of apoptotic cell number in a BAC dose-dependent manner. BAC at 0.010% and 0.020% significantly increased the number of TUNEL-positive cells at 24h, and at 24h+24h-recovery compared to control. The antiglaucoma eye drops travoprost and latanoprost, which contain higher concentrations of BAC, showed much higher expression of apoptotic cells than did PF-tafluprost or control. * P<0.02 for 0.010% BAC and P <0.01 for other solutions compared to PBS at the same time point; # P <0.005 compared to 0.010% BAC at the same time point; $ P <0.001 compared to 0.020% BAC at the same time point; ♦ P <0.003 compared to 0.020% BAC-latanoprost at the same time point; • P <0.03 compared to 0.015% BAC-travoprost at the same time point

Fig. 3: Inflammation analysis

Immunolocalization of ICAM-1 (CD54) positive cells (green) on 3D-HCE samples after PBS (A), PF-tafluprost (B), 0.005%BAC-bimatoprost (C), 0.010%BAC (D), 0.015%BACtravoprost (E), 0.02%BAC-latanoprost (F) or 0.020% BAC (G) treatments after 24 h of treatment without or with the 24-h post-incubation period. Quantification of ICAM-1 (H) was performed after 24 h of treatment with all solutions followed or not by a 24-h post-incubation period.

ICAM-1-expressing cells were already found in the PBS (A) and PF-tafluprost (B) groups with the same level of expression after 24h or 24h+24h. The increase of ICAM-1 expressions was observed after the treatments with 0.005%BAC-bimatoprost (C) and 0.010%BAC (D) with and without the 24-h post-incubation period. However, few remaining ICAM-1 positive cells were observed after treatments with 0.015%BAC-travoprost (E), 0.020%BAClatanoprost (F) and 0.020% BAC (G) on the 3D-HCE samples due to the erosion. The nuclei were stained with propidium iodide (PI, red). The quantification of TUNEL-positive cells confirmed the observations. * P <0.001 compared to PBS at the same time point; # P<0.005 compared to 0.010% BAC at the same time point; $ P <0.0001 compared to 0.020% BAC at the same time point; ♦ P <0.002 compared to 0.020%BAC-latanoprost at the same time point; • P <0.008 compared to 0.015%BAC-travoprost at the same time point.

Fig. 4: Proliferation analysis

Immunolocalization of Ki-67 positive cells (green) on 3D-HCE samples after PBS (A), PFtafluprost (B), 0.005%BAC-bimatoprost (C), 0.01%BAC (D), 0.015%BAC-travoprost (E), 0.020%BAC-latanoprost (F) or 0.020% BAC (G) treatments after 24h of treatment without (left column) or with the 24h recovery period (right column).

Ki67 immunostaining showed an epithelium with proliferating cells not restricted to the basal layer after PBS treatment (A) at 24h. The same aspect was found after the treatment of PFtafluprost (B). Numerous proliferating cells, with a great number located in all the layers were found after the treatments of 0.005%BAC-bimatoprost (C), 0.010%BAC (D), 0.015%BACtravoprost (E), 0.020%BAC-latanoprost (F) and 0.020% BAC (G) at 24h. At 24h+24hrecovery, we observed the persistence of some proliferative cells after the treatments of PBS or PF-tafluprost. However, at this time, there were no more proliferative cells after the treatments with all BAC-containing solutions. Immunofluorescence analysis of occludin (green: A to G) and ZO-1 (green: H to N) expressions using en-face confocal microscopy after treatment with PBS (A), PF-tafluprost (B), 0.005%BAC-bimatoprost (C), 0.010%BAC (D), 0.015%BAC-travoprost (E), 0.020%BAC-latanoprost (F) and 0.020% BAC (G) treatments after 24 h (left column) or 24h+24h-recovery (right column). The nuclei were stained with DAPI (blue).

After the treatment with PBS or PF-tafluprost, the 3D-HCE disclosed a thin occludin/ZO-1 immunostaining in the most superficial cells of samples by forming a green ring around, leaving a diffuse cytoplasmic staining. We could still find a slight tight junction expression in the apical cells with 0.005%BAC-bimatoprost (C). The expressions disappeared completely after the treatments with 0.010%BAC (D), 0.015%BAC-travoprost (E), 0.020%BAClatanoprost (F) and 0.020% BAC (G) at 24h or 24h+24h-recovery, leaving dense green patches. Note that the nuclei became smaller after the treatments of high dose BACcontaining solutions.

Fig. 5 :

 5 Fig. 5: Tight junction: occludin and ZO-1

2

The study was supported by an unrestricted grant from Santen Oy (Tempere, Finland). The authors have no commercial interest in the drugs mentioned in the manuscript.

Competing Interest: None declared.