D Gandolfo
email: gandolfo@cpt.univ-mrs.fr

M M Rahmatullaev
email: mrahmatullaev@rambler.rurozikovu@yandex.ru

Boundary conditions for translation-invariant Gibbs measures of the Potts model on Cayley trees

Keywords: Mathematics Subject Classifications (2010). 82B26; 60K35 Cayley tree, Potts model, boundary condition, Gibbs measure

We consider translation-invariant splitting Gibbs measures (TISGMs) for the qstate Potts model on a Cayley tree of order two. Recently a full description of the TISGMs was obtained, and it was shown in particular that at sufficiently low temperatures their number is 2 q -1. In this paper for each TISGM µ we explicitly give the set of boundary conditions such that limiting Gibbs measures with respect to these boundary conditions coincide with µ.

1 Introduction.

The analysis of translational invariant splitting Gibbs measures of the q-state Potts model on Cayley trees is based on the classification of translation-invariant boundary laws which are in one-to-one correspondence with the TISGMs. Recall that boundary laws are lengthq vectors which satisfy a non-linear fixed-point equation (tree recursion).

It has been known for a long time that for the anti-ferromagnetic Potts model there exists a unique TISGM [START_REF] Rozikov | On pure phase of the anti-ferromagnetic Potts model on the Cayley tree[END_REF] and for the ferromagnetic Potts model at sufficiently low temperatures there are at least q + 1 translation-invariant Gibbs measures [START_REF] Ganikhodjaev | On pure phases of the three-state ferromagnetic Potts model on the second-order Bethe lattice[END_REF], [START_REF] Ganikhodzhaev | On pure phases of the ferromagnetic Potts model Bethe lattices[END_REF].

One of the q + 1 well-known measures mentioned above is obtained as infinite-volume limit of the finite-dimensional Gibbs measures with free boundary condition and each of the remaining q measures is obtained as the corresponding limit with the boundary conditions of homogeneous (constant) spin-configurations. While the q measures with homogeneous boundary conditions are always extremal in the set of all Gibbs measures [START_REF] Ganikhodjaev | On pure phases of the three-state ferromagnetic Potts model on the second-order Bethe lattice[END_REF], [START_REF] Ganikhodzhaev | On pure phases of the ferromagnetic Potts model Bethe lattices[END_REF], for the free boundary condition measure there is a temperature, denoted by T 0 , which is below the transition temperature, such that the measure is an extremal Gibbs measure if T ≥ T 0 and loses its extremality for even lower temperatures [START_REF] Rozikov | Gibbs measures on Cayley trees[END_REF]Theorem 5.6.].

Recently, in [START_REF] Kuelske | Description of all translation-invariant (splitting) Gibbs measures for the Potts model on a Cayley tree[END_REF] all TISGMs for the Potts model were found on the Cayley tree of order k ≥ 2, and it is shown that at sufficiently low temperatures their number is 2 q -1. We note that the number of TISGMs does not depend on k ≥ 2. In the case k = 2 the explicit formulae for the critical temperatures and all TISGMs are given.

In [START_REF] Kuelske | Fuzzy transformations and extremaity of Gibbs measures for the Potts model on a Cayley tree[END_REF] some regions for the temperature ensuring that a given TISGM is (non-)extreme in the set of all Gibbs measures are found. In particular the existence of a temperature interval is shown for which there are at least 2 q-1 + q extremal TISGMs.

The fact that these measures can never be nontrivial convex combinations of each other (i.e., they are extremal in the set of all TIGMs) is almost automatic (see [START_REF] Kuelske | Description of all translation-invariant (splitting) Gibbs measures for the Potts model on a Cayley tree[END_REF]Theorem 2]). However it is not clear what kind of boundary conditions are needed to get the remaining 2 q -q -2 TISGMs as corresponding limits with the boundary conditions. In this paper we shall answer this question. It is non-trivial problem since the number of TISGMs (i.e., 2 q -1) is larger than the number (i.e., q) of translation-invariant configurations. Therefore one expects to need non-translation-invariant boundary conditions for some TISGMs. Concerning the Ising model, the dependence of TISGMs on boundary conditions has been studied in [START_REF] Higuchi | Remarks on the limiting Gibbs states on a (d+1)-tree[END_REF].

The paper is organized as follows. Section 2 contains preliminaries (necessary definitions and facts). In section 3 we will show how to connect boundary laws with boundary conditions, moreover we shall give the list of known TISGMs. Section 4 contains our main result, namely given any TISGM µ, we show how to compute explicitly a set of boundary conditions such that the limiting Gibbs measures with respect to these boundary conditions coincide with µ. In the last section we construct concrete boundary conditions.

Definitions

Let Γ k = (V, L) be the regular Cayley tree, where each vertex has k + 1 neighbors with V being the set of vertices and L the set of edges.

Two vertices t, s ∈ V, (t = s) are called neighbors if they are connected by an edge. In this case we write s, t . Each vertex of Γ k has k + 1 neighbors.

Fix an origin 0 of Γ k . We write s → t, if t = s and the path connecting 0 and t passes through s. If s → t and s, t are neighbors, then t is called a direct successor of s and this we write as s → 1 t.

For any finite A ⊂ V , the boundary ∂A of A is

∂A = {t ∈ V \ A : ∃x ∈ A, x, t }.
For every A ⊂ V , let Ω A = {1, 2, ..., q} A be the set of all possible spin configurations on A. For brevity we write Ω instead of Ω V .

For every A ⊂ V we define the σ-algebra B A by

B A = the σ -algebra generated by {X t , t ∈ A},
where X t (σ) = σ(t) for all t ∈ A, σ ∈ Ω. For brevity we write B instead of B V .

Let A be a finite subset of V , ω ∈ Ω and σ ∈ Ω A . We define Potts interaction energy on A given the inner configuration σ and the boundary condition ω by

E ω A (σ) = -J t,s∈A: t,s δ σ(t)σ(s) -J t∈A, s∈∂A: t,s δ σ(t)ω(s) , (2.1)
where J ∈ R and δ is the Kronecker's delta.

A finite Gibbs measure P ω A on Ω A corresponding to E ω A is defined by

P ω A (σ) = [Z ω A ] -1 exp[-E ω A (σ)], σ ∈ Ω A , (2.2)
where

Z ω A = σ∈Ω A exp[-E ω A ( σ)].
As usual P ω A can be considered as a probability measure on (Ω, B).

For fixed J, if there is an increasing sequence of finite subsets {V n } such that V n V as n → ∞ and P ω = w -lim n→∞ P ω Vn (the weak convergence of measures) exists for suitable fixed ω ∈ Ω, then P ω is called a limiting Gibbs measure with boundary condition ω for J. On the other hand, a Gibbs measure P for J is defined as a probability measure on (Ω, B) such that for every M in

B A P (M |B A c )(ω) = P ω A (M ). a.s.(P ) (2.3)
It is known ( [START_REF] Georgii | Gibbs measures and phase transitions[END_REF], [START_REF] Higuchi | Remarks on the limiting Gibbs states on a (d+1)-tree[END_REF]) that the set (J) of all Gibbs measures for a fixed J is a nonempty, compact convex set. A limiting Gibbs measure is a Gibbs measure for the same J. Conversely, every extremal point of (J) is a limiting Gibbs measure with a suitable boundary condition for the same J. It is known (see page 241 of [START_REF] Georgii | Gibbs measures and phase transitions[END_REF]) that any extreme Gibbs measure of a Hamiltonian with nearest-neighbor interactions is a splitting Gibbs measure (which is equivalently called a tree-indexed Markov chain [START_REF] Georgii | Gibbs measures and phase transitions[END_REF]). Consequently, any non-splitting Gibbs measure is not extreme. However, any splitting Gibbs measure (not necessary extreme) is a limiting Gibbs measure, because it corresponds to a (generalized)1 boundary condition satisfying a compatibility (tree recursion) condition of Kolmogorov's theorem. In [START_REF] Coquille | Examples of DLR states which are not weak limits of finite volume Gibbs measures with deterministic boundary conditions[END_REF] it was shown that for non-extremal Gibbs measures on Z d a Gibbs measure need not be a limiting Gibbs measure (see [START_REF] Coquille | Examples of DLR states which are not weak limits of finite volume Gibbs measures with deterministic boundary conditions[END_REF] and [START_REF] Friedli | Statistical Mechanics of Lattice Systems: a Concrete Mathematical Introduction[END_REF] for more details).

Translation-invariant limiting Gibbs measures

Let |t| denote the distance between 0 and t ∈ V , i.e. |t| = n if there exists a chain 0

→ 1 u 1 → 1 u 2 → 1 u 3 → 1 ... → 1 u n-1 → 1 t.
We only consider the sequence of boxes

V n = {t ∈ V : |t| ≤ n}, n ≥ 1.
For every s ∈ V we define

Γ k s = {s} ∪ {t ∈ V : s → t}, and V n,s = Γ k s ∩ V n , n ≥ 1.
Elements of the set V 1 \ {0} are the nearest neighbors of the origin 0, since this set contains k + 1 elements, we number them by 1, 2, . . . , k + 1. We note that the subtrees V n,i , i = 1, 2, . . . , k + 1 are similar to each other, i.e., for any i, j ∈ {1, . . . , k + 1} the subtree V n,i can be obtained from the V n,j by a rotation around 0. Moreover we have

V n = {0} ∪ k+1 i=1 V n,i , and V n,i = {i} ∪ j:i→ 1 j V n-1,j . (3.1)
Using (3.1), from (2.1) we get

E ω Vn (σ) = k+1 i=1 E ω V n,i (σ) -Jδ σ(0)σ(i) . (3.2)
For every ω ∈ Ω, s ∈ V \ {0} and n ≥ |s| define

W ω n,s (l) = σ∈Ω Vn,s : σ(s)=l exp[-E ω Vn,s (σ) -Jδ lω(t) ], l = 1, 2, ..., q, (3.3) R l n,s (ω) = W ω n,s (l) W ω n,s (q)
, l = 1, 2, ..., q, (3.4) here t is the unique vertex such that t → 1 s. Now by (2.2) and (3.2)-(3.4) we get

P ω Vn (σ(0) = l) P ω Vn (σ(0) = q) = σ:σ(0)=l exp -E ω Vn (σ) σ:σ(0)=q exp -E ω Vn (σ) = k+1 i=1 (exp(J) -1)R l n,i (ω) + q-1 p=1 R p n,i (ω) + 1 exp(J) + q-1 p=1 R p n,i (ω) , l = 1, 2, ..., q, (3.5)
By (3.1) and (3.3) we obtain

W ω n,s (l) = u:s→ 1 u [(exp(J) -1)W ω n,u (l) + q p=1 W ω n,u (p)], l = 1, 2, ..., q, (3.6)
and for n > m, η ∈ Ω Vm , we get

P ω Vn ({σ(s) = η(s), s ∈ V m }) = exp[-E η V m-1 (η)] s∈∂V m-1 W ω n,s (η(s)) ξ∈Ω Vm exp[-E ξ V m-1 (ξ)] s∈∂V m-1 W ω n,s (ξ(s)) . (3.7)
From the above equalities we obtain the following Lemma 1. Let ω ∈ Ω be given. If there is N > 0 such that R l n,s (ω) converges as n → ∞ for every s ∈ V \ V N and for every l = 1, ..., q, then P ω = w -lim n→∞ P ω Vn exists. For n ≥ 1, p = 1, 2, ...q, i = 1, ..., k + 1 we denote

A n = {t ∈ V : |t| = n}, N (p) n,i (σ) = |{x ∈ A n ∩ V n,i : σ(x) = p}|.
Lemma 2. Let l = 1, . . . , q and ω be a configuration such that2

c l (ω) = s:t→ 1 s δ lω(s) is independent of t ∈ V \ {0}. Then R l n,i (ω) = R l n,j (ω) for any i, j = 1, 2, ..., k + 1. Proof. Since R l n,i (ω) = W ω n,i (l) W ω n,i (q)
, it suffices to prove that W ω n,i (l) = W ω n,j (l) for any i, j = 1, 2, ..., k + 1.

For the Hamiltonian we have

E ω V n,i (σ) = E σ V n-1,i (σ) -J x∈An∩V n,i x→ 1 y δ σ(x)ω(y) . (3.8)
By the condition of lemma 2 we obtain

E ω V n,i (σ) = E σ V n-1,i (σ) -J x∈An∩V n,i c σ(x) (ω) = E σ V n-1,i (σ) -J q p=1 N (p) n,i (σ)c p (ω).
For W ω n,i (l) we have

W ω n,i (l) = σ∈V n,i :σ(i)=l exp[E σ V n-1,i (σ) -J q p=1 N (p) n,i (σ)c p (ω) -Jδ lω(t) ].
Since V n,i is similar to V n,j , for any i, j ∈ {1, . . . , k + 1}, n ≥ 1, there is an one-to-one correspondence γ between sets Ω V n,i and Ω V n,j , which can be obtained by a rotation of the V n,i on the set V n,j . We note that the Potts interaction energy (2.1) is translation-invariant and by the condition of the lemma the quantity c l (ω) also does not depend on vertices of the tree. Therefore, if γ(σ) = ϕ then

E σ V n-1,i (σ) = E γ(σ) V n-1,j (γ(σ)) = E ϕ V n-1,j (ϕ), N (p) n,i (σ) = N (p)
n,j (ϕ), ∀i, j = 1, . . . , k + 1, p = 1, . . . , q. Using these equalities we get

W ω n,i (l) = σ∈Ω V n,i :σ(i)=l exp[E σ V n-1,i (σ) -J q p=1 N (p) n,i (σ)c p (ω) -Jδ lω(t) ] = ϕ∈Ω V n,j :ϕ(j)=l exp[E ϕ V n-1,j (ϕ) -J q p=1 N (p)
n,j (ϕ)c p (ω) -Jδ lω(t) ] = W ω n,j (l).

Thus R l n,i (ω) = R l n,j (ω) for any i, j = 1, 2, ..., k + 1.

From the above proof it follows that R l n,s (ω) depends only on n -|s|, i.e., we have

R l n,s (ω) = R l n-|s|+1 (ω), l = 1, 2, ...q. ( 3.9)
Then from (3.6) we get the following

Y l n (ω) = kF l (Y 1 n-1 (ω), Y 2 n-1 (ω), ..., Y q-1 n-1 (ω)), (3.10)
where l = 1, ..., q -1, n ≥ 2, Y l n (ω) = ln R l n (ω) and F = (F 1 , ..., F q-1 ) with coordinates F l (x 1 , x 2 , ..., x q-1 ) = ln (exp (J) -1) exp (x l ) + q-1 p=1 exp (x p ) + 1 exp (J) + q-1 p=1 exp (x p )

. (3.11) It is clear that if Y i (ω) is a limit point of Y i n (ω), as n → ∞ then by (3.10) we get Y l (ω) = kF l (Y 1 (ω), Y 2 (ω), .
.., Y q-1 (ω)), l = 1, 2, .., q -1.

(3.12)

For convenience we denote θ = exp(J), h l = Y l (ω), l = 1, 2, ...q -1.

Then the system (3.12) becomes h i = k ln (θ -1)e h i + q-1 j=1 e h j + 1 θ + q-1 j=1 e h j , i = 1, . . . , q -1.

(3.13)

In [START_REF] Kuelske | Description of all translation-invariant (splitting) Gibbs measures for the Potts model on a Cayley tree[END_REF] it is proven that to each solution of (3.13) corresponds a unique Gibbs measure which is called a translation-invariant splitting Gibbs measure (TISGM).

In [START_REF] Kuelske | Description of all translation-invariant (splitting) Gibbs measures for the Potts model on a Cayley tree[END_REF] all solutions of the equation (3.13) are given. By these solutions the full set of TISGMs is described. In particular, it is shown that any TISGM of the Potts model corresponds to a solution of the following equation

h = f m (h) ≡ k ln (θ + m -1)e h + q -m me h + q -m -1 + θ , (3.14)
for some m = 1, . . . , q -1.

Denote

θ m = 1 + 2 m(q -m), m = 1, . . . , q -1. (3.15)
It is easy to see that

θ m = θ q-m and θ 1 < θ 2 < • • • < θ q 2 -1 < θ q 2 ≤ q + 1. (3.16) Proposition 1. [7] Let k = 2, J > 0.
1. If θ < θ 1 then there exists a unique TISGM;

2. If θ m < θ < θ m+1 for some m = 1, . . . , q 2 -1 then there are 1 + 2 m s=1 q s TISGMs which correspond to the solutions h i ≡ h i (θ, s) = 2 ln[x i (s, θ)], i = 1, 2 s = 1, . . . , m of (3.14), where

x 1 (s, θ) = θ-1- √ (θ-1) 2 -4s(q-s) 2s , x 2 (s, θ) = θ-1+ √ (θ-1) 2 -4s(q-s) 2s
.

(3.17)

3. If θ q 2 < θ = q + 1 then there are 2 q -1 TISGMs;

4 If θ = q + 1 the number of TISGMs is as follows

2 q-1 , if q is odd 2 q-1 -q-1 q/2 , if q is even; 5. If θ = θ m , m = 1, . . . , q 2
, (θ q 2 = q + 1) then the number of TISGMs is

1 + q m + 2 m-1 s=1 q s .
The number of TISGMs does not depend on k ≥ 2 see [7, Theorem 1]. But for k ≥ 3 explicit formulas for the solutions are not known. Therefore in this paper we consider only the case k = 2.

Following [START_REF] Kuelske | Description of all translation-invariant (splitting) Gibbs measures for the Potts model on a Cayley tree[END_REF] we note that each TISGM corresponds to a solution of (3.14) with some m ≤ q 2 . Moreover, for a given m ≤ q 2 , a fixed solution h i (θ, m) to (3.14) generates q m vectors by permuting coordinates of the vector (h i , h i , . . . , h i m , 0, 0, . . . , 0 q-m ) and giving q m TISGMs. Thus without loss of generality we can only consider the measure µ i (θ, m) corresponding to vector h(m, i) = (h i , h i , . . . , h i m , 0, 0, . . . , 0 q-m-1

), i.e., normalized on qth coordinate (see Remark 2 and Corollary 2 of [START_REF] Kuelske | Description of all translation-invariant (splitting) Gibbs measures for the Potts model on a Cayley tree[END_REF]). Denote by µ 0 ≡ µ 0 (θ) the TISGM corresponding to solution h i ≡ 0 and by µ i ≡ µ i (θ, m) the TISGM corresponding to the solution h i (θ, m), i = 1, 2, m = 1, . . . , q 2 (given in Proposition 1). In this paper our aim is to obtain measures µ i by changing boundary conditions.

Boundary conditions for TISGMs

The following lemma can be proved by simple analysis. Lemma 3.

i. For k ≥ 2 and θ > 1 the function f m (h), h ∈ R defined in (3.14) has the following properties:

a) {h : f m (h) = h} = {0, h 1 , h 2 }, if θ > θ m , m ≤ [q/2]; b) a < f m (h) < A, with a = k ln q-m q+θ-m-1 , A = k ln θ+m-1 m ; c) d dh f m (h) = k(θ-1)(θ+q-1)e h
(me h +θ+q-m-1)((θ+m-1)e h +q-m) > 0;

ii. If k = 2 and m ≤ q/2 then for solutions h 1 and h 2 mentioned in Proposition 1 the following statements hold

0 < h 1 = h 2 , if θ = θ m 0 < h 1 < h 2 , if θ m < θ < θ c , with θ c = q + 1 0 = h 1 < h 2 , if θ = q + 1 h 1 < 0 < h 2 , if q + 1 < θ.
For each solution h i (θ, m) we want to find ω = ω(h i ) ∈ Ω, such that µ i (θ, m) = P ω , where P ω is defined in Lemma 1.

Consider the dynamical system (3.10) for k = 2. Denote G(h) = 2F (h). For a given initial vector v (0) = (v

(0) 1 , . . . , v (0)
q-1 ), we shall study the limit

lim n→∞ G (n) (v (0) ), (4.1)
here

G (n) (v) = G(G(...G(v))...) n .
Figures 1234show the streamlines of the vector field G (n) (v) for k = 2, q = 3. These figures also illustrate the limit points of (4.1).

Denote

I m = {v ∈ R q-1 : v 1 = • • • = v m , v m+1 = • • • = v q-1 = 0}.
It is easy to see that the set I m is invariant with respect to G, i.e. G(I m ) ⊂ I m .

The following lemma gives the limits of (4.1) on the invariant I m (compare with Fig. 1234).

G (n) (v (0) ) =    h(m, 1), if v (0) ∈ I m and v (0) 1 ≥ h 1 (0, . . . , 0), if v (0) ∈ I m and v (0) 1 < h 1 (4.2) 2) If θ m < θ < θ c = q + 1 then lim n→∞ G (n) (v (0) ) =          h(m, 2), if v (0) ∈ I m and v (0) 1 > h 1 h(m, 1), if v (0) ∈ I m and v (0) 1 = h 1 (0, . . . , 0), if v (0) ∈ I m and v (0) 1 < h 1 (4.3) 3) If θ = θ c then lim n→∞ G (n) (v (0) ) =    h(m, 2), if v (0) ∈ I m and v (0) 1 > 0 (0, . . . , 0), if v (0) ∈ I m and v (0) 1 ≤ 0 (4.4) 4) If θ > θ c then lim n→∞ G (n) (v (0) ) =          h(m, 2), if v (0) ∈ I m and v (0) 1 > 0 h(m, 1), if v (0) ∈ I m and v (0) 1 < 0 (0, . . . , 0), if v (0) ∈ I m and v (0) 1 = 0 (4.5)
Proof. Restrict function G(h) to I m , then we get the jth coordinate of G(h) (for any j = 1, . . . , m) is equal to f m (h) which is introduced in (3.14). Other coordinates of G(h) are equal to 0. By Lemma 3 we have that f m is an increasing function. Here we consider the case when the function f m has three fixed points 0, h 1 , h 2 . This proof is more simple for cases when f m has two fixed points. We prove the part 2), other parts are similar. In case 2), by Lemma 3 we have that 0 < h 1 < h 2 and the point h 1 is a repeller, i.e., f m (h 1 ) > 1 and the points 0, h 2 are attractive, i.e., f m (0) < 1, f m (h 2 ) < 1. Now we shall take arbitrary x 0 > 0 and prove that

x n = f m (x n-1 ), n ≥ 1 converges as n → ∞. Consider the following partition (-∞, +∞) = (-∞, 0) ∪ {0} ∪ (0, h 1 ) ∪ {h 1 } ∪ (h 1 , h 2 ) ∪ {h 2 } ∪ (h 2 , +∞). For any x ∈ (-∞, 0) we have x < f m (x) < 0, since f m is an increasing function, from the last inequalities we get x < f m (x) < f 2 m (x) < f m (0) = 0.
Iterating this argument we obtain f n-1 m (x) < f n m (x) < 0, which for any x 0 ∈ (-∞, 0) gives x n-1 < x n < 0, i.e., x n converges and its limit is a fixed point of f m , since f m has unique fixed point 0 in (-∞, 0] we conclude that the limit is 0. For x ∈ (0, h 1 ) we have h 1 > x > f (x) > 0, consequently x n > x n+1 , i.e., x n converges and its limit is again 0. Similarly, one can show that if x 0 > h 1 then x n → h 2 as n → ∞. .

Lemma 4 . 1 )

 41 If θ = θ m , for some m = 1, . . . , q 2 then lim n→∞

Figure 1 : 2 √ 2 .

 122 Figure 1: Here and in Figures2-4the streamlines of the vector field G (n) (v) for k = 2, q = 3 are shown. These figures also illustrate the limit points of (4.1). The plane is formed by a horizontal v 1 -axis, and a vertical v 2 -axis. This figure applies for the caseθ = θ 1 = 1 + 2 √ 2.Four fixed points. The origin is an attractor. There are 3 saddle fixed points.

Figure 2 :

 2 Figure 2: θ = 3.9 > θ 1 . Seven fixed points. Four of them (black dots) are attractors. Three (rectangular dots) points are saddles coming from the saddle points of Fig.1.

Figure 3 :

 3 Figure 3: θ = 4. Four fixed points. The origin is a repeller point. Other fixed points are attractors.

Figure 4 :

 4 Figure 4: θ = 4.5. Seven fixed points. The origin is a repeller, other rectangular dots are saddles. The black dots are attractors.

Figure 5 :

 5 Figure 5: An example of a boundary condition for the TISGM µ 0, q = 5, m = 2. Here c 1 (ω) = c 2 (ω) = 0, c 3 (ω) = c 4 (ω) = c 5 (ω) = 1.Note that the values 1, 2 occur sufficiently often keeping the conditions c 1 (ω) = c 2 (ω) = 0..

Figure 6 :

 6 Figure 6: An example of a boundary condition for the TISGM µ 0 , for q = 15, m = 3. Here we have c i (ω) = 1, i = 1, 2, 3 and c j (ω) = 2, j = 4, 5, . . . , 15..

Adding a boundary field at each site of the boundary is called a generalized boundary condition[START_REF] Gandolfo | On free energies of the Ising model on the Cayley tree[END_REF] or boundary law[START_REF] Georgii | Gibbs measures and phase transitions[END_REF]

The sum in the RHS of c l (ω) is taken over all direct successors of t, it should not be confused with the sum over all neighbors of t, i.e., s: t,s .

This remark and some examples below are added corresponding to a suggestion of a reviewer.

Acknowledgements

U.Rozikov thanks Aix-Marseille University Institute for Advanced Study IMéRA (Marseille, France) for support by a residency scheme. We thank both referees for their helpful suggestions.

By (3.10) the asymptotic behavior of the vector Y n (ω) = (Y 1 n (ω), . . . , Y q-1 n (ω)) depends only on the vector Y 1 (ω), where Y l 1 (ω) = J c l (ω) -c q (ω) , l = 1, . . . , q -1. (4.6)

For a given m ∈ {1, . . . , q 2 } and J > 0 we introduce the following sets of configurations:

. Now taking the coordinates of an initial vector as in (4.6) by Lemma 1 and Lemma 4 we get the following Theorem 1. 1) If θ = θ m , for some m = 1, . . . , q 2 then

2) If θ m < θ < θ c = q + 1 then

3) If θ = θ c then

4) If θ > θ c then

In the next section we use Theorem 1 to construct some concrete boundary conditions.

Construction of boundary conditions

In this section for k = 2, J > 0, q ≥ 3 and m ∈ {1, . . . , q 2 } we shall give examples of boundary configurations.

For k = 2 by Proposition 1 there are up to 2 q -1 TISGMs . We shall consider only µ 0 (θ) corresponding to h = (0, 0, . . . , 0) and µ i (θ, m) corresponding to vector h(m, i) = (h i , h i , . . . , h i m , 0, 0, . . . , 0 q-m-1

), i = 1, 2 with

.

(5.1)

Using Theorem 1 we shall give some boundary conditions for each measure µ i . Boundary conditions for the remaining measures can be obtained by using the permutation symmetry of the Potts model. Case µ 0 . If θ < θ 1 then µ 0 is a unique measure, and one can take any boundary configuration ω to have P ω = µ 0 . But for θ ≥ θ 1 one has to check the conditions of Theorem 1 to have the limiting measure equal to µ 0 .

For example, if θ = θ m < q + 1 for some m ∈ {1, . . . , q 2 } then we must take

Remark 1. For a given TISGM to find its boundary condition one has to construct configurations ω which satisfy the system (like (5.2)) derived by corresponding sufficient conditions of Theorem 1. Below we give several examples of such configurations. It will be clear from our examples that some TISGM may have an infinite set of boundary configurations 3 .

Since h 1 > 0 (see Lemma 3), the system (5.2) is satisfied, for example, if ω satisfies one of the following

one has ω(x 1 ) = i, ω(x 2 ) = i; and c j (ω) = 1 for each j = m + 1, . . . , q, i.e., if ω(x) = j then one has ω(x 1 ) = j but ω(x 2 ) = j. See Fig. 5 for an example of such a configuration for q = 5 and m = 2.

• c i (ω) = 1, i = 1, 2, . . . , m, i.e., if ω(x) = i then ω(x 1 ) = i but ω(x 2 ) = i and it contains j ∈ {m + 1, . . . , q} in such a way that if ω(x) = j then on direct successors x 1 , x 2 of x one has ω(x 1 ) = ω(x 2 ) = j. In this case c i (ω) = 2 for each i = m+1, . . . , q. See Fig. 6 for an example of such a configuration for q = 15 and m = 3.

Case µ 1 . This measure exists for θ ≥ θ m . Subcase θ = θ m . By Theorem 1 for µ 1 we have the condition ω ∈ B + m,1 ∪ B 0 m,1 , i.e.

Note that J = ln θ. Assume ln θ m ≥ h 1 = 2 ln θm-1 2m which is equivalent to the following

(5.4)

Under condition (5.4) the system (5.3) is satisfied for example, if m ≥ 1 and ω is such that c i (ω) = 1, i = 1, . . . , m and c j (ω) = 0, i = m + 1, . . . , q. See Fig. 7 for an example of such a configuration for q = 5 and m = 2 and Fig. 8 for configuration in case q = 10, m = 4. Subcase θ m < θ < q + 1. From Theorem 1 for µ 1 we have the condition

Assume θ is a solution to the equation ln θ = h 1 . Computer analysis shows that this equation has a solution if for example q = 17, m = 1 or q = 55, m = 2. So assuming existence of such a solution θ = θ * one can take a boundary condition configuration as in the previous case (like in Fig. 7) Subcase θ = q + 1. In this case we have µ 1 = µ 0 . Therefore the boundary condition can be taken as in Case µ 0 .

Subcase θ > q + 1. For µ 1 we should have

we can take a configuration ω such that c i (ω) = 0, i = 1, . . . , m and c j (ω) = 1, j = m + 1, . . . , q. (See Fig. 5 for such a configuration).

Remark 2.

From above examples one can see that depending on the temperature (equivalently depending on the parameter θ) a configuration may be the boundary condition for different TISGMs. For example, the configuration given in Fig. 5 is a boundary condition for TISGM µ 0 if θ = θ 2 < q + 1, but the same configuration is the boundary condition for TISGM µ 1 if θ > q + 1.

Case µ 2 . Check the conditions of Theorem 1: Subcase θ = θ m . In this case we have µ 2 = µ 1 , i.e. the boundary condition is constructed in the previous case.

Subcase θ m < θ < q + 1. From Theorem 1 for µ 2 we have the condition

If ln θ > h 1 then it is easy to see that ω satisfies the condition (5.7) if m ≥ 1 and ω is such that c i (ω) = 1, i ∈ {1, . . . , m} and c j (ω) = 0, i ∈ {m + 1, . . . , q} (like in Fig. 7) Subcase θ ≥ q + 1. For µ 2 we should have

(5.8) Condition (5.8) is easily checkable. For example, configurations shown in Fig. 7 and Fig. 8 satisfy this condition.


Daniel Gandolfo
email: gandolfo@cpt.univ-mrs.fr

Roger Rodriguez

Henry C Tuckwell
email: henry.tuckwell@adelaide.edu.au

Mean field analysis of large-scale interacting populations of stochastic conductance-based spiking neurons using the Klimontovich method

Keywords: Computational neuroscience, Conductance-based neural models, Fitzhugh-Nagumo model, Stochastic differential equations, Klimontovich method, Mean-field limits, Neural networks

Mean field analysis of large-scale interacting populations

of stochastic conductance-based spiking neurons using the Klimontovich method

Introduction

A quantitative understanding of brain neuronal activity requires the construction of detailed models of neuronal networks, the development and application of powerful numerical tools and mathematical techniques to analyze them. The large number of cells and their complex anatomy and physiology present challenges which must be met in order to achieve these goals. Thus even though there is much experimental detail available on the electrophysiology and neurochemistry of individual cells, incorporating these data at the macroscopic level is generally intractable. In order to make some progress in this endeavor, it is therefore useful to make several simplifications at the single cell and population levels. We find it useful to employ mathematical methods from the kinetic theory of plasmas [START_REF] Yu | Kinetic Theory of Nonideal Gases and Nonideal Plasmas[END_REF][START_REF] Ichimaru | Statistical Plasma Physics[END_REF][START_REF] Nicholson | Introduction to Plasma Theory[END_REF], whose concepts will be developed here within the framework of brain nervous activity. Such an approach was used recently for deterministic neural networks [START_REF] Buice | Dynamic Finite Size Effects in Spiking Neural Networks[END_REF].

With large populations of neurons the goal is a probabilistic description which predicts, not the behavior of any given neuron in a particular state, but is able to give sufficiently precise estimates of the probability that any individual neuron is in a particular state, or equivalently, the proportion of neurons in a specified state. The probabilistic approach of the network dynamics has a long history in computational neuroscience. It has been found to be essential in the analysis and description of many aspects of brain functioning [START_REF] Tuckwell | Linear Cable Theory and Dendritic Structure[END_REF][START_REF] Tuckwell | Stochastic Processes in the Neurosciences[END_REF][START_REF] Tuckwell | Stochastic Partial Differential Equation Models in Neurobiology: linear and nonlinear models for spiking neurons[END_REF][START_REF] Tuckwell | Stochastic Modeling of Spreading Cortical Depression[END_REF][START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF][START_REF] Rodriguez | Statistical properties of stochastic nonlinear dynamical models of single spiking neurons and neural networks[END_REF], at the levels of ion channel, single neuron and neuronal population. Probabilistic descriptions are also important in studies of the mechanisms of synchronization appearing in many fields of physics as shown by Kuramoto [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF][START_REF] Kuramoto | Cooperative Dynamics of Oscillator Community A Study Based on Lattice of Rings[END_REF][START_REF] Kuramoto | Cooperative Dynamics in Complex Physical Systems[END_REF]. Using this approach, neural populations, which can be viewed as weakly connected oscillators, have often been analyzed [START_REF] Mirollo | The spectrum of the locked state for the Kuramoto model of coupled oscillators[END_REF][START_REF] Abbott | Asynchronous States in Networks of Pulse-Coupled Oscillators[END_REF][START_REF] Daido | Onset of cooperative entrainment in limit-cycle oscillators with uniform all-to-all interactions : bifurcation of the order function[END_REF][START_REF] Hildebrand | Kinetic Theory of Coupled Oscillators[END_REF][START_REF] Buice | Dynamic Finite Size Effects in Spiking Neural Networks[END_REF]. Other neural systems including integrate and fire neuron models, and their generalizations, have also been described within probabilistic frameworks [START_REF] Tuckwell | Cortical network modeling: Analytical methods for firing rates and some properties of networks of LIF neurons[END_REF][START_REF] Ly | Critical Analysis of Dimension Reduction by a Moment Closure Method in a Population Density Approach to Neural Network Modeling[END_REF][START_REF] Nicola | Bifurcations of large networks of two-dimensional integrate and fire neurons[END_REF][START_REF] Nicola | Non-smooth bifurcations of mean field systems of two-dimensional integrate and fire neurons[END_REF][START_REF] Knight | Dynamics of Encoding in Neuron Populations: Some General Mathematical Features[END_REF][START_REF] Apfaltrer | Population density methods for stochastic neurons with realistic synaptic kinetics: Firing rate dynamics and fast computational methods[END_REF][START_REF] Hansel | Existence and stability of persistent states in large neural networks[END_REF].

One of the main features of our approach is to incorporate, at the level of the individual neurons, a large degree of biological plausibility by describing them according to extended Hodgkin-Huxley (HH) frameworks [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF][START_REF] Rodriguez | Noisy spiking neurons and networks: useful approximations for firing probabilities and global behavior description[END_REF][START_REF] Tuckwell | Analytical and simulation results for stochastic Fitzhugh-Nagumo neurons and neural networks[END_REF]. About connectivity between cells, our model incorporates neuroplasticity properties through the introduction of connection variables that are potential dependent. A stochastic dynamical system is built for such cells which are organized in different interacting populations. Then, the concept of set variables which is adapted from Klimontovich's approach [START_REF] Yu | Kinetic Theory of Nonideal Gases and Nonideal Plasmas[END_REF][START_REF] Ichimaru | Statistical Plasma Physics[END_REF][START_REF] Nicholson | Introduction to Plasma Theory[END_REF] of the kinetic theory of plasmas is introduced. Neural population probability distributions (PPD) are defined as expectations of these set variables with respect to the probability distribution of the full interacting system which evolves through a Fokker Planck equation. A hierarchical set of non closed exact coupled equations for the different moments of the set variables can be derived. However, these Fokker-Planck type equations are rather difficult to analyze and numerically simulate. To overcome this difficulty, our approach uses the rather traditional framework of mean field theory. If the number of cells in the network is made arbitrarily large, these cells being distributed in different populations, a set of non local, nonlinear integro-partial differential equations (IPDE) can be derived for the PPDs. It is interesting to observe that the resulting equations we obtain, are of the same type as the McKean-Vlasov-Fokker-Planck (MVFP) equations which have been recently derived [START_REF] Baladron | Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons[END_REF][START_REF] Touboul | Limits and Dynamics of Stochastic Neuronal Networks with Random Heterogeneous Delays[END_REF][START_REF] Mischler | On a kinetic Fitzhugh Nagumo model of neuronal network[END_REF][START_REF] Faugeras | A constructive mean-field analysis of multi-population neural networks with random synaptic weights and stochastic inputs[END_REF] using rigorous techniques of probability theory (see also [START_REF] De Masi | Hydrodynamic Limit for Interacting Neurons[END_REF]). In the case of Fitzhugh-Nagumo (FN) neural systems, we have computed the solution of the PPDs in both cases of globally coupled and time dependent inhomogeneously connected networks.

In section 2, we present the derivation of the IPDE suitable for the analysis of interacting populations of conductance-based spiking model neurons of the (HH) type. In section 3, the special case of networks of Fitzhugh-Nagumo model neurons is considered. Various statistical measures are defined which are suitable for the analysis of the dynamical behavior of these neuronal populations from a direct integration of the stochastic system and from the evaluation of solutions of IPDE. Excitatory and inhibitory uniformly connected networks are analyzed as well as excitatory non uniformly time dependent connected networks. Moreover, the case of two large-scale interacting Fitzhugh-Nagumo populations is considered. Section 4 describes the methods and parameters used in the simulations presented in section 3. Conclusion follows in section 5.

2 Integro-partial differential equations for population dynamical systems of coupled large-scale conductance-based neural networks.

We consider networks of conductance-based model neurons with non-uniform synaptic connectivity. These networks are organized into a set of populations with interactions both within and between them. Let us call P γ , γ = 1, 2, . . . , P a typical population, which is composed of K γ cells.

2.1 The dynamical system for a given population P γ

The state variables of each cell i = 1, 2, . . . , K γ in P γ are its membrane potential V γ i and R γ i which is an m-dimensional set of auxiliary variables. The dynamical laws for these auxiliary variables are of the Hodgkin Huxley type

d R γ i dt = Ψ γ (V γ i , R γ i ). (1)
More precisely, the j th component of the m-dimensional vector R γ i evolves according to laws of the form

d dt ( R γ i ) j = Ψ γ (V γ i , R γ i ) j = A 1 j (V γ i )(1 -( R γ i ) j ) + A 2 j (V γ i )( R γ i ) j (2)
where A 1 j , A 2 j : R → R, j = 1, 2, . . . , m are nonlinear functions of the membrane potential. We call I γ (V γ i , R γ i ) the set of ionic, active channels and passive leak, currents across the membrane of the i th cell. External currents I γ i,ext (t) may also be applied on all cells of the network. They are composed of a deterministic part I γ i,ext (t) and a stochastic part η γ i,t ,

I γ i,ext (t) = I γ i,ext (t) + η γ i,t (3)
where η γ i,t is a white noise such that

η γ i,s η γ j,t = δ ij δ(s -t) β γ i,V , i, j = 1, 2, . . . , K γ , (4)
where δ ij is the Kronecker symbol, δ(•) is the delta distribution and β γ i,V i = 1, 2, . . . , K γ are noise parameters. As in [START_REF] Baladron | Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons[END_REF], we assume that the deterministic part of the current and the noise parameters are the same for all neurons in each given population, I γ i,ext (t) = I γ ext (t) and β γ i,V = β γ V , i = 1, 2, . . . , K γ . Let us now introduce the connectivity properties of the model. Highly nonuniform connectivity has been observed in many biological neural networks. For example, in [START_REF] Song | Highly Nonrandom Features of Synaptic Connectivity in Local Cortical Circuits[END_REF], using statistical techniques, structural features of synaptic connectivity have been shown to be quite different for local cortical circuitry from those of random networks. In [START_REF] Mclaughlin | Large-scale modeling of the primary visual cortex: influence of cortical architecture upon neuronal response[END_REF], a computational model has been developed for the primary visual cortex where emphasis is given to the role of cortical architecture, particularly through synaptic connectivity, in the functional cortical activity.

Our approach to the variability of the connection strengths between cells is made in terms of new variables entering in the description of synaptic currents. More precisely, let us call I γ syn,i (t) the current entering the i th cell as a result of the firing activity of all other cells in the network,

I γ syn,i (t) = 1 K γ Kγ j=1 M (V γ i , R γ i , φ γ i , V γ j , R γ j , φ γ j ) (5)
where the bounded coupling functions M will be given explicitly below.

The connection coefficients between pairs (i, j) are scaled by the factor 1 Kγ , so any divergence of this interaction term is prevented when the network size increases (see e.g. [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF], [START_REF] Baladron | Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons[END_REF]). In our approach, the synaptic current I γ syn,i (t) depends also on real functions φ γ i (t) and φ γ j (t) which characterize the time dependent non homogeneous coupling J γ ij (t) between cells i and j through a suitable function Γ γ , namely

J γ ij (t) = Γ γ (φ γ i (t), φ γ j (t)). (6)
Neurons in neocortical and other networks in the brain usually have specific sets of afferent inputs with a variety of strengths (see for example [START_REF] Mclaughlin | Large-scale modeling of the primary visual cortex: influence of cortical architecture upon neuronal response[END_REF]). Moreover, it is assumed that there is a dynamical law for these φ i variables of the form

dφ γ i dt = Ω γ (φ γ i , V γ i , R γ i ) (7)
for a suitable function Ω γ . In what follows, the functions Γ γ and Ω γ will be chosen in such a way to ensure positivity of couplings and the bounded variation of φ-variables. In that way, the variables φ γ i (resp φ γ j ) which build the connection matrix J γ ij (t) as time evolves, depend on postsynaptic (resp. presynaptic) activity. In the following, we shall consider homogeneous networks for which these φ variables are constant and inhomogeneous networks for which are considered specific simple expressions for the functions Γ γ and Ω γ . In both cases, the forms of I γ and I γ syn,i will be specified. Finally, let us summarize the structure of the network dynamical system (omitting the capacitance parameter). For i = 1, 2, . . . , K γ , the system can be written

dV γ i dt = I γ (V γ i , R γ i ) + I γ ext (t) + I γ syn,i (t) + η γ i,t d R γ i dt = Ψ γ (V γ i , R γ i ) (8)
dφ γ i dt = Ω γ (φ γ i , V γ i , R γ i )
which can be put in the following vector form

dZ γ i dt = F γ (Z γ i ) + 1 K γ Kγ j=1 M γ (Z γ i , Z γ j ) + ζ γ i,t (9)
where the vectors

Z γ i , F γ (Z γ i ), ζ γ i,t , M γ (Z γ i , Z γ j ) in R m+2 are given by Z γ i = (V γ i , R γ i , φ γ i ) (10)
F γ (Z γ i ) = (I γ (V γ i , R γ i ) + I γ ext (t), Ψ γ (V γ i , R γ i ), Ω γ (φ γ i , V γ i , R γ i )) (11)
M γ (Z γ i , Z γ j ) = (M γ (V γ i , R γ i , φ γ i , V γ j , R γ j , φ γ j ), 0 , 0 ) (12)
ζ γ i,t = (η γ i,t , 0 , 0 ) (13)
0 being the null vector in R m . Note that in [START_REF] Kuramoto | Cooperative Dynamics in Complex Physical Systems[END_REF], M γ includes the coupling term Γ γ (φ γ i (t), φ γ j (t)) (see [START_REF] Tuckwell | Stochastic Partial Differential Equation Models in Neurobiology: linear and nonlinear models for spiking neurons[END_REF]) accounting for the amplitude of the synaptic conductance involved in the connection between the presynaptic cell j and postsynaptic cell i.

The structure established in [START_REF] Rodriguez | Statistical properties of stochastic nonlinear dynamical models of single spiking neurons and neural networks[END_REF] is now generalized to a set of coupled neural populations.

The mean field approach for coupled conductance-based neural populations

One considers non uniformly connected HH type neurons which are organized in a set of P populations. The dynamical system governing the state of this set, at each time t, assumes the following form

dZ γ i dt = F γ (Z γ i ) + ζ γ i,t + P α=1 1 K α Kα j=1 M γα (Z γ i , Z α j ), i = 1, 2, . . . , K γ , γ = 1, 2, . . . , P. (14)
We now derive mean field population equations for [START_REF] Sakaguchi | Cooperative Phenomena in Coupled Oscillator Systems under External Fields[END_REF].

Let us denote by

p t ((Z α i ) i,α ) = p t ((Z 1 i ) i=1,2,...,K1 , . . . , (Z P j ) j=1,2,••• ,K P ) (15)
the joint probability distribution of the stochastic variables

(Z α i ) i,α = (Z α i )
α=1,2,...,P i=1,2,...,Kα and introduce the notation

∂ ∂Z (H(Z)) = m+2 j=1 ∂ ∂Z j (H(Z) j )
where Z = (Z j ) j=1,2,...,m+2 ∈ R m+2 and H(Z) = (H(Z) j ) j=1,2,...,m+2 ∈ R m+2 .

The Fokker Planck equation for the system (14) is

∂ ∂t p t ((Z α i ) i,α ) = - P γ=1 Kγ i=1 ∂ ∂Z γ i (F γ (Z γ i ) p t ((Z α i ) i,α )) (16)
- P γ=1 Kγ i=1 ∂ ∂Z γ i P α=1 1 K α Kα j=1 M γα (Z γ i , Z α j ) p t ((Z α i ) i,α ) + 1 2 P γ=1 Kγ i=1 (β γ V ) 2 ∂ 2 ∂(V γ i ) 2 p t ((Z α i ) i,α )
We now adopt the Klimontovich approach [START_REF] Ichimaru | Statistical Plasma Physics[END_REF] which has been successfully developed for the kinetic theory of gases and plasmas. The method is here adapted to the derivation of a probabilistic description of the system (14) of noisy interacting spiking neural populations (see also [START_REF] Buice | Dynamic Finite Size Effects in Spiking Neural Networks[END_REF]). One defines the following set variable

n µ (U ) = 1 K µ Kµ i=1 δ(Z µ i -U ) (17)
where δ(•) is the Dirac distribution and Z µ i is the solution of (14) written for the population P µ . Z µ i and U (see [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF]) are in R m+2 . In what follows, the following notation is used: U = (u, ẑ, s), u, s in R, ẑ in R m . The expectation value of the stochastic variables n µ (U ) with respect to the probability distribution p t is denoted by n µ (U, t), so that

n µ (U, t) = n µ (U ) pt . ( 18
)
We call n µ (U, t) the neural population probability distribution (PPD) for the population P µ . We now derive an equation for n µ (U, t). The time derivative

∂ ∂t n µ (U, t) is composed of 3 terms ∂ ∂t n µ (U, t) = γ µ 1 + γ µ 2 + γ µ 3 . (19)
Let us consider each of them separately. The first one, γ µ 1 is given by

γ µ 1 = - R m+2 P δ=1 K δ l=1 dZ δ l P γ=1 Kγ i=1 ∂ ∂Z γ i {F γ (Z γ i )p t ((Z α i ) i,α )} 1 K µ Kµ j=1 δ(Z µ j -U ). (20)
Being a probability distribution, p t has nice vanishing properties for sufficiently large values of the variables Z µ i . Thus, a simple integration by parts on [START_REF] Buice | Dynamic Finite Size Effects in Spiking Neural Networks[END_REF] enables us to deduce the following expression for γ µ

1 γ µ 1 = R m+2 P δ=1 K δ l=1 dZ δ l P γ=1 Kγ i=1 {F γ (Z γ i )p t ((Z α i ) i,α )} ∂ ∂Z γ i 1 K µ Kµ j=1 δ(Z µ j -U ). (21)
Clearly,

∂ ∂Z γ i Kµ j=1 δ(Z µ j -U ) = 0 if γ = µ, so that ∂ ∂Z µ i Kµ j=1 δ(Z µ j -U ) = ∂ ∂Z µ i δ(Z µ i -U ), i = 1, 2, . . . , K µ .
Hence ( 21) can be rewritten

γ µ 1 = -∂ ∂U    R m+2 P δ=1 K δ l=1 dZ δ l Kµ i=1 {F µ (U )p t ((Z α i ) i,α )} 1 K µ δ(Z µ i -U )    (22) = -∂ ∂U R m+2 P δ=1 K δ l=1 dZ δ l {F µ (U ) n µ (U )p t ((Z α i ) i,α )} . (23)
Finally, using [START_REF] Hildebrand | Kinetic Theory of Coupled Oscillators[END_REF], γ µ 1 is given by

γ µ 1 = -∂ ∂U (F µ (U )n µ (U, t)). (24)
Let us now consider the second term γ µ 2 in ( 19)

γ µ 2 = - R m+2 P δ=1 K δ l=1 dZ δ l P γ=1 Kγ i=1 ∂ ∂Z γ i P α=1 1 K α Kα l=1 M γα (Z γ i , Z α l )p t ((Z α i ) i,α ) 1 K µ Kµ j=1 δ(Z µ j -U ). (25)
By the same argument used for the derivation of ( 24), an integration by parts in [START_REF] Ly | Critical Analysis of Dimension Reduction by a Moment Closure Method in a Population Density Approach to Neural Network Modeling[END_REF] implies

γ µ 2 = R m+2 P δ=1 K δ l=1 dZ δ l P γ=1 Kγ i=1 P α=1 1 K α Kα l=1 M γα (Z γ i , Z α l )p t ((Z α i ) i,α ) ∂ ∂Z γ i 1 K µ Kµ j=1 δ(Z µ j -U ) (26)
= - ∂ ∂U R m+2 P δ=1 K δ l=1 dZ δ l Kµ i=1 P α=1 1 K α Kα l=1 M µα (U, Z α l )p t ((Z α i ) i,α ) 1 K µ δ(Z µ i -U ) = - ∂ ∂U R m+2 dU R m+2 P δ=1 K δ l=1 dZ δ l Kµ i=1 P α=1 1 K α 1 K µ Kα l=1 M µα (U, U ) p t ((Z α i ) i,α )δ(Z µ i -U )δ(Z α l -U ) = - ∂ ∂U R m+2 dU R m+2 P δ=1 K δ l=1 dZ δ l P α=1 M µα (U, U )p t ((Z α i ) i,α ) n µ (U ) n α (U )
which again leads to

γ µ 2 = - ∂ ∂U R m+2 dU P α=1 M µα (U, U ) n µ (U ) n α (U ) pt . (27)
The term γ µ 2 has been obtained without the use of approximations and hence is exact. However, it is not really useful in applications. A way to go further consists in making the so called mean field estimate (see e.g. [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF])

n µ (U ) n α (U ) pt ≈ n µ (U ) pt n α (U ) pt . ( 28
)
This approximation is valid because the fluctuations of

n µ (U ) (resp. n α (U )) are small for K µ (resp. K α ) large and are O( 1 √ Kµ ) (resp. O( 1 √ Kα )
). Accordingly, the coupling term γ µ 2 takes the form :

γ µ 2 = - ∂ ∂U R m+2 dU P α=1 M µα (U, U )n µ (U, t)n α (U , t). (29)
The last diffusive term γ µ 3 can be derived in a similar fashion. It is given by

γ µ 3 = 1 2 R m+2 P δ=1 K δ l=1 dZ δ l P γ=1 Kγ i=1 (β γ V ) 2 ∂ 2 ∂(V γ i ) 2 p t ((Z α i ) i,α ) 1 K µ Kµ j=1 δ(Z µ j -U ). (30)
Integration by parts gives

γ µ 3 = 1 2 R m+2 P δ=1 K δ l=1 dZ δ l P γ=1 Kγ i=1 (β γ V ) 2 p t ((Z α i ) i,α ) 1 K µ ∂ 2 ∂(V γ i ) 2 Kµ j=1 δ(Z µ j -U ). (31)
Here also, one has

∂ 2 ∂(V γ i ) 2 Kµ j=1 δ(Z µ j -U ) = 0 if γ = µ (32) while ∂ 2 ∂(V µ i ) 2 Kµ j=1 δ(Z µ j -U ) = ∂ 2 ∂(V µ i ) 2 δ(Z µ i -U ). (33)
Thus

γ µ 3 = 1 2 R m+2 P δ=1 K δ l=1 dZ δ l p t ((Z α i ) i,α ) 1 K µ Kµ i=1 (β µ V ) 2 ∂ 2 ∂(V µ i ) 2 δ(Z µ i -U ). (34)
So

γ µ 3 = (β µ V ) 2 2 ∂ 2 ∂u 2 R m+2 P δ=1 K δ l=1 dZ δ l p t ((Z α i ) i,α ) 1 K µ Kµ i=1 δ(Z µ i -U ) (35)
and hence

γ µ 3 = (β µ V ) 2 2 ∂ 2 ∂u 2 n µ (U, t). (36)
Finally, the neural population probability distributions n µ (U, t) which are the expectation values of the set variables n µ (U ) µ = 1, 2, . . . , P satisfy the following system of (nonlinear)

IPDE ∂ ∂t n µ (U, t) = - ∂ ∂U (F µ (U )n µ (U, t)) - ∂ ∂U R m+2 dU P α=1 M µα (U, U )n µ (U, t)n α (U , t) + 1 2 (β µ V ) 2 ∂ 2 ∂U 2 n µ (U, t) (37)
µ = 1, 2, . . . , P, U = (u, z, s), u, s ∈ R, z ∈ R m , U = (u , z , s ), u , s ∈ R, z ∈ R m .
Let us recall U = (U l ) l=1,2,...,m+2 . We write explicitly the components of U as U = (u, z, s), u, s ∈ R, z ∈ R m , where u (resp. z, s) has potential (resp. recovery, synaptic) meaning.

We have thus obtained a system of equations for the study of interacting large-scale neural populations using a mean field approach. The dynamical behavior of the neurons in these populations is developed in terms of HH concepts which are usually capable of giving a quantitative description of ionic currents across the neural membranes. Interneuronal connections in these networks are not necessarily constant over time. A synaptic plasticity mechanism is proposed (see ( 5)-( 7)). For such systems, equation [START_REF] Faugeras | A constructive mean-field analysis of multi-population neural networks with random synaptic weights and stochastic inputs[END_REF], which describes the evolution over time of the PPD for the population P µ which interacts with populations P ν , ν = 1, 2, . . . , P is of the same type as that obtained in [START_REF] Baladron | Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons[END_REF] where it was shown, using methods of probability theory, that in the limit of large populations, the network may display the property of propagation of chaos. In such a way, all neurons which belong to a given population, asymptotically have the same probability distribution which characterizes a mean field stochastic process. The equation [START_REF] Faugeras | A constructive mean-field analysis of multi-population neural networks with random synaptic weights and stochastic inputs[END_REF] is the (MVFP) equation for this process.

3 Application to coupled large-scale Fitzhugh-Nagumo networks

Population probability distributions for coupled Fitzhugh-Nagumo networks

In what follows, we give an illustration of the general result [START_REF] Faugeras | A constructive mean-field analysis of multi-population neural networks with random synaptic weights and stochastic inputs[END_REF] in the case of FN model neurons which leads to a simple example of the stochastic system [START_REF] Ermentrout | Mathematical Foundations of Neuroscience[END_REF]. The basic excitability properties of many neurons are approximately described by the FN model and the mathematical tractability of the resultant system makes it suitable as a first application of the theory.

The dynamical system is now given by (for all i = 1, 2, . . . , K γ , γ = 1, 2, . . . , P )

dV γ i dt = F γ (V γ i , X γ i ) + I γ ext (t) + η γ i,t + P α=1 1 K α Kα j=1 M γα (V γ i , X γ i , φ γ i , V α j , X α j , φ α j ) (38)
dX γ i dt = G γ (V γ i , X γ i ) (39)
dφ γ i dt = Ω γ (φ γ i , V γ i , X γ i ) (40)
where

F γ (V, X) = -k γ V (V -a γ )(V -1) -X (41) G γ (V, X) = b γ (V -m γ X). (42)
The variables V γ i represent the membrane potential of the i th cell in P γ while, in FN networks, the variables Rγ i which have been introduced in (1) are called X γ i and are one-dimensional. The parameters k γ , a γ , b γ , m γ govern the dynamics of the FN neural model in the population P γ (see e.g. [START_REF] Tuckwell | Linear Cable Theory and Dendritic Structure[END_REF] for typical values of these parameters). The quantities I γ ext (t), η γ i,t have the same meaning as in section 2.1. In the following sections, we shall give a more precise form for the coupling term M γα (V γ i , X γ i , φ γ i , V α j , X α j , φ α j ) in the case of uniformly and non uniformly distributed sets of connection parameters between the cells of these populations. According to the general result (37), in the mean field limit, the PPDs n µ (V, X, φ, t) µ = 1, 2, . . . , P , which are built for a set of interacting populations of FN neurons, satisfy the following system of (nonlinear) integro-partial differential equations

∂ ∂t n µ (V, X, φ, t) = (43)
- ∂ ∂V ((F µ (V, X) + I µ ext (t))n µ (V, X, φ, t)) - ∂ ∂X (G µ (V, X)n µ (V, X, φ, t)) - ∂ ∂φ (Ω µ (φ, V, X)n µ (V, X, φ, t)) - ∂ ∂V R 3 dV dX d φ P α=1 M µα (V, X, φ, V , X , φ )n µ (V, X, φ, t)n α (V , X , φ , t) + 1 2 (β µ V ) 2 ∂ 2 ∂V 2 n µ (V, X, φ, t) µ = 1, 2, . . . , P.
We shall give examples of applications of this mean-field approach to neural networks dynamics by considering the cases of firstly one population and then two populations of FN cells. Equation (43) will be solved numerically and its solutions compared to the numerical solutions of the stochastic system (38)-(40).

Statistical measures from the stochastic and integro-partial differential systems

Our goal is to have, at least numerically, a reasonably good understanding of the dynamics of the above considered set of neural populations. From the solution of the system of IPDE (43), it is in principle possible to have insights into this dynamics, in the limiting case of infinitely many cells. Since the real situation corresponds generally to a large (but finite) number of cells, the solutions of ( 43) may be viewed only as an approximation.

Another approach for the study of the network dynamical behavior is based on a direct numerical evaluation of solutions of the (stochastic) coupled differential equations ( 38)- [START_REF] Achard | Perinatal maturation of the respiratory rhythm generator in mammals: from experimental results to computational simulation[END_REF]. Let us call this system SCDE. In [START_REF] Brette | Simulation of networks of spiking neurons: a review of tools and strategies[END_REF], one can find a review of various methods of numerical solution of these systems, and in [START_REF] Djurfeldt | Large-scale modeling-a tool for conquering the complexity of the brain[END_REF] are given appropriate techniques for solving such very large systems.

Clearly, the solution of IPDE has a great advantage over the direct solutions of the K γ dimensional SCDE, γ = 1, 2, . . . , P since only P equations must be solved.

Both systems (IPDE and SCDE) have been numerically implemented over the same domain

Ω × Λ × Φ, where Ω = [V min , V max ] (resp. Λ = [X min , X max ], Φ = [φ min , φ max ]
) are the bounded domains of variation of potential (resp. recovery, synaptic coupling variables).

Moreover, similar initial conditions have been chosen in both cases. More precisely, a normalized Gaussian function n(V, X, φ, t = 0) has been considered for IPDE with mean and standard deviation parameters (V 0 , σ V ), (X 0 , σ X ), (φ 0 , σ φ ). For SCDE, all cells initial data (V γ i (t = 0), X γ i (t = 0), φ γ i (t = 0)), i = 1, 2, . . . , K γ , γ = 1, 2, . . . , P were selected through the use of i.i.d. Gaussian random variables with the same moments as the ones used for the solution of IPDE.

Concerning the numerical aspects of the problem, a discretization procedure is necessary for IPDE, requiring a partition of the parameter space Ω × Λ × Φ. Obviously, this method of solution takes much less computing time than the possibly large numbers of trials required for SCDE.

In order to compare the results obtained by the solutions of IPDE and SCDE, it is necessary to show how one can extract information on each dynamical variable V, X or φ, for any cell in the full system, at any time, from a probabilistic point of view. If one concentrates on the potential variable, for example, a neural population potential distribution has to be defined.

Concerning IPDE, we have at our disposal the PPD n µ (V, X, φ, t), solution of ( 43) and this potential distribution can be defined as the marginal

ρ µ (V, t) = R 2 n µ (V, X, φ, t) dXdφ. (44)
In that way, any statistical measure may be obtained through ρ µ (V, t). In the following, we call V (t) ρ µ , Var (V (t)) ρ µ the mean and variance of the membrane potential with respect to ρ µ . We define a firing measure ϕ ρ µ (t), as the probability that, at a given time t, the potential V (t) > θ where θ is some prescribed threshold

ϕ ρ µ (t) = +∞ θ ρ µ (V , t)dV . ( 45
)
Regarding SCDE, we propose a numerical measure of the PPD in the following way. Rather than following the trajectories of individual neurons (each starting from initial conditions (V, X, φ) µ i (t = 0), i = 1, 2, . . . , K µ , µ = 1, 2, . . . , P ), we set up a counting process about the presence (at a given time t, for any cell i in P µ ) of the potential V µ i (t) in a given subinterval [a j , b j ] = bin j . The intervals [a j , b j ], j = 1, 2, . . . , N build a partition of Ω, all having the same length λ = (V max -V min )/N (see section 4 for domain extensions).

Because of the presence of noise, one has to perform sufficiently many trials. Accordingly, the PPD for SCDE will be defined as the discrete distribution ξ µ (t) given by

ξ µ (t) = { p µ j (t), j = 1, 2, . . . , N } = N (V µ i (t) ∈ [a j , b j ]) K µ * N trials , 1 < i < K µ , j = 1, 2, . . . , N (46)
N denotes the number of occurrences for all trials and N trials the total number of trials which have been used in the simulations. Clearly, ξ µ (t) is a piecewise constant function, each p µ j (t) giving the average number of events in each bin j , j = 1, 2, . . . , N . Therefore, for SCDE, the mean and variance of V (t) with respect to this discrete probability are given by

V (t) ξ µ = N j=1 (V min + λj) p µ j (t) (47)
Var (V (t)) ξ µ = N j=1 (V min + λj) 2 p µ j (t) - N j=1 (V min + λj) p µ j (t) 2 . (48)
Finally, we evaluate a discrete firing measure ϕ ξ µ (t) as the probability that the potential

V (t) > θ ϕ ξ µ (t) = j > (θ-Vmin )/λ p µ j (t). (49)
In the next section, we compare the two types of probability distributions, the continuous one ρ µ (V, t) and the discrete one ξ µ (t), as well as averages, standard deviations and firing measures derived from these two distributions.

Numerical results for a single population

Let us first give the form of the equation [START_REF] Lindsey | Computational Models and Emergent Properties of Respiratory Neural Networks[END_REF] in the case of only one population, namely P = 1. In the system ( 38)- [START_REF] Achard | Perinatal maturation of the respiratory rhythm generator in mammals: from experimental results to computational simulation[END_REF], the matrix

M γα (V γ i , X γ i , φ γ i , V α j , X α j , φ α j ) reduces to the term M 11 (V 1 i , X 1 i , φ 1 i , V 1 j , X 1 j , φ 1 j )), i = 1, 2, . . . , K 1 . The triples (V 1 i , X 1 i , φ 1 i ) are rewritten (V i , X i , φ i ) and M 11 , K 1 , n 1 (V, X, φ, t
) are rewritten M , K and n(V, X, φ, t), respectively. Similarly, the densities ρ 1 (V, t) and ξ 1 (t) will be denoted by ρ(V, t) and ξ(t).

Excitatory synapses, uniform connectivity

In the present section, we assume that the population has all to all connectivity as in [START_REF] Buice | Dynamic Finite Size Effects in Spiking Neural Networks[END_REF][START_REF] Fasoli | A formalism for evaluating analytically the cross-correlation structure of a firing-rate network model[END_REF]. For discussions on this topic see [START_REF] Amit | Dynamics of a recurrent network of spiking neurons before and following learning[END_REF][START_REF] Fasoli | The complexity of dynamics in small neural circuits[END_REF][START_REF] Tuckwell | Cortical network modeling: Analytical methods for firing rates and some properties of networks of LIF neurons[END_REF].

We give here an explicit form of the synaptic current I syn,i (t) introduced in (5), in the case where all cells in the network are connected with excitatory synapses. Because of the homogeneous character of the coupling considered here, the variables φ i , i = 1, 2, . . . , K assume constant values. We call J E the constant coupling term Γ(φ i , φ j ) (see [START_REF] Tuckwell | Stochastic Partial Differential Equation Models in Neurobiology: linear and nonlinear models for spiking neurons[END_REF]) and V E S the corresponding parameters for resting values. Moreover, a simple form has been considered for the detection of presynaptic events using sigmoidal functions σ E Θ , defined by

σ E Θ (V ) = (1 + e -β(V -V E Θ ) ) -1 (50)
for some prescribed threshold values V E Θ , and parameter β. One has

M (V i , X i , V j , X j ) = J E (V E S -V i )σ E Θ (V j ) (51)
where the indices i (resp. j) refer to post (resp. pre) synaptic cells. Finally, the synaptic input current I syn,i (t) for each cell i is given by

I syn,i (t) = 1 K K j=1 M (V i , X i , V j , X j ). (52)
Accordingly, the stochastic system (SCDE) ( 38)-( 40) takes the form

dV i dt = F (V i , X i ) + I ext (t) + η i,t + I syn,i (t) dX i dt = G(V i , X i ) i = 1, 2, . . . , K (53)
while the integro-partial differential equation (IPDE) [START_REF] Lindsey | Computational Models and Emergent Properties of Respiratory Neural Networks[END_REF] for the PPD n(V, X, t) is

∂ ∂t n(V, X, t) = -∂ ∂V {(F (V, X) + I ext (t)) n(V, X, t)} - ∂ ∂X {G(V, X) n(V, X, t)} (54)
-

J E ∂ ∂V (V E S -V ) n(V, X, t) R 2 dV dX σ E Θ (V ) n(V , X , t) + β 2 V 2 ∂ 2 ∂V 2 n(V, X, t).
Numerical integration of (53), (54) have been performed, with suitable initial conditions. The results are shown in Fig. 1 below. A 3-dimensional graph of the time behavior of the discrete probability density ξ(t) is represented in Fig. 1(a) which may be compared to the behavior of the continuous one ρ(V, t) shown in Fig. 1

(b).

Similarly, in Fig. 2(a) and Fig. 2(b), the mean values and variances of the potential variables obtained from the solutions of SCDE and IPDE, are compared. These quantities give a global description of the network activity. It is seen that there is close agreement between the two methods of solution.

In addition, one can obtain insight into the activity of cells from a local point of view. This is illustrated in Fig. 2(c) showing the superimposed time development, over several trials, of the potential of some arbitrarily chosen cells in the network. There appears significant variability in the activity of individual cells in the network. However, as shown in Fig. 2(a), an order emerges in the mean behavior due to synchronizations between cells.

Excitatory synapses, time dependent non uniform connectivity

In section 3.3.1, we considered the case of one population with uniform and constant connections over time. In this section, we present the structure of the mean field equation and the nature of the numerically obtained solutions, in the case of one population (i.e. P = 1) whose non-uniform interconnections can vary over time.

In section 2.1, the functions Γ and Ω have been introduced for the definition of the coupling terms. As in section 3.3.1, the dynamical behavior of synaptic excitatory conductances is viewed in terms of the sigmoidal dependence σ E Θ (V ) with respect to presynaptic depolarization. The connection term M (see (52)), which we consider now, is M

(V i , X i , φ i , V j , X j , φ j ) = (V E S -V i ) J E Γ(φ i , φ j ) σ E Θ (V j )
where V E S and J E are constant. Thus, the postsynaptic current which is applied to each cell i = 1, 2, . . . , K is Finally, the structure of the network dynamical system in the case of non uniform, time dependent couplings between cells, which we wish to analyze from the kinetic point of view, is the following

I syn,i (t) = (V E S -V i ) J E K K j=1 Γ(φ i , φ j )σ E Θ (V j ). (55)
dV i dt = F (V i , X i ) + I ext (t) + η i,t + I syn,i (t) dX i dt = G(V i , X i ) dφ i dt = Ω(φ i , X i ). (56)
The general result [START_REF] Faugeras | A constructive mean-field analysis of multi-population neural networks with random synaptic weights and stochastic inputs[END_REF] can now be applied to produce the structure of the mean field equation which can be derived for this kind of neural network system. One looks for the probability distribution n(V, X, φ, t) of depolarization, recovery and connection variables (V, X, φ) ∈ R 3 at time t. From (37), n(V, X, φ, t) must satisfy

∂ ∂t n(V, X, φ, t) = -∂ ∂V {(F (V, X) + I ext (t)) n(V, X, φ, t)} - ∂ ∂X {G(V, X) n(V, X, φ, t)} (57) -J E ∂ ∂V (V E S -V ) n(V, X, φ, t) R 3 dV dX d φ σ E Θ (V ) Γ(φ, φ ) n(V , X , φ , t) - ∂ ∂φ {Ω (φ, X) n(V, X, φ, t)} + β 2 V 2 ∂ 2 ∂V 2 n(V, X, φ, t).
Our approach here is the same as in previous sections where results obtained for SCDE and IPDE have been compared. The new variables φ i , i = 1, 2, . . . , K have random character, in the same way as potential V i and recovery variables X i . As has been described for potential variables in section 3.2, statistical measures for variables φ i can be defined in both cases. For IPDE, the probability distribution of φ variables is given by the

χ(φ, t) = R 2 n(V, X, φ, t)dV dX . ( 58
)
Thus, for example, the (continuous) mean of the variable φ is obtained from

φ(t) χ = R φ χ(φ , t) d φ . ( 59
)
Concerning SCDE, the discrete version for the probability distribution of φ-variables and related moments is obtained following the same lines as those developed for membrane potential variables (see formulas (46)-( 49)). For all simulations which are described in this section, the function Ω has been taken to be of the form

Ω (φ, X) = -α φ φ + β φ σ φ (X) (60)
where α φ and β φ are parameters and σ φ is a sigmoidal function

σ φ (X) = (1 + e -ν(X-X thresh ) ) -1 (61)
for some prescribed threshold value X thresh and parameter ν. The function Γ appearing in (55) has been taken to be of the Hebbian form

Γ(φ i , φ j ) = φ i φ j , (62)
In Figs. 3(a) and 3(b) are shown the mean values V (t) ξ , φ(t) ξ according to SCDE and V (t) ρ , φ(t) χ according to IPDE, in the absence and in the presence of external current, for a network with excitatory connections. When there is no applied input current, V (t) ξ and V (t) ρ decrease after a sufficiently long time. If the coupling parameter takes increasing values, the network activity is increasingly synchronized, even in the absence of a stimulus. The proposed mechanism (56), together with (55) and (62) leads to significantly increased values of these couplings, especially when an external current is applied to the network during a short period. This is illustrated in Fig. 3 important spiking synchronization occurs that causes, according to the relaxation process in (56), a net increase of φ i variables and their average value. This synchronization is maintained periodically (see also [START_REF] Achard | Perinatal maturation of the respiratory rhythm generator in mammals: from experimental results to computational simulation[END_REF][START_REF] Harris-Warrick | General Principles of Rhythmogenesis in Central Pattern Networks[END_REF][START_REF] Dayan | Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems[END_REF][START_REF] Lindsey | Computational Models and Emergent Properties of Respiratory Neural Networks[END_REF]).

As shown in Fig. 3(c), this approach is part of a widely accepted rule of reinforcement of (maximum values of ) synaptic conductance for a given connection. This strengthening is linked to a causal relationship between spiking activity of pre and postsynaptic neurons in this connection, see [START_REF] Dayan | Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems[END_REF]. In this figure, the temporal behavior of the membrane activities V pre (t), V post (t) and strength variables W pre,post (t) = J E K Γ(φ pre (t), φ post (t)) are represented in a superimposed manner for two arbitrarily chosen cells in the network, for a small number of trials. Before the start of the stimulation, during the background activity, action potentials (APs) which are due to noise may be emitted in an uncoordinated way. The production of these APs does not cause a significant change in connection weights. During and after stimulation, and because of increased synchronization of preand postsynaptic activities, φ pre (t), φ post (t) variables grow significantly and remain at a higher value.

On the other hand, for a network with inhibitory connections, the same mechanisms (56) with ( 60) and (62), do not produce such sharp synchronization as in the previous case where the connections are only excitatory. Outside the period of stimulation, during which the φ i variables are growing, they decrease rapidly causing loss of rhythmicity, see Fig. 3(d).

Two interacting Fitzhugh-Nagumo populations

One cell connected to a network

In section 2.2, a method was presented for the mean field derivation of the dynamical activity of a set of large-scale interconnected neuronal populations. The general form of the interaction terms that allows one to deduce the system of integro-partial differential equations that must give the neuronal probability densities was determined in [START_REF] Nicola | Non-smooth bifurcations of mean field systems of two-dimensional integrate and fire neurons[END_REF]. It may be useful to analyze this relationship in a simple case where some populations consist not of a large number of cells but, in contrast, are composed of only a small number, and even consist of only one cell. The goal here is to show that a mean field derivation may be obtained in the case of neural systems where only one cell is submitted to external stimulation, this cell being connected to a large network consisting of cells whose parameters may be different from those of the stimulated cell. In the notations of section 2.2, Z 1 1 is the dynamical variable of Cell 1 on which is applied an external stimulation, {Z 2 i } i=1,2,...,K 2 denote the dynamical variables of cells which form a large-scale network P 2 , all of them being connected to Cell 1 . The stochastic differential systems which are attached to such neuronal assemblies are

dZ 1 1 dt = F 1 (Z 1 1 ) + ζ 1 1,t + M 11 (Z 1 1 , Z 1 1 ) + 1 K 2 K 2 j=1 M 12 (Z 1 1 , Z 2 j ) (63)
dZ 2 i dt = F 2 (Z 2 i ) + ζ 2 i,t + M 21 (Z 2 i , Z 1 1 ) + 1 K 2 K 2 j=1 M 22 (Z 2 i , Z 2 j ), i = 1, 2, . . . , K 2 . ( 64
)
In our approach, the P 2 population is not submitted to an external stimulation, the term M 11 (Z 1 1 , Z 1 1 ) has the meaning of a self-interaction for Cell 1 , M 21 (Z 2 i , Z are noise terms, the functions F j , j = 1, 2 are given in Section 2.1. For such systems, the set variable of Cell 1 and P 2 are respectively n

1 (U ) = δ(Z 1 1 -U ) and n 2 (U ) = 1 K 2 K 2 i=1 δ(Z 2 i -U ). The terms γ µ k , k = 1, 2, 3, µ = 1, 2 (see section 2.
2) were obtained in a general form. Let us examine more specifically the terms γ µ 2 , as expressed in [START_REF] Nicola | Non-smooth bifurcations of mean field systems of two-dimensional integrate and fire neurons[END_REF], which have been obtained without any approximation. For µ = 1, the Cell 1 term is :

γ 1 2 = - ∂ ∂U R m+2 dU M 11 (U, U ) n 1 (U ) n 1 (U ) pt - ∂ ∂U R m+2 dU M 12 (U, U ) n 1 (U ) n 2 (U ) pt (65)
whereas for µ = 2, the P 2 term is :

γ 2 2 = - ∂ ∂U R m+2 dU M 21 (U, U ) n 2 (U ) n 1 (U ) pt - ∂ ∂U R m+2 dU M 22 (U, U ) n 2 (U ) n 2 (U ) pt . (66)
Since P 1 consists of only one cell, it is clear that in general,

n 1 (U ) n 1 (U ) pt = n 1 (U ) pt n 1 (U ) pt .
However, this term only appears if a self-interaction is considered on Cell 1 . Examination of the other terms in (65), (66) leads to the same conclusions as those that were deduced in section 2.2. Thus, if P 2 is large-scale, admitting fluctuations in n 2 (U ) of small amplitude, n 2 (U )

n 2 (U ) pt and n 1 (U ) n 2 (U ) pt n 1 (U ) pt n 2 (U ) pt . Similarly, n 2 (U ) n 2 (U ) pt n 2 (U ) pt n 2 (U ) pt .
Finally, when one considers a system consisting of a cell Cell 1 and a population P 2 , the structure of mean field integro-partial differential equations is the same as that which has been obtained previously except for the terms of self-interaction of Cell 1 which must be evaluated without approximation.

In the numerical illustration of "fill-in" which is shown in Fig. 4, Cell 1 has no self-interaction and for clarity of the presentation, the backward coupling terms are set to zero. The parameter values used in these simulations are given in section 4.

Two large-scale excitatory and inhibitory connected populations

In this section, we consider two networks P E and P I of neural populations with different parameters k α , a α , b α , m α , α = E, I. The connection matrix function Ψ = M EE M EI M IE M II , which describes the connections inside each population and between the two populations, has its elements M αγ (V α i , X α i , φ α i , V γ j , X γ j , φ γ j ), α, γ = E, I, i = 1, 2, . . . , K E , j = 1, 2, . . . , K I where K E and K I are the cell numbers of each population. As it has been shown in section 2, these matrix elements may have a uniform character or a non uniform one which incorporates pre and postsynaptic aspects. In this section, we have considered the simplest form for these parameters, the main objective here being to show how the integro-partial differential equations system can help in the description and control of different coupled neural systems which are built of a great number of cells. More precisely, all connections inside P E and P I are supposed uniformly distributed. As it has been considered in section 3.3.1, the variables φ α i , i = 1, 2, . . . , K α , α = E, I, assume constant values and synaptic currents between presynaptic cell j and postsynaptic cell i in P E (resp. P I ), are described by the term

M EE (V E i , X E i , V E j , X E j ) = J E (V E S -V E i ) σ E Θ (V E j ) (resp. M II (V I i , X I i , V I j , X I j ) = J I (V I S -V I i ) σ I Θ (V I j )
). The connectivity between presynaptic cell j in P I (resp. P E ) and postsynaptic cell i in P E (resp. P I ) is described by

M EI (V E i , X E i , V I j , X I j ) = J EI (V EI S -V E i ) σ EI Θ (V I j ) (resp. M IE (V I i , X I i , V E j , X E j ) = J IE (V IE S - V I i ) σ IE Θ (V E j )).
As in previous sections, the numerical work is concentrated on a comparison between solutions of noisy dynamical systems, taking sufficiently many cells, and solutions which can be obtained from partial differential equations. The structure of the latter is made of a coupled system of equations for n E (V, X, t) and n I (V, X, t) which are the probability distributions for the excitatory and the inhibitory populations

∂ ∂t n E (V, X, t) = -∂ ∂V (F E (V, X) + I E ext (t))n E (V, X, t) - ∂ ∂X G E (V, X)n E (V, X, t) (67)
-

J E ∂ ∂V (V E S -V )n E (V, X, t) R 2 dV dX σ E Θ (V )n E (V , X , t) -J EI ∂ ∂V (V EI S -V )n E (V, X, t) R 2 dV dX σ EI Θ (V )n I (V , X , t) + 1 2 (β E V ) 2 ∂ 2 ∂V 2 n E (V, X, t)
For n I (V, X, t) the equation is of the same type by making an exchange between the indices I and E. An example of numerical results is shown in Fig. 5.

• Firstly, we have considered the time development of the probability densities in the case where the populations P E and P I are uncoupled (J EI = J IE = 0). An external input I E ext (t) (resp. I I ext (t)) has been applied on P E (resp. P I ). The response of P E , which develops after applying the stimulus, appears as spiking activity in a temporal window (Fig. 5(a)) where the firing probability is large -see Fig. 5(e). In P I , after application of the stimulus, the activity returns to an equilibrium (Fig. 5(b)), in which the firing probability is small -see Fig. 5(f).

• When the coupling between P E and P I is turned on (J EI = 0, J IE = 0) and the same forms of stimuli are applied to the 2 populations, the response of P E , after application of the stimulus, does not change much as seen in Fig. 5(c) and, similarly for the firing probability (Fig. 5(g)). However, in the response of P I , we can see a significant difference from the uncoupled case. Indeed, the response in P I is far from being a simple return to equilibrium. Significant spiking activity appears well after applying the stimulus (Fig. 5(d)), which can be measured precisely in terms of the firing probability -see Fig. 5(h). Here also, we note good agreement between the solutions of SCDE and IPDE.

Numerics and parameter settings

In this section, methods and parameters used in the various simulations of sections 3.3 and 3.4 are presented.

The numerical solutions of systems SCDE and IPDE were sought on the bounded domain

D = Ω × Γ × Φ, such that Ω = [V min , V max ] = [-1.0, 1.8], Γ = [X min , X max ] = [-0.4, 0.6], Φ = [φ min , φ max ] = [-1.0, 3.0]
. The domain amplitudes have been chosen sufficiently large so that, in the numerical simulations, during the temporal evolution, the probability density and its derivatives are negligible outside D.

• SCDE. An Euler method was used for the solution of systems SCDE with a time step δt = 0.01. The required number of cells in each population was K = 200. In each case N = 500 Monte Carlo simulations of the network equations were performed. The same Gaussian distribution was chosen for the random initial conditions of these systems. In the case of homogeneous and constant connections, the mean and dispersion parameters of this distribution for potential and recovery variables were set to V 0 = 0, σ V = 0.15, X 0 = 0, σ X = 0.15.

When dealing with the inhomogeneous case (section 3.3.2), the corresponding Gaussian parameters for the synaptic connections variables were set to φ 0 = 1.0, σ φ = 0.01.

For the determination of the discrete probability densities for the potential and synaptic connections variables (see (46) for ξ µ (t)), a partitioning of both Ω and Φ domains in N = 200 subintervals (bins) was performed.

In our network model with excitatory synapses (section 3.3.1), detections of presynaptic activity causing the appearance of incoming synaptic currents in each cell was achieved by means of sigmoidal functions (see ( 50)) with β = 20, V E Θ = 0.5. Other parameters were introduced to describe these currents, J E = 0.1, V E S = 0.8. In simulations involving 2 different populations (section 3. while internal coupling parameters for each population are J E = 0.25, J I = 0.2. Regarding the couplings between the two populations, one has

V EI S = 1.0, V IE S = -1.0, J EI = J IE = 0.2.
Single populations of FN systems are expressed in terms of the following parameters a = 0.1, γ = 0.2, k = 1.0, b = 0.015, which are defined in ( 41), ( 42) both for homogeneous (section 3.3.1) or inhomogeneous (section 3.3.2) connections. In the case of two coupled homogeneous populations (section 3.4.2), these parameters have the following values: for P E , a E = 0.1, γ E = 0.2, k E = 1.0, b E = 0.015, and for P I , a I = 0.15, γ I = 0.18, k I = 0.9, b I = 0.017.

In section 3.3.2 we introduced a time-dependent inhomogeneous interneuronal connection model. For both excitatory and inhibitory networks, the dynamics of connection variables φ i , i = 1, 2, . . . , K = 200 depends on the parameters α φ = 0.003, β φ = 0.01 and on a sigmoidal function σ φ (X) = (1 + e -ν(X-X thresh ) ) -1 , where ν = 20.0 and X S = 0.1. In both cases J = 0.1. In the excitatory case, V E S = 0.8 while in the inhibitory case, V I S = -0.2. In section 3.4.1 a cell (Cell 1 ) is submitted to an external current I(t). Internal parameters for Cell 1 are a 1 = 0.1, γ 1 = 0.2, k 1 = 1.0, b 1 = 0.015 . Cell 1 is connected to a network P 2 of cells whose parameters are a 2 = 0.15, γ 2 = 0.18, k 2 = 0.9, b 2 = 0.017. P 2 cells are not subjected to external stimulation. The coupling is of forward type from Cell 1 to P 2 with J 12 = 0.02. The coupling from P 2 to Cell 1 is set equal to zero. Internal connections are of excitatory type, settings for these connections are identical to those considered in the previous sections. All these systems undergo the action of white noise (see ( 4)) for which the (diffusion) parameter is set to β V = 0.044 (β 2 V /2 0.001). • IPDE. Solving the nonlinear IPDE of the (MVFP) type for n(V, X, φ, t) requires the knowledge of an initial condition n(V, X, φ, t = 0). We took a Gaussian with the same distribution parameters as the distribution of initial conditions used to solve the stochastic systems. In the general case of inhomogeneous connections, one has

n(V, X, φ, t = 0) = 1 (2π) 3/2 σ V σ X σ φ exp - (V -V 0 ) 2 2σ 2 V - (X -X 0 ) 2 2σ 2 X - (φ -φ 0 ) 2 2σ 2 φ . (68)
Dirichlet boundary conditions were imposed on the bounded domain D. For computer simulation, variable dimension V (resp. X, φ) of the 3-D domain D has been discretized in N 1 (resp. N 2 , N 3 ) points, with

N 1 = N 2 = N 3 = 150.
We have used an explicit midpoint Euler (order 2) algorithm with time step δt = 0.001 which has been found sufficient to ensure stability of the numerical solution of the partial differential equation. During the simulation, regular monitoring of the normalization condition of the solution was made.

Finally, the detection threshold θ of action potentials for the firing probabilities of both SCDE and IPDE (see [START_REF] Djurfeldt | Large-scale modeling-a tool for conquering the complexity of the brain[END_REF] and (49)) was set at θ = 0.8.

Numerical simulations were performed using double precision floating point arithmetic on a Intel Xeon computer with 32 processors.

Conclusion

In recent work on the study of large-scale populations of neurons (see e.g. [START_REF] Baladron | Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons[END_REF]), mathematically rigorous methods from probability theory were developed. Mean field equations have been derived for a set of probability distributions governing the dynamical behavior of a set of noisy populations of neurons. These so called McKean-Vlasov-Fokker-Planck equations appear to bring a particularly interesting perspective in the field of neural networks.

We have developed a method which falls within the same framework of description of the macroscopic behavior of sets of conductance based neuronal populations in terms of the microscopic characteristics of the constituent nerve cells. This method uses some tools of kinetic theory of gases and plasmas [START_REF] Ichimaru | Statistical Plasma Physics[END_REF]. Here, the dynamics on the movements of the particles is replaced by the dynamics of the membrane depolarizations while the synaptic connections substitute for collisions between particles.

This approach of the neural behavior has been the subject of several studies in which cells are seen as coupled oscillators (see e.g. [START_REF] Buice | Dynamic Finite Size Effects in Spiking Neural Networks[END_REF]). This requires making simplifying assumptions (basic variables reduced to phases, weak coupling). The aim of this work is to propose a generalization of these methods to biologically plausible systems which are distributed in several interacting populations.

To achieve this goal, we have introduced stochastic set variables, as does Klimontovich [START_REF] Ichimaru | Statistical Plasma Physics[END_REF], for cells belonging to all populations of the neural ensemble. Fokker Planck techniques thus lead to derivation of a hierarchical set of non closed non local integro-partial differential equations for the probability distributions of the neuronal dynamical variables.

However, these equations remain intractable and require, as in the kinetic theory of plasmas, suitable approximations for further analysis. Developments in statistical moments can be considered. We limited ourselves to a mean field approach that allows closure of the resulting hierarchy and is based on conventional estimates of the smallness of the fluctuations of set variables considered in the case of large neuronal populations.

Models of neurons that we have considered are of Hodgkin Huxley type, whose synaptic connection strengths are not necessarily constant over time, the proposed plasticity model being Hebbian. For such neuronal assemblies of interacting populations, we deduce McKean-Vlasov-Fokker-Planck type systems of equations similar to those obtained in [START_REF] Baladron | Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons[END_REF].

Numerical integration of these equations has been made in the case of Fitzhugh-Nagumo neural networks and showed that statistical measurements obtained in this scheme are in good agreement with those obtained by direct simulations of finite size stochastic neural dynamic systems.

Kinetic theories have been developed in the framework of neural populations of oscillators, allowing analytical development for the study of statistical mechanics concepts, such as phase transitions and bifurcation analysis [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF][START_REF] Sakaguchi | Cooperative Phenomena in Coupled Oscillator Systems under External Fields[END_REF][START_REF] Strogatz | From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators[END_REF]. The theory presented in this paper is a generalization for populations of conductance-based neurons.

While it seems impossible to consider such developments starting from the initial neural stochastic systems, it is hoped that the MVFP equations which have been obtained here, could bring new responses related to these statistical mechanics aspects.

( a )

 a Evolution of the normalized discrete probability density ξ(t), solution of the stochastic dynamical system (51)-(53). (b) Evolution of the normalized continuous probability density ρ(V, t), solution of the MVFP equation (54).

Figure 1 :

 1 Figure 1: Time behavior of probability densities. Initial conditions and simulation parameters, see section 4.

( a )

 a Mean values of potential obtained from SCDE and IPDE. (b) Variance of potential obtained from SCDE and IPDE. (c) Depolarizations versus time for 5 randomly selected cells (SCDE).

Figure 2 :

 2 Figure 2: Excitatory synapses, uniform connectivity. Parameters for the FN system and noise, see section 4.

  marginal

  (b) where we see the emergence of a rhythmic response of the network in the form of trains associated with changes in the average values of the connections. During the stimulation, an (a) Excitatory synapses. Mean values of potential and connection variables obtained from IPDE and SCDE. No applied current. (b) Excitatory synapses. Mean values of potential and connection variables obtained from IPDE and SCDE. A step current is applied: rhythm generation. (c) Excitatory synapses. Superimposed time development, over a sample of 5 trials, of potential and strength variables of 2 arbitrarily chosen cells. (d) Inhibitory synapses. Superimposed time development, over a sample of 5 trials, of potential and strength variables of 2 arbitrarily chosen cells .

Figure 3 :

 3 Figure 3: Non uniform, time dependent connectivity. In Figs.3(c) and 3(d), the strength variables are defined as J i1i2 = Γ(φ i1 , φ i2 ) (see (62)) for 2 arbitrarily chosen cells, Cell i1 and Cell i2 , i 1 = i 2 ∈ {1...K}. For parameters of FN system and noise, see section 4.

  (a) Response of Cell 1 to an external stimulation. (b) Response of P 2 without coupling with Cell 1 . (c) Response of P 2 with coupling with Cell 1 .

Figure 4 :

 4 Figure 4: Excitatory uniform connectivity in P 2 with no external input. External input on P 1 = Cell 1 . No backward connectivity from P 2 to Cell 1 . The stimulation of a single cell leads to a complex response of the P 2 network which was analyzed by both methods.

  4.2), these parameters are V E S = 0.8, V I S = -0.2, No coupling between P E and P I (a) Mean values of potential for P E . No coupling between P E and P I (b) Mean values of potential for P I . With coupling between P E and P I (c) Mean values of potential for P E . With coupling between P E and P I (d) Mean values of potential for P I .No coupling between P E and P I (e) Firing measures for P E , see (45) and (49)No coupling between P E and P I (f) Firing measures for P I , see (45) and (49)With coupling between P E and P I (g) Firing measures for P E , see (45) and (49).With coupling between P E and P I (h) Firing measures for P I , see[START_REF] Djurfeldt | Large-scale modeling-a tool for conquering the complexity of the brain[END_REF] and (49).

Figure 5 :

 5 Figure 5: Behavior of two coupled FN populations, one excitatory (P E ) and one inhibitory (P I ).

  1 1 ) refers to couplings of Cell 1 to P 2 , M 12 (Z 1 1 , Z 2 j ) represents the backward couplings from P 2 to Cell 1 and M 22 (Z 2 i , Z 2 j ) are internal interactions within P 2 . ζ 1 1,t and ζ 2

	i,t