Olivier Faugeras

Jonathan Touboul

Bruno Cessac

A constructive mean field analysis of multi population neural networks with random synaptic weights and stochastic inputs

Keywords: Mean field analysis, stochastic processes, stochastic differential equations, stochastic networks, stochastic functional equations, random connectivities, multi populations networks, neural mass models Analyse champs moyen, processus stochastiques, équations différentielles stochastiques, réseaux stochastiques, équations stochastiques fonctionnelles, connectivités aléatoires, réseaux multi populations, modèles de masses neuronales

come

A constructive mean field analysis of multi population neural networks with random synaptic weights and stochastic inputs

Introduction

Modeling neural activity at scales integrating the effect of thousands of neurons is of central importance for several reasons. First, most imaging techniques are not able to measure individual neuron activity ("microscopic" scale), but are instead measuring mesoscopic effects resulting from the activity of several hundreds to several hundreds of thousands of neurons. Second, anatomical data recorded in the cortex reveal the existence of structures, such as the cortical columns, with a diameter of about 50µm to 1mm, containing of the order of one hundred to one hundred thousand neurons belonging to a few different species. These columns have specific functions. For example, in the visual cortex V1, they respond to preferential orientations of bar-shaped visual stimuli. In this case, information processing does not occur at the scale of individual neurons but rather corresponds to an activity integrating the collective dynamics of many interacting neurons and resulting in a mesoscopic signal. The description of this collective dynamics requires models which are different from individual neurons models. In particular, if the accurate description of one neuron requires "m" parameters (such as sodium, potassium, calcium conductances, membrane capacitance, etc...), it is not necessarily true that an accurate mesoscopic description of an assembly of N neurons requires N m parameters. Indeed, when N is large enough averaging effects appear, and the collective dynamics is well described by an effective mean field, summarizing the effect of the interactions of a neuron with the other neurons, and depending on a few effective control parameters. This vision, inherited from statistical physics requires that the space scale be large enough to include a large number of microscopic components (here neurons) and small enough so that the region considered is homogeneous. This is in effect the case of cortical columns.

However, obtaining the equations of evolution of the effective mean field from microscopic dynamics is far from being evident. In simple physical models this can be achieved via the law of large numbers and the central limit theorem, provided that time correlations decrease sufficiently fast. This type of approach has been generalized to such fields as quantum field theory or non equilibrium statistical mechanics. To the best of our knowledge, the idea of applying mean field methods to neural networks dates back to Amari [START_REF] Amari | Characteristics of random nets of analog neuron-like elements[END_REF][START_REF] Amari | A mathematical foundation for statistical neurodynamics[END_REF]. In his approach, the author uses an assumption that he called the "local chaos hypothesis", reminiscent of Boltzmann's "molecular chaos hypothesis", that postulates the vanishing of individual correlations between neurons, when the number N of neurons tends to infinity. Later on, Crisanti, Sompolinsky and coworkers [START_REF] Sompolinsky | Chaos in Random Neural Networks[END_REF]) used a dynamic mean field approach to conjecture the existence of chaos in an homogeneous neural network with random independent synaptic weights. This approach was formerly developed by Sompolinsky and coworkers for spin-glasses (Crisanti and Sompolinsky, 1987a,b;[START_REF] Sompolinsky | Relaxational dynamics of the Edwards-Anderson model and the mean-field theory of spin-glasses[END_REF], where complex effects such as aging or coexistence of a diverging number of metastable states, renders the mean field analysis delicate in the long time limit [START_REF] Houghton | Role of initial conditions in the mean-field theory of spinglass dynamics[END_REF].

On the opposite, these effects do not appear in the neural network considered in [START_REF] Sompolinsky | Chaos in Random Neural Networks[END_REF] because the synaptic weights are independent [START_REF] Cessac | Increase in complexity in random neural networks[END_REF] (and especially non symmetric, in opposition to spin glasses). In this case, the Amari approach and the dynamic mean field approach lead to the same mean field equations. Later on, the mean field equations derived by Sompolinsky and Zippelius [START_REF] Sompolinsky | Relaxational dynamics of the Edwards-Anderson model and the mean-field theory of spin-glasses[END_REF] for spin-glasses were rigorously obtained by Ben Arous and Guionnet (Ben-Arous andGuionnet, 1995, 1997;[START_REF] Guionnet | Averaged and quenched propagation of chaos for spin glass dynamics[END_REF]. The application of their method to a discrete time version of the neural network considered in [START_REF] Sompolinsky | Chaos in Random Neural Networks[END_REF] and in [START_REF] Molgedey | Supressing chaos in neural networks by noise[END_REF] was done by [START_REF] Moynot | Large deviations and mean-field theory for asymmetric random recurrent neural networks[END_REF].

Mean field methods are often used in the neural network community but there are only a few rigorous results using the dynamic mean field method. The main advantage of dynamic mean field techniques is that they allow one to consider neural networks where synaptic weights are random (and independent). The mean field approach allows one to state general and generic results about the dynamics as a function of the statistical parameters controlling the probability distribution of the synaptic weights [START_REF] Samuelides | Random recurrent neural networks[END_REF]. It does not only provide the evolution of the mean activity of the network but, because it is an equation on the law of the mean field, it also provides informations on the fluctuations around the mean and their correlations. These correlations are of crucial importance as revealed in the paper by Sompolinsky and coworkers [START_REF] Sompolinsky | Chaos in Random Neural Networks[END_REF]. Indeed, in their work, the analysis of correlations allows them to discriminate between two distinct regimes: a dynamics with a stable fixed point and a chaotic dynamics, while the mean is identically zero in the two regimes.

However, this approach has also several drawbacks explaining why it is so seldom used. First, this method uses a generating function approach that requires heavy computations and some "art" for obtaining the mean field equations. Second, it is hard to generalize to models including several populations. Their approach consists in considering that dynamic mean field equations characterize in fine a stationary process. It is then natural to search for stationary solutions. This considerably simplifies the dynamic mean field equations by reducing them to a set of differential equations (see section 5) but the price to pay is the unavoidable occurrence in the equations of a non free parameter, the initial condition, that can only be characterized through the investigation of the non stationary case. Hence it is not clear whether such a stationary solution exists, and, if it is the case, how to characterize it. To the best of our knowledge, this difficult question has only been investigated for neural networks in one paper by Crisanti and coworkers [START_REF] Crisanti | chaos in neural networks : chaotic solutions Dayan P, Abbott LF[END_REF].

Different alternative approaches have been used to get a mean field description of a given neural network and to find its solutions. In the neuroscience community, a static mean field study of multi population network activity was developed by Treves in [START_REF] Treves | Mean-field analysis of neuronal spike dynamics[END_REF]. This author did not consider external inputs but incorporated dynamical synaptic currents and adaptation effects. His analysis was completed in [START_REF] Abbott | Asynchronous states in networks of pulse-coupled neuron[END_REF], where the authors considered a unique population of nonlinear oscillators subject to a noisy input current. They proved, using a stationary Fokker-Planck formalism, the stability of an asynchronous state in the network. Later on, Gerstner in [START_REF] Gerstner | Time structure of the activity in neural network models[END_REF] built a new approach to characterize the mean field dynamics for the Spike Response Model, via the introduction of suitable kernels propagating the collective activity of a neural population in time.

Brunel and Hakim considered a network composed of integrate-and-fire neurons connected with constant synaptic weights [START_REF] Brunel | Fast global oscillations in networks of integrate-and-fire neurons with low firing rates[END_REF]. In the case of sparse connectivity, stationarity, and considering a regime where individual neurons emit spikes at low rate, they were able to study analytically the dynamics of the network and to show that the network exhibited a sharp transition between a stationary regime and a regime of fast collective oscillations weakly synchronized. Their INRIA approach was based on a perturbative analysis of the Fokker-Planck equation. A similar formalism was used in [START_REF] Mattia | Population dynamics of interacting spiking neurons[END_REF] which, when complemented with self-consistency equations, resulted in the dynamical description of the mean field equations of the network, and was extended to a multi population network.

In the present paper, we investigate this question using a new and rigorous approach based on stochastic analysis.

The article is organized as follows. In section 2 we derive from first principles the equations relating the membrane potential of each of a set of neurons as function of the external injected current and noise and of the shapes and intensities of the postsynaptic potentials in the case where these shapes depend only on the post-synaptic neuron (the so-called voltage-based model) and in the case where they depend only on the nature of the presynaptic neurons (the so-called activity-based model). Assuming that the shapes of the postsynaptic potentials can be described by linear (possibly timedependent) differential equations we express the dynamics of the neurons as a set of stochastic differential equations and give sufficient conditions for the equivalence of the voltage-and activity based descriptions. This allows us to obtain the mean field equations when the neurons belong to P populations whose sizes grow to infinity and the intensities of the postsynaptic potentials are independent Gaussian random variables whose law depend on the populations of the pre-and post-synaptic neurons and not on the individual neurons themselves. These equations can be derived in several ways, either heuristically as in the work of Amari [START_REF] Amari | Characteristics of random nets of analog neuron-like elements[END_REF][START_REF] Amari | A mathematical foundation for statistical neurodynamics[END_REF], Sompolinsky [START_REF] Crisanti | chaos in neural networks : chaotic solutions Dayan P, Abbott LF[END_REF][START_REF] Sompolinsky | Chaos in Random Neural Networks[END_REF], and Cessac [START_REF] Cessac | Increase in complexity in random neural networks[END_REF][START_REF] Samuelides | Random recurrent neural networks[END_REF], or rigorously as in the work of Benarous andGuionnet (Ben-Arous andGuionnet, 1995, 1997;[START_REF] Guionnet | Averaged and quenched propagation of chaos for spin glass dynamics[END_REF]. Our purpose in this article is not their derivation but to prove that they are wellposed and to provide an algorithm for computing their solution. Before we do this we provide the reader with two important examples of such mean field equations. The first example is what we call the simple model, a straightforward generalization of the case studied by Amari and Sompolinsky.

The second example is a neuronal assembly model, or neural mass model, as introduced by Freeman [START_REF] Freeman | Mass action in the nervous system[END_REF] and examplified in Jansen and Rit's cortical column model [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF].

In section 3 we consider the problem of solutions over a finite time interval [t 0 , T ]. We prove, under some mild assumptions, the existence and uniqueness of a solution of the dynamic mean field equations given an initial condition at time t 0 . The proof consists in showing that a nonlinear equation defined on the set of multidimensional Gaussian random processes defined on [t 0 , T ] has a fixed point. We extend this proof in section 4 to the case of stationary solutions over the time interval [-∞, T ] for the simple model. Both proofs are constructive and provide an algorithm for computing numerically the solutions of the mean field equations.

We then study in section 5 the complexity and the convergence rate of this algorithm and put it to good use: We first compare our numerical results to the theoretical results of Sompolinsky and coworkers [START_REF] Crisanti | chaos in neural networks : chaotic solutions Dayan P, Abbott LF[END_REF][START_REF] Sompolinsky | Chaos in Random Neural Networks[END_REF]. We then provide an example of numerical experiments in the case of two populations of neurons where the role of the mean field fluctuations is emphasized.

Along the paper we introduce several constants. To help the reader we have collected in table 1 the most important ones and the place where they are defined in the text.

Mean field equations for multi-populations neural network models

In this section we introduce the classical neural mass models and compute the related mean field equations they satisfy in the limit of an infinite number of neurons

The general model

General framework

We consider a network composed of N neurons indexed by i ∈ {1, . . . , N } belonging to P populations indexed by α ∈ {1, . . . , P } such as those shown in figure 1. Let N α be the number of neurons in population α. We have N = P α=1 N α . In the following we are interested in the limit N → ∞. We assume that the proportions of neurons in each population are non-trivial, i.e. :

lim N →∞ N α N = n α ∈ (0, 1) ∀α ∈ {1, . . . , P }.
If it were not the case the corresponding population would not affect the global behavior of the system, would not contribute to the mean field equation, and could be neglected. We introduce the function p : {1, . . . , N } → {1, . . . , P } such that p(i) is the index of the population which the neuron i belongs to.

The following derivation is built after Ermentrout's review [START_REF] Ermentrout | Neural networks as spatio-temporal pattern-forming systems[END_REF]. We consider that each neuron i is described by its membrane potential V i (t) or by its instantaneous firing rate ν i (t), the relation between the two quantities being of the form ν i (t) = S i (V i (t)) (Dayan and Abbott, 2001;[START_REF] Gerstner | Mathematical formulations of hebbian learning[END_REF], where S i is sigmoidal. A single action potential from neuron j is seen as a post-synaptic potential P SP ij (ts) by neuron i, where s is the time of the spike hitting the synapse and t the time after the spike. We neglect the delays due to the distance travelled down the axon by the spikes.

Assuming that the post-synaptic potentials sum linearly, the average membrane potential of neuron i is

V i (t) = j,k P SP ij (t -t k ),
where the sum is taken over the arrival times of the spikes produced by the neurons j. The number of spikes arriving between t and t + dt is ν j (t)dt. Therefore we have

V i (t) = j t -∞ P SP ij (t -s)ν j (s) ds = j t -∞ P SP ij (t -s)S j (V j (s)) ds,
or, equivalently

ν i (t) = S i   j t t0 P SP ij (t -s)ν j (s) ds   . (1)

INRIA

The P SP ij s can depend on several variables in order to account for instance for adaptation or learning.

The voltage-based model

The assumption, made in [START_REF] Hopfield | Neurons with graded response have collective computational properties like those of two-state neurons[END_REF], is that the post-synaptic potential has the same shape no matter which presynaptic population caused it, the sign and amplitude may vary though. This leads to the relation

P SP ij (t) = J ij g i (t).
g i represents the unweighted shape (called a g-shape) of the postsynaptic potentials and J ij is the strength of the postsynaptic potentials elicited by neuron j on neuron i. Thus we have

V i (t) = t t0 g i (t -s)   j J ij ν j (s)   ds.
So far we have only considered the synaptic inputs to the neurons. We also assume that neuron i receives an external current density I i (t) and some noise n i (t) so that

V i (t) = t t0 g i (t -s)   j J ij ν j (s) + I i (s) + n i (s)   ds. (2)
We assume that the external current and the g-shapes satisfy I i = I p(i) , g i = g p(i) , S i = S p(i) , i.e. they only depend upon the neuron population. The noise model is described later. Finally we assume that g i = g α (where α = p(i)) is the Green function of a linear differential equation of order k, i.e. satisfies k l=0 a l α (t)

d l g α dt l (t) = δ(t). (3)
We assume that the functions a l α (t) are continuous for l = 0,

• • • , k and α = 1, • • • , P . We also assume a k α (t) ≥ c > 0 for all t ∈ Ê, α = 1, • • • , P . Known examples of g-shapes, see section 2.2.3 below, are g α (t) = Ke -t/τ Y (t) (k = 1, a 1 (t) = 1 K , a 0 (t) = 1 K τ )or g α (t) = Kte -t/τ Y (t) (k = 2, a 2 (t) = 1 K , a 1 (t) = 2 K τ , a 0 (t) = 1 K τ )
, where Y is the Heaviside function.

We note D k α the corresponding differential operator, D k α g α = δ, and

D k N the N -dimensional differential operator containing N α copies of D k α , α = 1, • • • , P . We write (2) in vector form V (N ) = J (N ) diag(g α ) * S (N ) (V (N ) ) + diag(g α ) * I (N ) + diag(g α ) * n (N ) , where diag(g α ) is the N -dimensional diagonal matrix containing N α copies of g α , α = 1, • • • , P and * indicates the convolution operator. S (N ) is the mapping Ê N → Ê N such that S (N ) (V (N ) ) i = S p(i) (V (N ) i which is a stochastic differential equation d d k-1 V (N ) dt k-1 = -D k-1 N V (N ) + J (N ) • S (N ) (V (N ) ) + I (N ) dt + dn (N ) t , where D k-1 N is obtained from the P differential operators of order k -1 D k-1 α = k-2 l=0 a l α (t) a k α (t) d l dt l .
The activity-based model If we make the assumption that the shape of a PSP depends only on the nature of the presynaptic cell, that is

P SP ij = J ij g j ,
and define the activity as

A j (t) = t -∞ g j (t -s)ν j (s) ds,
multiplying both sides of equation ( 1) by g i (ts) and integrating with respect to s, we obtain

A i (t) = t -∞ g i (t -s)S i   j J ij A j (s) + I i (s) + n i (s)   ds,
where we have added an external current and a noise. If p(i) = α, this yields

D k α A i = S i   j J ij A j (t) + I i (t) + n i (t)   ,
and in terms of the N -dimensional vector A (N )

D k N A (N ) = S (N ) (J (N ) A (N ) + I (N ) A + n (N ) A ). (5)
Equivalence of the two models As a matter of fact these two equations are equivalent provided that J (N ) is invertible1 . Indeed, let us use the change of variable

V (N ) = J (N ) A (N ) + I (N ) A + n (N )
A . We have, because J (N ) is not a function of time,

D k N V (N ) = J (N ) D k N A (N ) + D k N I (N ) A + D k N n (N )
A .

Replacing D k N V (N ) by this value in (4) we obtain

J (N ) D k N A (N ) + D k N I (N ) A + D k N n (N ) A = J (N ) • S (N ) (J (N ) A (N ) + I (N ) A + n (N ) A ) + I (N ) V + n (N )
V .

Assuming that the matrix J (N ) is invertible yields

D k N A (N ) = S (N ) (J (N ) A (N ) + I (N ) A + n (N ) A ) + (J (N ) ) -1 I (N ) V -D k N I (N ) A + n (N ) V -D k N n (N ) A .
Given the current I (N ) V

(respectively the noise n (N )

V ), we can choose the current

I (N ) A (respec- tively the noise n (N ) A ) solution of the linear differential equation D k N I (N ) A = I (N ) V (respectively D k N n (N ) A = n (N ) V ). Using the Green functions g α , α = 1, • • • , P this is equivalent to I (N ) A = diag(g α ) * I (N ) V (respectively n (N ) A = diag(g α ) * n (N ) V ).

The dynamics We introduce the

k -1 N -dimensional vectors V (l) (t) = [V (l) 1 , • • • , V (l) N ] T , l = 1, • • • , k -1 of the lth-order derivative of V (N ) (t)
, and the N k-dimensional vector

V (N ) (t) =      V (N ) (t) V (N ) (1) (t)
. . .

V (N ) (k-1) (t)      .
The N -neurons network is described by the N k-dimensional vector V (N ) (t). We consider the 1) , where each E (l) = Ê N , l = 0, • • • , k -1 and introduce the following notation: if x is a vector of Ê N k , x l is its component in E (l) , l = 0, • • • , k -1, an N -dimensional vector. In particular we have N ) .

direct sum Ê N k = E (0) ⊕ • • • ⊕ E (k-
V (N ) l (t) = V (N ) (l) (t) for l = 0, • • • , k -1 with the convention that V (N ) (0) = V (
We now write the equations governing the time variation of the first k -1 vectors of V (N ) (t), i.e. the derivatives of order 0, . . . k -2 of V (N ) (t). These equation in effect determine the noise model. We write

d V (N ) l (t) = V (N ) l+1 (t) dt + Λ (N ) l • dW (N ) t l = 0, • • • , k -2, (6)
where Λ

(N ) l is the N × N diagonal matrix diag(s l α ), where s l α , α = 1, • • • , P is repeated N α times, and W (N ) t an N -dimensional standard Brownian process. The equation governing the (k -1)th differential of the membrane potential has a linear part determined by the differential operator D k-1 and must account for the external inputs (deterministic and stochastic) and the activity of the neighbors, see (4). Keeping the same notations as before for the inputs and denoting by L (N ) the N × N k matrix describing the action of the neurons membrane potentials and their derivatives on the (k -1)th derivative of V, we have:

d V (N ) k-1 (t) = L (N ) (t) • V (N ) (t) + J (N ) • S (N ) ( V (N ) 0 (t)) + I (N ) (t) dt + Λ (N ) k-1 (t) • dW (N ) t , (7)
where

L (N ) = diag(a 0 α (t)) • • • diag(a k-1 α (t)) . We define L (N ) (t) =       0 Id N • • • 0 0 0 . . . 0 . . . . . . Id N diag(a 0 α (t)) diag(a 1 α (t)) • • • diag(a k-1 α (t))       ,
where Id N is the N × N identity matrix. We also denote by:

U (N ) t =      0 . . . 0 J (N ) • S (N ) ( V 0 (t))      and I (N ) t =      0 . . . 0 I (N ) (t)      .
The full equation satisfied by V (N ) can be written:

d V (N ) (t) = L (N ) (t) V (N ) (t) + U (N ) t + I (N ) t dt + Λ (N ) (t) • dW (N ) t , (8)
where the kN × kN matrix

Λ (N ) (t) is equal to diag(Λ (N ) 0 , • • • , Λ (N )
k-1 ). Note that the kth-order differential equation describing the time variation of the membrane potential of each neuron contains a noise term which is a linear combination of various integrated Brownian processes (up to the order k -1) as shown in the following formula which is derived from (6) and (7).

d V (N ) k-1 (t) = k-1 l=0 L (N ) l (t) V (N ) l dt+ k-2 l=0 L (N ) l (t) k-l-2 h=0 t 0 s1 0 • • • s h-1 0 Λ (N ) l+h (s h )dW s h d s h-1 • • • ds 0 dt+ Λ (N ) k-1 (t)dW t .
Comparing with equation (4) we see that the noise n (N ) dt is a weighted sum of Brownian and integrated Brownian processes.

The Mean Field equations

General derivation of the mean field equation

The connectivity weight J ij are modeled as independent Gaussian random variables. Their distribution depends only on the population pair α = p(i), β = p(j), and on the total number of neurons N β of population β:

INRIA J ij ∼ N Jαβ N β , σ αβ N β .
We are interested in the limit law when N → ∞ of the vector V (N ) under the joint law of the connectivities and the Brownian motions, which we call the mean field limit. This law can be described by a set of P equations, the mean field equations. As mentioned in the introduction these equations can be derived in several ways, either heuristically as in the work of Amari [START_REF] Amari | Characteristics of random nets of analog neuron-like elements[END_REF][START_REF] Amari | A mathematical foundation for statistical neurodynamics[END_REF], Sompolinsky [START_REF] Crisanti | chaos in neural networks : chaotic solutions Dayan P, Abbott LF[END_REF][START_REF] Sompolinsky | Chaos in Random Neural Networks[END_REF], and Cessac [START_REF] Cessac | Increase in complexity in random neural networks[END_REF][START_REF] Samuelides | Random recurrent neural networks[END_REF], or rigorously as in the work of Benarous and Guionnet (Ben-Arous and [START_REF] Ben-Arous | Large deviations for Langevin spin glass dynamics[END_REF][START_REF] Guionnet | Averaged and quenched propagation of chaos for spin glass dynamics[END_REF][START_REF] Guionnet | Averaged and quenched propagation of chaos for spin glass dynamics[END_REF]. We derive them here in a pedestrian way, prove that they are well-posed, and provide an algorithm for computing their solution.

The effective description of the network population by population is possible because the neurons in each population are interchangeable, i.e. have the same probability distribution under the joint law of the multidimensional Brownian motion and the connectivity weights. This is the case because of the form of equation ( 8).

The Mean Field equations We note C([t 0 , T ], Ê P ) (respectively C((-∞, T ], Ê P )) the set of continuous functions from the real interval [t 0 , T ] (respectively (-∞, T ]) to Ê P . By assigning a probability to subsets of such functions, a continuous stochastic process X defines a positive measure of unit mass on C([t 0 , T ], Ê P ) (respectively C((-∞, T ], Ê P )). This set of positive measures of unit

mass is noted M + 1 (C([t 0 , T ], Ê P )) (respectively M + 1 (C((-∞, T ], Ê P )).
We now define a process of particular importance for describing the limit process: the effective interaction process. Definition 2.1 (Effective Interaction Process). Let X ∈ M + 1 (C([t 0 , T ], Ê P )) (resp. M + 1 (C((-∞, T ], Ê P ))

be a given stochastic process. The effective interaction term is the Gaussian process U X ∈ M + 1 (C([t 0 , T ], Ê P ×P )),

(resp. M + 1 (C((-∞, T ], Ê P ×P )) statistically independent of the external noise (W t ) t≥t0 and of the initial condition X t0 (when t 0 > -∞), defined by:

           E U X αβ (t) = Jαβ m X αβ (t) where m X αβ (t) def = E[S αβ (X β (t))]; Cov(U X αβ (t), U X αβ (s)) = σ 2 αβ ∆ X αβ (t, s) where ∆ X αβ (t, s) def = E S αβ (X β (t))S αβ (X β (s)) ; Cov(U X αβ (t), U X γδ (s)) = 0 if α = γ or β = δ. (9)
Choose P neurons i 1 , . . . , i P , one in each population (neuron i α belongs to the population α). Then it can be shown, using either a heuristic argument or large deviations techniques (see appendix A), that the sequence of processes

V (N ) (t) = [ Ṽ (N ) i1 (t), . . . , Ṽ (N ) iP (t)] T t≥t0 N ≥1
converges in law to the process V(t) = [ Ṽ1 (t), . . . , ṼP (t)] T t≥t0 solution of the following mean field equation:

d V(t) = L(t) V(t) + U V t + I(t) dt + Λ(t) • dW t , ( 10
)
where V is a kP -dimensional vector containing the P -dimensional vector V and its k-1 derivatives, and L is the P k × P k matrix

L(t) =       0 P ×P Id P • • • 0 P ×P 0 P ×P 0 P ×P . . . 0 P ×P . . . . . . Id P L 0 (t) L 1 (t) • • • L k-1 (t)       ,
where Id P is the P × P identity matrix and 0 P ×P the null P × P matrix. (W t ) is a kP -dimensional standard Brownian process and:

U V t =      0 P . . . 0 P U V t • 1      I(t) =      0 P . . . 0 P I(t)      Λ(t) = diag(Λ 0 (t), • • • , Λ k-1 (t)). The matrices L 0 , • • • , L k-1 (respectively Λ 0 , • • • , Λ k-1
) are obtained by selecting the same P rows and P columns of the matrices L

(N ) 0 , • • • , L (N ) k-1 (respectively Λ (N ) 0 , • • • , Λ (N ) k-1 ) corresponding to P neurons in different populations, (U V t )
is the effective interaction process associated with V, and I(•) is the P -dimensional external current.

To proceed further we formally integrate the equation using the flow, or resolvent, of the equation, noted Φ L (t, t 0 ) (see appendix B), and we obtain, since we assumed L continuous, an implicit representation of V:

V(t) = Φ L (t, t 0 ) V(t 0 ) + t t0 Φ L (t, s) • U V s + I(s) ds + t t0 Φ L (t, s) • Λ(s)dW s (11)
We now introduce for future reference a simpler model which is quite frequently used in the description on neural networks.

Example: The Simple Model

In the Simple Model, each neuron membrane potential decreases exponentially to its rest value if it receives no input, with a time constant τ α depending only on the population. The noise is modeled by an independent Brownian process per neuron whose standard deviation is the same for all neurons belonging to a given population.

Hence the dynamics of a given neuron i from population α of the network reads:

INRIA dV (N ) i (t) = - V (N ) i (t) τ α + P β=1 N β j=1 J ij S αβ V (N ) j (t) + I α (t) dt + s α dW (i) (t). ( 12
)
This is a special case of equation ( 10) where k = 1 and

L = -diag( 1 τ1 , • • • , 1 τP ), Φ L (t, t 0 ) = diag(e -(t-t0)/τ1 , • • • , e -(t-t0)/τP ), and Λ = diag(s 1 , • • • , s P ).
The corresponding mean field equation reads:

dV α (t) = - V α τ α (t) + P β=1 U V αβ (t) + I α (t) dt + s α dW α (t), ∀α ∈ {1, . . . , P }, (13)
where the processes (W α (t)) t≥0 are independent standard Brownian motions,

U V (t) = (U V αβ (t); α, β ∈ {1, . . . , P }) t is the effective interaction term.
This equation can be integrated implicitly and we obtain the following integral representation of the process V α (t):

V α (t) = e -(t-t0)/τα V α (t 0 ) + t t0 e -(t-s)/τα P β=1 U V αβ (s) + I α (s) ds + s α t t0 e -(t-s)/τα dW α (s) ( 14
)
where t 0 is the initial time. It is an implicit equation on the probability distribution of V(t), a special case of (11).

The Jansen and Rit's model

One of the motivations of this study is to characterize the global behavior of an assembly of neurons in particular to get a better understanding of non-invasive cortical signals like EEG or MEG. One of the classical models of neural masses is Jansen and Rit's mass model [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF], in short the JR model (see figure 1). The model features a population of pyramidal neurons (central part of figure 1.a.) that receives excitatory and inhibitory feedback from local inter-neurons and an excitatory input from neighboring cortical units and sub-cortical structures such as the thalamus. The excitatory input is represented by an arbitrary average firing rate p(t) that can be stochastic (accounting for a non specific background activity) or deterministic, accounting for some specific activity in other cortical units. The transfer functions h e and h i of figure 1 convert the average firing rate describing the input to a population into an average excitatory or inhibitory post-synaptic potential (EPSP or IPSP). They correspond to the synaptic integration.

In the model introduced originally by Jansen and Rit, the connectivity weights were assumed to be constant, equal to their mean value (it is the constants C i , i = 1 . . . 4 in figure 1). Nevertheless, there exists a variability on these coefficients, and as we will see in the sequel, the effect of the connectivity variability impacts the solution at the level of the neural mass. Statistical properties INRIA of the connectivities have been studied in details for instance in [START_REF] Braitenberg | Cortex: Statistics and Geometry of Neuronal Connectivity[END_REF]. In our model we consider these connectivities as independent Gaussian random variables of mean and standard deviation equal to the ones found in [START_REF] Braitenberg | Cortex: Statistics and Geometry of Neuronal Connectivity[END_REF].

We now use diagram 1 to derive the membrane potential expressions. We consider a network of N neurons belonging to the three populations described. We denote by P (resp E, I) the pyramidal (respectively excitatory, inhibitory) populations. We choose in population P (respectively populations E, I) a particular pyramidal neuron (respectively excitatory, inhibitory interneuron) indexed by i pyr (respectively i exc , i inh ). The equations of their activity variable read:

               A N ipyr = h e * S( j Exc J ij A N j + h e * p(•) + j Inh J ij A N j ) A N iexc = h e * S( j P yr J ij A N j ) A N i inh = h i * S( j P yr J ij A N j )
This is therefore an activity-based model. As stated before, it is equivalent via a change of variable to a voltage-based model, with the same connectivity matrix, the same intrinsic dynamics, and modified inputs (see section 2.1.1).

In the mean field limit, denoting by A P (respectively A E , A I ) the activity of the pyramidal neurons (resp excitatory, inhibitory interneurons), we obtain the following activity equations:

     A P = h e * S(U P E + h e * p + U P I ) A E = h e * S(U EP ) A I = h e * S(U IP ) (15)
where U = (U ij ) i,j∈{P,E,I} is the effective interaction process associated with this problem, i.e. a Gaussian process of means:

         E [U EP ] = JEP E [A E ] E [U IP ] = JIP E [A I ] E [U P I ] = JP I E [A P ] E [U P E ] = JP E E [A P ]
and whose covariance matrix can be deduced from ( 9). The voltage-based model can be deduced from this activity-based description using a simple change of variable as stated previously. Note that the change of variable is possible since the activity current I A is equal to h e * p and, as shown in section 2.1.1, I A is smooth enough so that we can apply to it the suitable differential operator. p is the corresponding voltage current I V .

Let us now instantiate the synaptic dynamics and compare the mean field equation with Jansen's population equations (sometimes improperly called also mean field equations).

The simplest model of synaptic integration is a first-order integration, which yields exponential post-synaptic potentials:

h(t) = αe -βt t ≥ 0 0 t < 0
that satisfies the following differential equations

ḣ(t) = -βh(t) + αδ(t),
In these equations β is the time constant of the synaptic integration and α the synaptic efficiency. The coefficients named α and β are the same for the pyramidal and the excitatory population, and different from the ones of the inhibitory synapse. In the pyramidal or excitatory (respectively the inhibitory) case we have α = A, β = a (respectively α = B, β = b). Eventually, the sigmoid functions are the same whatever the populations, and is given by v0-v) , ν max is the maximum firing rate, and v 0 is a voltage reference.

S(v) = ν max 1 + e r(
With this synaptic dynamics we obtain the first-order Jansen and Rit's equation:

     dAP dt (t) = -a A P (t) + A S(U P E + U P I + h e * p(t)) dAE dt (t) = -a A E (t) + A S(U EP ) dAI dt (t) = -b A I (t) + B S(U IP ) . ( 16
)
while the "original" Jansen and Rit's equation [START_REF] Grimbert | Bifurcation analysis of Jansen's neural mass model[END_REF][START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF] reads:

     dAP dt (t) = -a A P (t) + A S(C 2 A E (t) -C 4 A I (t) + h e * p(t)) dAE dt (t) = -a A E (t) + A S(C 1 A P (t)) dAI dt (t) = -b A I (t) + B S(C 3 A P (t)) . ( 17
)
Hence the original JR equation amounts to computing the expectation of the activity in each population and to assume that

E [S(U P E + U P I + h e * p(t))] = S(E [U P E + U P I + h e * p]),
which is a quite sharp assumption given that the sigmoidal function is nonlinear.

A higher order model was introduced to better account for the synaptic integration and to better reproduce the characteristics of real EPSPs and IPSPs by van Rotterdam and colleagues [START_REF] Van Rotterdam | A model of the spatial-temporal characteristics of the alpha rhythm[END_REF]. In this model the PSP satisfies a second order differential equation:

h(t) = αβte -βt t ≥ 0 0 t < 0 , INRIA solution of the differential equation ÿ(t) = αβδ(t) -2β ẏ(t) -β 2 y(t).
With this type of synaptic integration, we obtain the following mean field equations:

     d 2 AP dt 2 (t) = AaS(U P E + U P I + h e * p(t)) -2a dAP dt (t) -a 2 A P (t) d 2 AE dt 2 (t) = AaS(U EP ) -2a dAE dt (t) -a 2 A E (t) d 2 AI dt 2 (t) = BbS(U IP ) -2b dAI dt (t) -b 2 A I (t) (18)
while the original system satisfies the equations:

     d 2 AP dt 2 (t) = AaS(C 2 A E (t) -C 4 A I (t) + h e * p(t)) -2a dAP dt (t) -a 2 A P (t) d 2 AE dt 2 (t) = AaS(C 1 A P (t)) -2a dAE dt (t) -a 2 A E (t) d 2 AI dt 2 (t) = BbC 4 S(C 3 A P (t)) -2b d 2 AI dt 2 (t) -b 2 A I (t) (19)
Here again, going from the mean field equations ( 18) to the neural mass model ( 19) consists in studying the equation of the mean of the process given by ( 18) and commuting the sigmoidal function with the expectation. Note that the introduction of higher order synaptic integrations results in richer behaviors. For instance, Grimbert and Faugeras [START_REF] Grimbert | Bifurcation analysis of Jansen's neural mass model[END_REF] showed that some bifurcations can appear in the second-order JR model giving rise to epileptic like oscillations and alpha activity, that do not appear in the first order model.

Existence and uniqueness of solutions in finite time

The mean field equation ( 11) is an implicit equation of the stochastic process (V (t)) t≥t0 . We prove in this section that under some mild assumptions this implicit equation has a unique solution. This solution is a fixed point in the set M + 1 (C([t 0 , T ], Ê kP )) of kP -dimensional processes. We construct a sequence of processes and prove that it converges in distribution toward this fixed point.

We denote by X the set of random variables (r.v.) with values in Ê kP . We first recall some results on the convergence of random variables and stochastic processes.

Convergence of Gaussian processes

We recall the following result from [START_REF] Bogachev | Gaussian Measures[END_REF].

Theorem 3.1. Let {X n } ∞ n=1 be a sequence of kP -dimensional Gaussian processes defined on [t 0 , T ]
or on an unbounded interval of Ê2 . The sequence converges to a Gaussian process X if and only if the following three conditions are satisfied:

• The sequence {X n } ∞ n=1 is uniformly tight.

• The sequence µ n (t) of the mean functions converges for the uniform norm.

• The sequence C n of the covariance operators converges for the uniform norm.

We now define such a sequence of Gaussian processes.

Let us fix Z 0 , a kP -dimensional Gaussian random variable, independent of the Brownian.

Definition 3.1. Let X an element of M + 1 (C([t 0 , T ], Ê kP )). Let F k be the function such that F k (X) t = Φ L (t, t 0 ) • Z 0 + t t0 Φ L (t, s) • U X s + I(s) ds + t t0 Φ L (t, s) • Λ(s)dW s
where U X s and I(s) are defined in section 2. Note that, by definition, the random process (F k (X)) t∈[t0,T ] , k ≥ 1 is the sum of a deterministic function (defined by the external current) and three independent random processes defined by the initial condition, the interaction between neurons, and the external noise. Let X be a given stochastic process of

M + 1 (C([t 0 , T ], Ê kP )) such that X t0 = Z 0 . We define the sequence of processes {X n } ∞ n=0 ∈ M + 1 (C([t 0 , T ], Ê kP ) by: X 0 = X X n+1 = F k (X n ) = F (n) k (X 0 ). (20)
In the remaining of this section we show that the sequence of processes

{F (n) k (X)} ∞ n=0 converges in distribution toward the unique fixed-point Y of F k .

Existence and uniqueness of solution for the mean field equations

The following upper and lower bounds are used in the sequel. Lemma 3.2. We consider the Gaussian process ((U X t • 1) t ) t∈[t0,T ] . U X is defined in 2.1 and 1 is the P -dimensional vector with all coordinates equal to 1. We have

E U X t • 1 ∞ ≤ µ def = max α β | Jαβ | S αβ ∞ (21)
for all t 0 ≤ t ≤ T . The maximum eigenvalue of its covariance matrix is upperbounded by

σ 2 max def = max α β σ 2 αβ S αβ 2 ∞
where S αβ ∞ is the supremum of the absolute value of S αβ .

Proof. The proof is straightforward from definition 3.1. We also note σ 2 min def = min α,β σ 2 αβ . The proof of existence and uniqueness of solution, and of the convergence of the sequence (20) is in two main steps. We first prove that the sequence of Gaussian processes

{F (n) k (X)} ∞ n=0
, k ≥ 1 is uniformly tight by proving that Kolmogorov's criterion for tightness holds. This takes care of condition 1) in theorem 3.1. We next prove that the sequences of the mean functions and covariance operators are Cauchy sequences for the uniform norms, taking care of conditions 2) and 3).

INRIA

Theorem 3.3. The sequence of processes

F (n) k (X) ∞ n=0 , k ≥ 1 is uniformly tight.
Proof. We use Kolmogorov's criterion for tightness and do the proof for k = 1, the case k > 1 is similar. If we assume that n ≥ 1 and s < t we have

F (n) 1 (X) t -F (n) 1 (X) s = (Φ L (t, t 0 ) -Φ L (s, t 0 ))X t0 + s t0 (Φ L (t, s) -Id)Φ L (s, u)U F (n-1) 1 (X) u • 1 du + t s Φ L (t, u)U F (n-1) 1 (X) u • 1 du + s t0 (Φ L (t, s) -Id)Φ L (s, u)Λ(u) dW u + t s Φ L (t, u)Λ(u) dW u + s t0 (Φ L (t, s) -Id)Φ L (s, u)I(u) du + t s Φ L (t, u)I(u) du
and therefore (Cauchy-Schwarz and Jensen's inequalities):

1 7 F (n) 1 (X) t -F (n) 1 (X) s 2 ≤ Φ L (t, t 0 ) -Φ L (s, t 0 ) 2 X t0 2 + (s -t 0 ) Φ L (t, s) -Id 2 s t0 Φ L (s, u) 2 U F (n-1) 1 (X) u • 1 2 du + (t -s) t s Φ L (t, u) 2 U F (n-1) 1 (X) u • 1 2 du + s t0 Φ L (s, u)(Φ L (t, s) -Id)Λ(u) dW u 2 + t s Φ L (t, u)Λ(u) dW u 2 + (s -t 0 ) 2 Φ L (t, s) -Id 2 I 2 max sup u∈[t0,s] Φ L (s, u) 2 + (t -s) 2 I 2 max sup u∈[s,t] Φ L (t, u) 2 .
Because Φ L (t, t 0 ) -Φ L (s, t 0 ) ≤ |t -s| L we see that all terms in the righthand side of the inequality but the two involving the Brownian motion are of the order of (ts) 2 . We raise again both sides to the second power, use the Cauchy-Schwarz inequality again, and take the expected value:

1 7 3 E F (n) 1 (X) t -F (n) 1 (X) s 4 ≤ Φ L (t, t 0 ) -Φ L (s, t 0 ) 4 E X t0 4 + (s -t 0 ) 3 Φ L (t, s) -Id 4 s t0 Φ L (s, u) 4 E U F (n-1) 1 (X) u • 1 4 du + (t -s) 3 t s Φ L (t, u) 4 E U F (n-1) 1 (X) u • 1 4 du + E s t0 Φ L (s, u)(Φ L (t, s) -Id)Λ(u) dW u 4 + E t s Φ L (t, u)Λ(u) dW u 4 (22) + (s -t 0 ) 4 Φ L (t, s) -Id 4 sup u∈[t0,s] Φ L (s, u) 4 I 4 max + (t -s) 4 I 4 max sup u∈[s,t] Φ L (t, u) 4 .
Remember that U

F (n-1) 1 (X) u • 1 is a P -dimensional diagonal Gaussian process, noted Y u , there- fore: E Y u 4 = α E Y α (u) 4 + α1 =α2 E Y 2 α1 (u) E Y 2 α2 (u) .
The second order moments are upperbounded by some regular function of µ and σ max (defined in lemma 3.2) and, because of the properties of Gaussian integrals, so are the fourth order moments.

Let us now evaluate E b a A(u) dW u 4

for some P × P matrix A. We have

E   b a A(u) dW u 4   = E      b a A(u) dW u 2   2    = E      P i=1   P j=1 b a A ij (u) dW j u   P k=1 b a A ik (u) dW k u   2    = i1,i2,j1,j2,k1,k2 E b a A i1j1 (u) dW j1 u b a A i1k1 (u) dW k1 u b a A i2j2 (u) dW j2 u b a A i2k2 (u) dW k2 u .

INRIA

Because of the properties of the Brownian process, the last term is the sum of three types of terms:

i1,i2 E   b a A i1j (u) dW j u 2 b a A i2j (u) dW j u 2   ≤ i1,i2 E   b a A i1j (u) dW j u 4   1/2 E   b a A i2j (u) dW j u 4   1/2
, and

i1,i2,j1 =j2 E   b a A i1j1 (u) dW j1 u 2 b a A i2j2 (u) dW j2 u 2   = i1,i2,j1 =j2 E   b a A i1j1 (u) dW j1 u 2   E   b a A i2j2 (u) dW j2 u 2   ,
A i1j (u) dW j u = N (0, b a A 2 i1j (u) du 1/2
) hence, because of the properties of the Gaussian integrals

E   b a A i1j (u) dW j u 4   = k b a A 2 i1j (u) du 2 , for some positive constant k. Moreover E   b a A i1j1 (u) dW j1 u 2   = b a A 2 i1j1 (u) du, and
E b a A i1j1 (u) dW j1 u b a A i2j1 (u) dW j1 u = b a A i1j1 (u)A i2j1 (u) du.
This shows that the two terms 22) are of the order of (ts) 1+a where a ≥ 1. Therefore we have

E s t0 (Φ L (t, s) -Id)Φ L (s, u)Λ(u) dW u 4 and E t s Φ L (t, u)Λ(u) dW u 4 in (
E F (n) 1 (X) t -F (n) 1 (X) s 4 ≤ C|t -s| 1+a , a ≥ 1 for all s, t in [t 0 , T ],
where C is a constant independent of t, s. According to Kolmogorov criterion for tightness, the sequence of processes

F (n) 1 (X) ∞ n=0 is uniformly tight. The proof for F k , k > 1 is similar.
Let us note µ n (t) (respectively C n (t, s)) the mean (respectively the covariance matrix) function of X n = F k (X n-1 ), n ≥ 1. We have:

C n+1 (t, s) = Φ L (t, t 0 )Σ Z0 Φ L (s, t 0 ) T + t∧s t0 Φ L (t, u)Λ(u)Λ(u) T Φ L (s, u) T du+ t t0 s t0 Φ L (t, u)cov U Xn u , U Xn v Φ L (s, v) T du dv (23)
Note that the kP × kP covariance matrix cov U Xn u , U Xn v has only one nonzero P × P block:

cov U Xn u , U Xn v kk = cov U Xn u • 1, U Xn v • 1 , We have cov U Xn u • 1, U Xn v • 1 = diag   β σ 2 αβ E [S αβ (X nβ (u))S αβ (X nβ (v))]   ,
and

E [S αβ (X n (u))S αβ (X n (v))] = Ê 2 S αβ C n ββ (u, u)C n ββ (v, v) -C n ββ (u, v) 2 C n ββ (u, u) x + C n ββ (u, v) C n ββ (u, u) y + µ n β (v) S αβ y C n ββ (u, u) + µ n β (u) Dx Dy, (24)
where Dx = 1 √ 2π e -x 2 2 dx.

INRIA

Similarly we have

µ n (t) = Φ L (t, t 0 )µ Z0 + t t0 Φ L (t, u) E U Xn u + Ĩ(u) du = Φ L (t, t 0 )µ Z0 + t t0 Φ L (t, u)      0 T P , • • • , 0 T P ,   β Jαβ Ê S αβ x C n ββ (u, u) + µ n β (u)   α=1,••• ,P   T + Ĩ(u)    Dxdu
We require the following four lemmas.

Lemma 3.4. For all α = 1, • • • , P and n ≥ 1 the quantity C n αα (s, s)C n αα (t, t) -C n αα (t, s) 2 is lowerbounded by the positive symmetric function:

θ(s, t) = |t -s|λ 2 min λ Σ Z 0 min λ Γ min ,
where λ min is the smallest singular value of the positive symmetric definite matrix

Φ L (t, t 0 )Φ L (t, t 0 ) T for t ∈ [t 0 , T ], λ Σ Z 0
min is the smallest eigenvalue of the positive symmetric definite covariance matrix Σ Z0 , and λ Γ min is the smallest singular value of the matrix Λ(u) for u ∈ [t 0 , T ]. Proof. We use equation ( 23) which we rewrite as follows, using the group property of the resolvent Φ L :

C n+1 (t, s) = Φ L (t, t 0 ) Σ Z0 + t∧s t0 Φ L (t 0 , u)Λ(u)Λ(u) T Φ L (t 0 , u) T du+ t t0 s t0 Φ L (t 0 , u)cov U Xn u , U Xn v Φ L (t 0 , v) T du dv Φ L (t 0 , s) T .
We now assume s < t and introduce the following notations, dropping the index n for simplicity:

A(s) = Σ Z0 + s t0 Φ L (t 0 , u)Λ(u)Λ(u) T Φ L (t 0 , u) T du B(s, t) = t s Φ L (t 0 , u)Λ(u)Λ(u) T Φ L (t 0 , u) T du a(t, s) = t t0 s t0 Φ L (t 0 , u)cov U Xn u , U Xn v Φ L (t 0 , v) T du dv
Let e α , α = 1, • • • , kP , be the unit vector of the canonical basis whose coordinates are all equal to 0 except the αth one which is equal to 1. We note E α (t) the vector Φ L (t, t 0 ) T e α . We have

C αα (t, s) = E α (t) T (A(s) + a(t, s)) E α (s) C αα (s, s) = E α (s) T (A(s) + a(s, s))E α (s) C αα (t, t) = E α (t) T (A(s) + B(s, t) + a(t, t))E α (t).
Note that the last expression does not depend on s, since A(s) + B(s, t) = A(t), which is consistent with the first equality. The reason why we introduce s in this expression is to simplify the following calculations.

The expression C αα (s, s)C αα (t, t) -C αα (t, s) 2 is the sum of four sub-expressions:

e 1 (s, t) = E α (s) T A(s)E α (s) E α (t) T A(s)E α (t) -E α (t) T A(s)E α (s) 2 ,
which is greater than or equal to 0 because A(s) is a covariance matrix,

e 2 (s, t) = E α (s) T a(s, s)E α (s) E α (t) T a(t, t)E α (t) -E α (t) T a(t, s)E α (s) 2 ,
which is also greater than or equal to 0 because a(t, s) is a covariance matrix function,

e 3 (s, t) = E α (s) T A(s)E α (s) E α (t) T a(t, t)E α (t) + E α (t) T A(s)E α (t) E α (s) T a(s, s)E α (s) - 2 E α (t) T A(s)E α (s) E α (t) T a(t, s)E α (s)
Because a(t, s) is a covariance matrix function we have

E α (t) T a(t, t)E α (t) + E α (s) T a(s, s)E α (s) -2E α (t) T a(t, s)E α (s) ≥ 0,
and , as seen above, e 2 (s, t) ≥ 0. Because e 1 (s, t) ≥ 0 we also have

-E α (s) T A(s)E α (s) E α (t) T A(s)E α (t) ≤ E α (t) T A(s)E α (s) ≤ E α (s) T A(s)E α (s) E α (t) T A(s)E α (t),
and, as it can be readily verified, this implies e 3 (s, t) ≥ 0.

Therefore we can lowerbound C αα (s, s)C αα (t, t) -C αα (t, s) 2 by the fourth subexpression:

C αα (s, s)C αα (t, t) -C αα (t, s) 2 ≥ E α (s) T A(s)E α (s) E α (t) T B(s, t)E α (t) + E α (s) T a(s, s)E α (s) E α (t) T B(s, t)E α (t) ≥ E α (s) T A(s)E α (s) E α (t) T B(s, t)E α (t) ,
since B(s, t) and a(s, s) are covariance matrixes. We next have

E α (s) T A(s)E α (s) = E α (s) T A(s)E α (s) E α (s) T E α (s) e T α Φ L (s, t 0 )Φ L (s, t 0 ) T e α e T α e α , by definition of E α (s). Therefore E α (s) T A(s)E α (s) ≥ λ A(s) min λ ΦL(s,t0)ΦL(s,t0) T min ≥ λ Σ Z 0 min λ min ,
where λ C min is the smallest eigenvalue of the symmetric positive matrix C. Similarly we have

E α (t) T B(s, t)E α (t) ≥ λ B(s,t) min λ min . INRIA Let us write Γ(u) = Λ(u)Λ(u) T . We have λ B(s,t) min = min x ≤1 t s x T Φ L (t 0 , u)Γ(u)Φ L (t 0 , u) T x x T x du = min x ≤1 t s x T Φ L (t 0 , u)Γ(u)Φ L (t 0 , u) T x x T Φ L (t 0 , u)Φ L (t 0 , u)x x T Φ L (t 0 , u)Φ L (t 0 , u)x x T x du ≥ t s min x ≤1 x T Φ L (t 0 , u)Γ(u)Φ L (t 0 , u) T x x T Φ L (t 0 , u)Φ L (t 0 , u)x x T Φ L (t 0 , u)Φ L (t 0 , u)x x T x du ≥ (t -s)λ min λ Γ min .
Combining these results we have

C αα (s, s)C αα (t, t) -C αα (t, s) 2 ≥ |t -s|λ 2 min λ Σ Z 0 min λ Γ min Lemma 3.5. For all t ∈ [t 0 , T ] all α = 1, • • • , P , and n ≥ 1, we have C n αα (t, t) ≥ k 0 > 0. Proof. C n αα (t, t) is larger than (Φ L (t, t 0 )Σ Z0 Φ L (t, t 0 ) T )
αα which is larger than the smallest eigenvalue of the matrix Φ L (t, t 0 )Σ Z0 Φ L (t, t 0 ) T . This smallest eigenvalue is equal to

min x x T Φ L (t, t 0 )Σ Z0 Φ L (t, t 0 ) T x x T x = min x x T Φ L (t, t 0 )Σ Z0 Φ L (t, t 0 ) T x x T Φ L (t, t 0 )Φ L (t, t 0 ) T x x T Φ L (t, t 0 )Φ L (t, t 0 ) T x x T x ≥ min x x T Φ L (t, t 0 )Σ Z0 Φ L (t, t 0 ) T x x T Φ L (t, t 0 )Φ L (t, t 0 ) T x min x x T Φ L (t, t 0 )Φ L (t, t 0 ) T x x T x .
In the last expression the first term is larger than the smallest eigenvalue λ Σ Z 0 min of the matrix Σ Z0 which is positive definite since we have assumed the Gaussian random variable Z 0 nondegenerate. The second term is equal to the smallest singular value λ min of the matrix Φ L (t, t 0 ) which is also strictly positive for all t ∈ [t 0 , T ] by hypothesis, see appendix B, equation (47).

We also use the following lemma.

Lemma 3.6. The 2n-dimensional integral

I n = [t0,t∨s] 2 ρ 1 (u 1 , v 1 ) [t0,u1∨v1] 2 • • • [t0,un-2∨vn-2] 2 ρ n-1 (u n-1 , v n-1 ) [t0,un-1∨vn-1] 2 ρ n (u n , v n )du n dv n du n-1 dv n-1 • • • du 1 dv 1 , RR n°6454
where the functions ρ i (u i , v i ), i = 1, • • • , n are either equal to 1 or to 1/ θ(u i , v i ), is upperbounded by k n /(n -1)! for some positive constant k.

Proof. First note that the integral is well-defined because of lemma 3.4. Second, note that there exists a constant K such that K/ θ(u, v) ≥ 1 for all (u, v) ∈ [t 0 , t∨s] 2 , i.e. K = λ min λ Σ Z 0 min λ Γ min (Tt 0 ). Therefore the integral is upperbounded by K n 0 , where

K 0 = max(1, K) times the integral obtained when ρ i (u i , v i ) = 1/ |u i -v i | for all i = 1, • • • , n.
Let us then consider this situation. Without loss of generality we assume t 0 = 0. The cases n = 1, 2, 3 allow one to understand the process.

I 1 ≤ K 0 [0,t∨s] 2 dudv |u -v| . ( 25
)
Let us rotate the axes by -π 4 by performing the change of variables

u = U + V √ 2 , v = V -U √ 2 .
Using the symmetry of the integrand in s and t and the change of variable, the integral in the righthand side of ( 25) is equal to (see figure 2):

2 1 2 1/4 t∨s √ 2 0 √ 2(t∨s)-U U dV dU √ U = 2 3/4 a/2 0 a -2U √ U dU = 2 3/4 α 1 a 3/2 ,
where a = √ 2(t ∨ s) and α 1 = 2 √ 2 3 . Let us now look at I 2 . It is upperbounded by the factor K 2 0 (2 3/4 ) 2 α 1 times the integral

a/2 0 a-U U ( √ 2(u ∨ v)) 3/2 √ U dU dV.
Since in the area of integration u ∨

v = v = V -U √ 2
we are led to the product of 2/5 by the onedimensional integral

a/2 0 (a -2U ) 5/2 √ U dU dV = α 2 a 3 , where α 2 = 5 √ 2π 32 . Similarly I 3 is upperbounded by the product of K 3 0 (2 3/4 ) 3 α 1 α 2 2 5 2 8 times the integral a/2 0 (a -2U ) 4 √ U dU dV = α 3 a 9/2 ,
where α 3 = 128 √ 2 315 . One easily shows then that:

I n ≤ K n 0 F (2 3/4 ) n 2 n n i=1 α i 1 n j=1 (2 + 3(j -1))
.

INRIA

RR n°6454

It can be verified by using a system for symbolic computation that 0 < α i < 1 for all i ≥ 1. One also notices that n j=1

(2 + 3(j -1)) ≥

3 n-1 2 (n -1)!, therefore I n ≤ K n 0 (2 3/4 ) n 2 n-1 3 -(n-1) 1 (n -1)! ,
and this finishes the proof.

We now prove the following proposition.

Proposition 3.7. The sequences of covariance matrix functions C n (t, s) and of mean functions µ n (t), s, t in [t 0 , T ] are Cauchy sequences for the uniform norms.

Proof. We have

C n+1 (t, s) -C n (t, s) = t t0 s t0 Φ L (t, u) cov U Xn u , U Xn v - cov U Xn-1 u , U Xn-1 v Φ L (s, v) T du dv.
We take the infinite matrix norm of both sides of this equality and use the upperbounds

Φ L (t, u) ∞ ≤ e L ∞ (T -t0) = k L and Φ L (t, u) T ∞ ≤ e L T ∞ (T -t0) = k L T (see appendix B) to obtain C n+1 (t, s) -C n (t, s) ∞ ≤ k L k L T t t0 s t0 cov U Xn u , U Xn v -cov U Xn-1 u , U Xn-1 v v ∞ du dv = k L k L T t t0 s t0 cov U Xn u • 1, U Xn v • 1 -cov U Xn-1 u • 1, U Xn-1 v • 1 v ∞ du dv. ( 26
)
According to equations (24) we are led to consider the difference A n -A n-1 , where:

A n def = S αβ C n ββ (u, u)C n ββ (v, v) -C n ββ (u, v) 2 C n ββ (u, u) x+ C n ββ (u, v) C n ββ (u, u) y+µ n β (v) S αβ y C n ββ (u, u) + µ n β (u) def = S αβ f n β (u, v)x + g n β (u, v)y + µ n β (v) S αβ h n β (u)y + µ n β (u) .

INRIA

We write next:

A n -A n-1 = S αβ f n β (u, v)x + g n β (u, v)y + µ n β (v) S αβ h n β (u)y + µ n β (u) -S αβ h n-1 β (u)y + µ n-1 β (u) + S αβ h n-1 β (u)y + µ n-1 β (u) S αβ f n β (u, v)x + g n β (u, v)y + µ n β (v) -S αβ f n-1 β (u, v)x + g n-1 β (u, v)y + µ n-1 β (v) .
The mean value theorem yields:

| A n -A n-1 |≤ S αβ ∞ S ′ αβ ∞ | x | | f n β (u, v) -f n-1 β (u, v) | + | y | | g n β (u, v) -g n-1 β (u, v) | + | µ n β (v) -µ n-1 β (v) | + | y | | h n β (u) -h n-1 β (u) | + | µ n β (u) -µ n-1 β (u) | .
Using the fact that

∞ -∞ | x | Dx = 2 π , we obtain: C n+1 (t, s) -C n (t, s) ∞ ≤ k L k L T k C 2 π t t0 s t0 f n (u, v) -f n-1 (u, v) ∞ dudv+ 2 π t t0 s t0 g n (u, v) -g n-1 (u, v) ∞ dudv+ (t -t 0 ) s t0 µ n (v) -µ n-1 (v) ∞ dv + (s -t 0 ) t t0 µ n (u) -µ n-1 (u) ∞ du+ 2 π (s -t 0 ) t t0 h n (u) -h n-1 (u) ∞ du ,
where

k C = max α β σ 2 αβ S αβ ∞ S ′ αβ ∞ . (27)
A similar process applied to the mean values yields:

µ n+1 (t) -µ n (t) ∞ ≤ k L µ t t0 h n (u) -h n-1 (u) ∞ du+ t t0 µ n (u) -µ n-1 (u) ∞ du .
We now use the mean value theorem and lemmas 3.5 and 3.4 to find upperbounds for

f n (u, v) -f n-1 (u, v) ∞ , g n (u, v) -g n-1 (u, v) ∞ and h n (u) -h n-1 (u) ∞ . We have |h n β (u) -h n-1 β (u)| = C n ββ (u, u) -C n-1 ββ (u, u) ≤ 1 2 √ k 0 C n ββ (u, u) -C n-1 ββ (u, u) ,
where k 0 is defined in lemma 3.5. Hence:

h n (u) -h n-1 (u) ∞ ≤ 1 2 √ k 0 C n (u, u) -C n-1 (u, u) ∞ .
Along the same lines we can show easily that:

g n (u, v) -g n-1 (u, v) ∞ ≤ k C n (u, v) -C n-1 (u, v) ∞ + C n (u, u) -C n-1 (u, u) ∞ ,
and that:

f n (u, v) -f n-1 (u, v) ∞ ≤ k θ(u, v) C n (u, v) -C n-1 (u, v) ∞ + C n (u, u) -C n-1 (u, u) ∞ + C n (v, v) -C n-1 (v, v) ∞ ,
where θ(u, v) is defined in lemma 3.4. Grouping terms together and using the fact that all integrated functions are positive, we write:

C n+1 (t, s) -C n (t, s) ∞ ≤ k [t0,t∨s] 2 1 θ(u, v) C n (u, v) -C n-1 (u, v) ∞ dudv+ [t0,t∨s] 2 1 θ(u, v) C n (u, u) -C n-1 (u, u) ∞ dudv+ [t0,t∨s] 2 C n (u, v) -C n-1 (u, v) ∞ dudv+ [t0,t∨s] 2 C n (u, u) -C n-1 (u, u) ∞ dudv+ [t0,t∨s] 2 µ n (u) -µ n-1 (u) ∞ dudv . ( 28
)
INRIA Note that, because of lemma 3.5, all integrals are well-defined. Regarding the mean functions, we write:

µ n+1 (t) -µ n (t) ∞ ≤ k [t0,t∨s] 2 C n (u, u) -C n-1 (u, u) ∞ dudv+ [t0,t∨s] 2 µ n (u) -µ n-1 (u) ∞ dudv . ( 29
)
Proceeding recursively until we reach C 0 and µ 0 we obtain an upperbound for C n+1 (t, s) -C n (t, s) ∞ (respectively for µ n+1 (t)µ n (t) ∞ ) which is the sum of less than 5 n terms each one being the product of k raised to a power less than or equal to n, times 2µ or 2Σ (upperbounds for the norms of the mean vector and the covariance matrix), times a 2n-dimensional integral I n given by

[t0,t∨s] 2 ρ 1 (u 1 , v 1 ) [t0,u1∨v1] 2 • • • [t0,un-2∨vn-2] 2 ρ n-1 (u n-1 , v n-1 ) [t0,un-1∨vn-1] 2 ρ n (u n , v n )du n dv n du n-1 dv n-1 • • • du 1 dv 1 ,
where the functions ρ i (u i , v i ), i = 1, • • • , n are either equal to 1 or to 1/ θ(u i , v i ). According to lemma 3.6, this integral is of the order of some positive constant raised to the power n divided by (n -1)!. Hence the sum is less than some positive constant k raised to the power n divided by (n -1)!. By taking the supremum with respect to t and s in [t 0 , T ] we obtain the same result for

C n+1 -C n ∞ (respectively for µ n+1 -µ n ∞ ). Since the series n≥1 k n n! is convergent, this implies that C n+p -C n ∞ (respectively µ n+p -µ n ∞
) can be made arbitrarily small for large n and p and the sequence C n (respectively µ n ) is a Cauchy sequence.

We can now prove the following theorem Theorem 3.8. For any nondegenerate kP -dimensional Gaussian random variable Z 0 and any initial process X such that X(t 0 ) = Z 0 , the map F k has a unique fixed point in

M + 1 (C([t 0 , T ], Ê kP )) towards which the sequence {F (n) k (X)} ∞ n=0 of Gaussian processes converges in law. Proof. Since C([t 0 , T ], Ê kP ) (respectively C([t 0 , T ] 2 , Ê kP ×kP )
) is a Banach space for the uniform norm, the Cauchy sequence µ n (respectively C n ) of proposition 3.7 converges to an element µ of

C([t 0 , T ], Ê kP ) (respectively an element C of C([t 0 , T ] 2 , Ê kP ×kP ))
. Therefore, according to theorem 3.1, the sequence {F Hence we know that there there exists at least one fixed point for the map F k . Assume there exist two distinct fixed points Y 1 and Y 2 of F k with mean functions µ i and covariance functions C i , i = 1, 2, with the same initial condition. Since for all n ≥ 1 we have

F (n) k (Y i ) = Y i , i = 1, 2, the proof of proposition 3.7 shows that µ n 1 -µ n 2 ∞ (respectively C n 1 -C 2 n ∞ ) is upper- bounded by the product of a positive number a n (respectively b n ) with µ 1 -µ 2 ∞ (respectively with C 1 -C 2 ∞ ). Since lim n→∞ a n = lim n→∞ b n = 0 and µ n i = µ i , i = 1, 2 (respectively C n i = C i , i = 1, 2)
, this shows that µ 1 = µ 2 and C 1 = C 2 , hence the two Gaussian processes Y 1 and Y 2 are indistinguishable.

Conclusion

We have proved that for any non degenerate initial condition Z 0 there exists a unique solution of the mean field equations. The proof of theorem 3.8 is constructive, and hence provides a way for computing the solution of the mean field equations by iterating the map F k , starting from any initial process X satisfying X(t 0 ) = Z 0 , for instance a Gaussian process such as an Ornstein-Uhlenbeck process. We build upon these facts in section 5.

Note that the existence and uniqueness is true whatever the initial time t 0 and the final time T .

Existence and uniqueness of stationary solutions

So far, we have investigated the existence and uniqueness of solutions of the mean field equation for a given initial condition. We are now interested in investigating stationary solutions, which allow for some simplifications of the formalism.

A stationary solution is a solution whose probability distribution does not change under the flow of the equation. These solutions have been already investigated by several authors (see [START_REF] Brunel | Fast global oscillations in networks of integrate-and-fire neurons with low firing rates[END_REF][START_REF] Sompolinsky | Chaos in Random Neural Networks[END_REF]). We propose a new framework to study and simulate these processes. Indeed we show in this section that under a certain contraction condition there exists a unique solution to the stationary problem. As in the previous section our proof is constructive and provides a way to simulate the solutions.

Remark.

The long-time mean field description of a network is still a great endeavor in mathematics and statistical physics. In this section we formally take the mean field equation we obtained and let t 0 → -∞. This way we obtain an equation which is the limit of the mean field equation when t 0 → -∞. It means that we consider first the limit N → ∞ and then t 0 → -∞. These two limits do not necessarily commute and there are known examples, for instance in spin glasses, where they do not.

It is clear that in order to get stationary solutions, we have to assume that the leak matrix L(t) does not depend upon t. Therefore, the resolvent Φ L (t, s) is equal to e L(t-s) . To ensure stability of the solutions and the existence of a stationary process we also assume that the real parts of its eigenvalues are negative:

Re(λ) < -λ L λ L > 0 (30)
for all eigenvalues λ of L. This implies that we only consider first-order system since otherwise the matrix L has eigenvalues equal to 0.

INRIA

For the same reason, we assume that the noise matrix Λ(t) and the input currents I(t) are constant in time. We further assume that the matrix Λ has full rank. Proposition 4.1. Under the previous assumptions we have

       lim t0→-∞ e L(t-t0) = 0, t -∞ e L(t-s) ds = ∞ 0 e Lu ∞ du def = M L < ∞, t -∞ e L T (t-s) ∞ ds = ∞ 0 e L T u ∞ du def = M L T < ∞,
and the process Y t0 t = t t0 e L(t-s) Λ•dW s is well-defined, Gaussian and stationary when t 0 → -∞. Proof. The first property follows from the fact that Re(λ) < -λ L for all eigenvalues λ of L. This assumption also implies that there exists a norm on Ê P such that e Lt ≤ e -λLt ∀t ≥ 0,

and hence e Lt ∞ ≤ ke -λLt ∀t ≥ 0, (31)
for some positive constant k. This implies the remaining two properties.

The stochastic integral t t0 e L(t-s) Λ • dW s is well-defined ∀t ≤ T and is Gaussian with zeromean. We note Y t0 t the corresponding process. Its covariance matrix reads:

Σ Y t 0 t Y t 0 t ′ = t∧t ′ t0
e L(t-s) ΛΛ T e L T (t ′ -s) ds.

Let us assume for instance that t ′ < t and perform the change of variable u = ts to obtain

Σ Y t 0 t Y t 0 t ′ = t-t0 t-t ′ e Lu ΛΛ T e L T u du e L T (t ′ -t) .
Under the previous assumptions this matrix integral is defined when t 0 → -∞ (dominated convergence theorem) and we have

Σ Y -∞ t Y -∞ t ′ = +∞ t-t ′ e Lu ΛΛ T e L T u du e L T (t ′ -t) , which is a function of t ′ -t.
This guarantees that there exists a stationary distribution of the equation:

dX 0 (t) = L • X 0 (t) dt + Λ • dW t , (32)
such that E [X 0 (t)] = 0. We have

X 0 (t) = t -∞ e L(t-s) Λ • dW s . Its covariance matrix Σ 0 is equal to Σ Y -∞ t Y -∞ t
and is independent of t.

We call long term mean field equation (LTMFE) the implicit equation:

V(t) = t -∞ e L(t-s) U V s • 1 + I ds + X 0 (t) (33)
where X 0 is the stationary process defined by equation ( 32) and where U V (t) is the effective interaction process introduced previously.

We next define the long term function

F stat : M + 1 (C((-∞, T ], Ê P ) → M + 1 (C((-∞, T ], Ê P ): F stat (X) t = t -∞ e L(t-s) U X s • 1 + I ds + X 0 (t).
Proposition 4.2. The function F stat is well defined on M + 1 (C((-∞, T ], Ê P ).

Proof. We have already seen that the process X 0 is well defined. The term -s) ds I is also well defined because of the assumptions on L.

t -∞ e L(t-s) I ds = t -∞ e L(t
Let X be a given process in M + 1 (C((-∞, T ], Ê P ). To prove the proposition we just have to ensure that the Gaussian process t -∞ e L(t-s) U X s • 1 ds is well defined. This results from the contraction assumption on L and the fact that the functions S αβ are bounded. We decompose this process into a "long memory" term 0 -∞ e L(t-s) U X s • 1 ds and the interaction term from time t = 0, namely t 0 e L(t-s) U X s • 1 ds. This latter term is clearly well defined. We show that the memory term is also well defined as a Gaussian random variable.

We write this term e Lt 0

-∞ e -Ls U X s • 1 ds and consider the second factor. This random variable is Gaussian, its mean reads

∞ 0 e Ls µ U X -s • 1 ds where µ U X -s =   P β=1 Jαβ E [S αβ (X β (-s))] + I α   α=1...P
The integral defining the mean is well-defined because of ( 31) and the fact that the functions S αβ are bounded. A similar reasoning shows that the corresponding covariance matrix is well-defined.

Hence the Gaussian process 0 -∞ e -Ls U X s • 1ds is well defined, and hence for any process X ∈ M + 1 (C((-∞, T ], Ê P ), the process F stat (X) is well defined.

We can now prove the following proposition.

Proposition 4.3. The mean vectors and the covariance matrices of the processes in the image of F stat are bounded.

INRIA

Proof. Indeed, since E [X 0 (t)] = 0, we have:

E [F stat (X) t ] ∞ = t -∞ e L(t-s) µ U X s ds ∞ ≤ M L (µ + I ∞ ) def = µ LT .
In a similar fashion the covariance matrices of the processes in the image of F stat are bounded. Indeed we have:

E F stat (X) t F stat (X) T t = Σ 0 + t -∞ t -∞ e L(t-s1) diag   β σ 2 αβ E [S αβ (X β (s 1 ))S αβ (X β (s 2 ))]   e L T (t-s2) ds 1 ds 2 ,
resulting in

E F stat (X) t F stat (X) T t ∞ ≤ Σ 0 ∞ + k 2 σ max λ L 2 def = Σ LT .
Lemma 4.4. The set of stationary processes is invariant by F stat .

Proof. Since the processes in the image of F stat are Gaussian processes, one just needs to check that the mean of the process is constant in time and that its covariance matrix C(s, t) only depends on ts.

Let Z be a stationary process and Y = F stat (Z). We denote by µ Z α the mean of the process Z α (t) and by C Z α (ts) its covariance function. The mean of the process U Z αβ reads:

m Z α,β (t) = E [S αβ (Z β (t))] = 1 2πC Z β (0) Ê S αβ (x)e (x-µ Z β ) 2 2C Z β (0) dx
and hence does not depends on time. We note µ Z the mean vector of the stationary process U Z • 1.

Similarly, its covariance function reads:

∆ Z αβ (t, s) = E [S αβ (Z β (t))S αβ (Z β (s))] = Ê 2 S αβ (x)S αβ (y) exp - 1 2 x -µ Z β y -µ Z β T C Z β (0) C Z β (t -s) C Z β (t -s) C Z β (0) -1 x -µ Z β y -µ Z β dx dy
which is clearly a function, noted ∆ Z αβ (ts), of ts. Hence U Z • 1 is stationary and we denote by C U Z (ts) its covariance function.

It follows that the mean of Y t reads:

µ Y (t) = E [F stat (Z) t ] = E [X 0 (t)] + E t -∞ e L(t-s) I + U Z s • 1 ds = t -∞ e L(t-s) I + E U Z s • 1 ds = 0 -∞ e Lu du I + µ Z Since we proved that E U Z s • 1 = µ Z
was not a function of s. Similarly, we compute the covariance function and check that it can be written as a function of (ts). Indeed, it reads:

C Y (t, s) = t -∞ s -∞ e L(t-u) Cov(U Z u • 1, U Z v • 1)e L T (s-v) du dv + Cov(X 0 (t), X 0 (s)) = 0 -∞ 0 -∞ e Lu C U Z (t -s + (u -v))e L T v du dv + C X0 (t -s) since the process X 0 is stationary. C Y (t, s) is clearly a function of t -s.
Hence Y is a stationary process, and the proposition is proved.

Theorem 4.5. The sequence of processes

{F (n) stat (X)} ∞ n=0 is uniformly tight.
Proof. The proof is essentially the same as the proof of theorem 3.3, since we can write

F stat (X) t = e Lt F stat (X) 0 + t 0 e L(t-s) (U X s • 1 + I) ds + t 0 e L(t-u) ΛdW s
F stat (X) t appears as the sum of the random variable F stat (X) 0 and the Gaussian process defined by

t 0 e L(t-s) (U X s • 1 + I) ds + t 0 e L(t-u) ΛdW s which is equal to F k (X) t defined in section 3 for t 0 = 0. Therefore F (n) stat (X) t = F (n) k (X) t for t > 0.
We have proved the uniform tightness of the sequence of processes

{F (n) k (X)} ∞
n=0 in theorem 3.3. Hence, according to Kolmogorov's criterion for tightness, we just have to prove that the sequence of Gaussian random variables:

F (n) stat (X) 0 = 0 -∞ Φ L (-u)(U F (n) stat (X) u • 1 + I)du + X 0 (0)
n≥0 is uniformly tight. Since it is a sequence of Gaussian random variables, it is sufficient to prove that their means and covariance matrices are upperbounded to obtain that for any ε > 0 there exists a compact K ε such that for any n ∈ AE, we have È(F (n) stat (X) 0 ∈ K ε ) ≥ 1ε. This is a consequence INRIA of proposition 4.3 for the first random variable and of the definition of X 0 for the second. By Kolmogorov's criterion the sequence of processes

{F (n) stat (X)} ∞ n=0 is uniformly tight
In order to apply theorem 3.1 we need to prove that the sequences of covariance and mean functions are convergent. Unlike the case of t 0 finite, this is not always true. Indeed, to ensure existence and uniqueness of solutions in the stationary case, the parameters of the system have to satisfy a contraction condition, and proposition 3.7 extends as follows.

Proposition 4.6. If λ L defined in (30) satisfies the conditions (34) defined in the proof, depending upon k C (defined in ( 27)), k 0 , µ LT and Σ LT (defined in proposition 4.3)then the sequences of covariance matrix functions C n (t, s) and of mean functions µ n (t), s, t in [t 0 , T ] are Cauchy sequences for the uniform norms.

Proof. The proof follows that of proposition 3.7 with a few modifications that we indicate. In establishing the equation corresponding to (26) we use the fact that Φ L (t, u) ∞ ≤ ke -λL(t-u) for some positive constant k and all u, t, u ≤ t. We therefore have:

C n+1 (t, s) -C n (t, s) ∞ ≤ k 2 e -λL(t+s) t -∞ s -∞ e λL(u+v) Cov U Xn u , U Xn v -Cov U Xn-1 u , U Xn-1 v v ∞ du dv
The rest of the proof proceeds the same way as in proposition 3.7. Equations ( 28) and ( 29) become:

C n+1 (t, s) -C n (t, s) ∞ ≤ Ke -λL(t+s) [-∞,t∨s] 2 e λL(u+v) f (u, v) C n (u, v) -C n-1 (u, v) ∞ dudv+ [-∞,t∨s] 2 e λL(u+v) f (u, v) C n (u, u) -C n-1 (u, u) ∞ dudv+ [-∞,t∨s] 2 e λL(u+v) C n (u, v) -C n-1 (u, v) ∞ dudv+ [-∞,t∨s] 2 e λL(u+v) C n (u, u) -C n-1 (u, u) ∞ dudv+ [-∞,t∨s] 2 e λL(u+v) µ n (u) -µ n-1 (u) ∞ dudv , RR n°6454
and

µ n+1 (t) -µ n (t) ∞ ≤ Ke -λL(t+s) [-∞,t∨s] 2 e λL(u+v) C n (u, u) -C n-1 (u, u) ∞ dudv+ [-∞,t∨s] 2 e λL(u+v) µ n (u) -µ n-1 (u) ∞ dudv ,
for some positive constant K, function of k, k C (defined in ( 27)), and k 0 .

Proceeding recursively until we reach C 0 and µ 0 we obtain an upperbound for C n+1 (t, s) -C n (t, s) ∞ (respectively for µ n+1 (t)µ n (t) ∞ ) which is the sum of less than 5 n terms each one being the product of K n , times 2µ LT or 2Σ LT , times a 2n-dimensional integral

I n given [-∞,t∨s] 2 ρ 1 (u 1 , v 1 ) [-∞,u1∨v1] 2 • • • [-∞,un-2∨vn-2] 2 ρ n-1 (u n-1 , v n-1 ) [-∞,un-1∨vn-1] 2 e λL(un+vn) ρ n (u n , v n )du n dv n du n-1 dv n-1 • • • du 1 dv 1 , where the functions ρ i (u i , v i ), i = 1, • • • , n are either equal to 1 or to 1/ θ(u i , v i ).
It can be shown by straightforward calculation that each sub-integral contributes at most either

K 0 λ 2 L if ρ i = 1 or π 2 K 0 λ 3/2 L ,
in the other case. Hence we obtain factors of the type

K n 0 1 λ 2 L p π 2 1 λ 3/2 L n-p = π 2 n-p 1 λ L (3n+p)/2 K n 0 , where 0 ≤ p ≤ n. If λ L < 1, (λ L ) (3n+p)/2 ≥ λ 2n L and else (λ L ) (3n+p)/2 ≥ λ 3n/2 L . Since π 2 n-p ≤ π 2
n we obtain the two conditions

1 > λ 2 L ≥ 5 π 2 KK 0 or λ 3/2 L > 5 π 2 KK 0 and λ L ≥ 1 (34)
Putting all these results together we obtain the following theorem of existence and uniqueness of solutions for the long term mean field equations: Theorem 4.7. Under the contraction conditions (34), the function F stat has a unique solution in M + 1 (C((-∞, T ], Ê P ) which is stationary, and for any process X, the sequence {F (n) stat (X)} ∞ n=0 of Gaussian processes converges in law toward the unique fixed point of the function F stat .

INRIA

Proof. The proof is essentially similar to the one of theorem 3.8. Indeed, the mean and the covariance matrixes converge since they are Cauchy sequences in the complete space of continuous functions equipped with the uniform norm. Using theorem 3.1, we obtain that the sequence converges to a process Y which is necessarily a fixed point of F stat . Hence we have existence of a fixed point for F stat . The uniqueness comes from the results obtained in the proof of proposition 4.6. The limiting process is necessarily stationary. Indeed, let X be a stationary process. Then for any n ∈ AE, the process F (n) stat (X) will be stationary by the virtue of lemma 4.4, and hence so will be the limiting process which is the only fixed point of F stat .

Hence in the stationary case, the existence and uniqueness of a solution is not always ensured. For instance if the leaks are too small (i.e. when the time constants of the decay of the membrane potentials are too long) then the sequence can diverge or have multiple fixed points.

Numerical experiments

Simulation algorithm

Beyond the mathematical results, the framework that we introduced in the previous sections gives us a strategy to compute numerically the solutions of the dynamic mean-field equations. Indeed, we proved in section 3 that under very moderate assumptions on the covariance matrix of the noise, the iterations of the map F k starting from any initial condition converge to the solution of the mean field equations.

This convergence result gives us a direct way to compute numerically the solution of the mean field equations. Since we are dealing with Gaussian processes, determining the law of the iterates of the map F k amounts to computing its mean and covariance functions. In this section we describe our numerical algorithm in the case of the Simple Model of section 2.2.2.

Computing

F k .
Let X be a P -dimensional Gaussian process of mean µ X = (µ X α (t)) α=1...P and covariance C X = (C X αβ (s, t)) α,β∈{1...P } . We fix a time interval [t 0 = 0, T ] and denote by Y the image of the process X under F 1 . In the case of the simple model, the covariance of Y is diagonal. Hence in this case the expressions we obtain in section 3 simply read:

µ Y α (t) = µ X α (0)e -t/τα + t 0 e -(t-s)/τα ( P β=1 Jαβ E [S α,β (X β (s))] + I α (s))ds = µ X α (0)e -t/τα + t 0 e -(t-s)/τα I α (s)ds + P β=1 Jαβ t 0 e -(t-s)/τα +∞ -∞ S αβ x v X β (s) + µ X β (s) Dxds.
where we denoted v X α (s) the standard deviation of X α at time s, instead of C X αα (s, s). Thus, knowing v X α (s), s ∈ [0, t] we can compute µ Y α (t) using a standard discretization scheme of the integral, with a small time step compared with τ α and the characteristic time of variation of the input current I α . Alternatively, we can use the fact that µ Y α satisfies the differential equation:

dµ Y α dt = - µ Y α τ α + P β=1 Jαβ +∞ -∞ S αβ x v X β (t) + µ X β (t) Dx + I α (t),
and compute the solution using a Runge-Kutta algorithm (which is faster and more accurate). Note that, when all the standard deviations of the process X are null for all time t ∈ [0, T ], we obtain a standard dynamical system. Nevertheless, in the general case, v X β (t) > 0 for some βs, and the dynamical evolution of µ Y α depends on the Gaussian fluctuations of the field X. These fluctuations must be computed via the complete equation of the covariance diagonal coefficient C Y αα (t, s), which reads:

C Y αα (t, s) = e -(t+s)/τα v X α (0) + τ α s 2 α 2 e 2s τα -1 + P β=1 σ 2 αβ t 0 s 0 e (u+v)/τα ∆ X αβ (u, v)dudv ,
where:

∆ X αβ (u, v) = IR 2 S αβ   x v X β (u)v X β (v) -C X ββ (u, v) 2 v X β (v) + y C X ββ (u, v) v X β (v) + µ X β (u)   × S αβ y v X β (v) + µ X β (v) Dx Dy.
Unless if we assume the stationarity of the process (see e.g. section 5.2), this equation cannot be written as an ordinary differential equation. We clearly observe here the non-Markovian nature of the problem: C X αα (t, s) depends on the whole past of the process until time t ∨ s.

This covariance can be split into the sum of two terms: the external noise contribution C OU αα (t, s) = e -(t+s)/τα v X α (0) +

ταs 2 α 2 e 2s
τα -1 and the interaction between the neurons. The external noise contribution is a simple function and can be computed directly. To compute the interactions contribution to the standard deviation we have to compute the symmetric two-variables function:

H X αβ (t, s) = e -(t+s)/τα t 0 s 0 e (u+v)/τα ∆ X αβ (u, v)dudv,
from which one obtains the standard deviation using the formula

C Y αα (t, s) = C OU αα (t, s) + P β=1 σ 2 αβ H X αβ (t, s).

INRIA

To compute the function H X αβ (t, s), we start from t = 0 and s = 0, where H X αβ (0, 0) = 0. We only compute H X αβ (t, s) for t > s because of the symmetry. It is straightforward to see that:

H X αβ (t + dt, s) = H X αβ (t, s) 1 - dt τ α + D X αβ (t, s)dt + o(dt), with D X αβ (t, s) = e -s/τα s 0 e v/τα ∆ X αβ (t, v)dv.
Hence computing H X αβ (t + dt, s) knowing H X αβ (t, s) amounts to computing D αβ (t, s). Fix t ≥ 0. We have D αβ (t, 0) = 0 and

D X αβ (t, s + ds) = D X αβ (t, s)(1 - ds τ α ) + ∆ X αβ (t, s)ds + o(ds).
This algorithm enables us to compute H X αβ (t, s) for t > s. We deduce H X αβ (t, s) for t < s using the symmetry of this function. Finally, to get the values of H X αβ (t, s) for t = s, we use the symmetry property of this function and get:

H X αβ (t + dt, t + dt) = H X αβ (t, t) 1 - 2dt τ α + 2D X αβ (t, t)dt + o(dt).
These numerical schemes provide an efficient way for computing the mean and the covariance functions of the Gaussian process F 1 (X) (hence its probability distribution) knowing the law of the Gaussian process X. The algorithm used to compute the solution of the mean field equations for the general models GM1 and GMk is a straightforward generalization.

Analysis of the algorithm

Convergence rate As proved in theorem 3.8, given Z 0 a nondegenerate kP -dimensional Gaussian random variable and X a Gaussian process such that X(0) = Z 0 , the sequences of means and covariance functions computed theoretically converge uniformly towards those of the unique fixed point of the map F k . It is clear that our algorithm converges uniformly towards the real function it emulates. Hence for a finite N , the algorithm will converge uniformly towards the mean and covariance matrix of the process F N k (X). Denote by X f the fixed point of F k in M + 1 (C([t 0 , T ], Ê kP )), of mean µ X f (t) and covariance matrix C X f (t, s), and by F N k (X) the numerical approximation of F N k (X) computed using the algorithm previously described, whose mean is noted (t, s). The uniform error between the simulated mean after N iterations with a time step dt and the fixed point's mean and covariance is the sum of the numerical error of the algorithm and the distance between the simulated process and the fixed point, is controlled by:

µ d F N k (X) (t) and whose covariance matrix is noted C d F N k (X)
µ d F N k (X) -µ X f ∞ + C d F N k (X) -C X f ∞ = O( (N + T ) dt + R N (k max ) ) (35)
where k max = max(k, k) and k and k) are the constants that appear in the proof of proposition 3.7 for the mean and covariance functions, and R N (x) is the exponential remainder, i.e. R N (x) = ∞ n=N x n /n!. Indeed, we have:

µ d F N k (X) -µ X f ∞ ≤ µ d F N k (X) -µ F N k (X) ∞ + µ F N k (X) -µ X f ∞ (36)
The discretization algorithm used converges in O(dt). Let us denote by C 1 the convergence constant, which depends on the sharpness of the function we approximate, which can be uniformly controlled over the iterations. Iterating the numerical algorithm has the effect of propagating the errors. Using these simple remarks we can bound the first term of the righthand side of (36), i.e. the approximation error at the N th iteration:

µ d F N k (X) -µ F N k (X) ∞ ≤ C 1 N dt
Because the sequence of means is a Cauchy sequence, we can also bound the second term of the righthand side of (36):

µ F N k (X) -µ X f ∞ ≤ ∞ n=N µ F n+1 k (X) -µ F n k (X) ∞ ≤ ∞ n=N k n n! =: R N (k)
for some positive constant k introduced in the proof of proposition 3.7. The remainders sequence (R n (k)) n≥0 converges fast towards 0 (an estimation of its convergence can be obtained using the fact that lim sup k→∞ (1/k!) 1/k = 0 by Stirling's formula).

Hence we have:

µ d F N k (X) -µ X f ∞ ≤ C 1 N dt + R N (k) (37)
For the covariance, the principle of the approximation is exactly the same:

C d F N k (X) -C X f ∞ ≤ C d F N k (X) -C F N k (X) ∞ + C F N k (X) -C X f ∞
The second term of the righthand side can be controlled using the same evaluation by R N ( k)

where k is the constant introduced in the proof of proposition 3.7, and the first term is controlled by the rate of convergence of the approximation of the double integral, which is bounded by C 2 (N + T ) dt where C 2 depends on the parameters of the system and the discretization algorithm used.

Hence we have:

C d F N k (X) -C X f ∞ ≤ C 2 (N + T -t 0 ) dt + R N ( k) (38)
The expressions (37) and ( 38) are the sum of two terms, one of which is increasing with N and T and decreasing with dt and the other one decreasing in N . If we want to obtain an estimation with an error bounded by some ε > 0, we can for instance fix N such that max(R N (k), R N ( k)) < ε/2 and then fix the time step dt smaller than min( ε/(2C 1 N ), ε/(2C 2 (N + Tt 0 )) ).

INRIA

Complexity

The complexity of the algorithm depends on the complexity of the computations of the integrals. The algorithm described hence has the complexity O(N ( T dt ) 2 ).

5.2

The importance of the covariance: Simple Model, one population.

As a first example and a benchmark for our numerical scheme we revisit the work of Sompolinsky and coworkers [START_REF] Sompolinsky | Chaos in Random Neural Networks[END_REF]. These authors studied the case of the simple model with one population (P = 1), with the centered sigmoidal function S(x) = tanh(gx), centered connectivity weights J = 0 of standard deviation σ = 1 and no input (I = 0, Λ = 0). Note therefore that there is no "noise" in the system, which therefore does not match the non degeneracy conditions of proposition 3.4 and of theorem 3.8 . This issue is discussed below. In this case, the mean equals 0 for all t. Nevertheless, the Gaussian process is non trivial as revealed by the study of the covariance C(t, s).

Stationary solutions

Assuming that the solution of the mean field equation is a stationary solution with C(t, s) ≡ C(ts) = C(τ ), Sompolinsky and his collaborators found that the covariance obeyed a second order differential equation :

d 2 C dτ 2 = - ∂V q ∂C . (39)
This form corresponds to the motion of a particle in a potential well and it is easy to draw the phase portrait of the corresponding dynamical system. However, there is a difficulty. The potential V q depends on a parameter q which is in fact precisely the covariance at τ = 0 (q = C(0)). In the stationary case, this covariance depends on the whole solution, and hence cannot be really considered as a parameter of the system. This is one of the main difficulties in this approach: mean field equations in the stationary regime are self-consistent.

Nevertheless, the study of the shape of V q , considering q as a free parameter gives us some informations. Indeed, V q has the following Taylor expansion (V q is even because S is odd):

V q (C) = λ 2 C 2 + γ 4 C 4 + O(C 6 )
where λ = (1g 2 J 2 S ′ 2 q ) and γ = 1 6 J 2 g 6 S (3) 2 q ), φ q being the average value of φ under the Gaussian distribution with mean zero and variance q = C(0).

If λ > 0, i.e. when g 2 J 2 S ′ 2 q < 1, then the dynamical system (39) has a unique solution C(t) = 0, ∀t ≥ 0. This corresponds to a stable fixed point (i.e. a deterministic trajectory, µ = 0 with no fluctuations) for the neural network dynamics. On the other hand, if g 2 J 2 S ′ 2 q ≥ 1 there is a homoclinic trajectory in (39) connecting the point q = C * > 0 where V q vanishes to the point C = 0. This solution is interpreted by the authors as a chaotic solution in the neural network. A stability analysis shows that this is the only stable 3 stationary solution [START_REF] Sompolinsky | Chaos in Random Neural Networks[END_REF].

The equation for the homoclinic solution is easily found using energy conservation and the fact that V q (q) = 0 and dVq dC (q) = 0. One finds:

u = dC dx = --V q (C).
At the fourth order in the Taylor expansion of V q this gives

C(τ ) = -2λ γ cosh( -λ 2 τ )
.

Though λ depends on q it can be used as a free parameter for interpolating the curve of C(τ ) obtained from numerical data.

Numerical experiments

This case is a good benchmark for our numerical procedure since we know analytically the solutions we are searching for. We expect to find two regimes. In one case the correlation function is identically zero in the stationary regime, for sufficiently small g values or for a sufficiently small q (trivial case). The other case corresponds to a regime where C(τ ) > 0 and C(τ ) → 0 has τ → +∞ ("chaotic" case). This regime requires that g be sufficiently large and that q be large too. We took τ α = 0.25, σ αα = 1. For these values, the change in dynamics predicted by Sompolinsky and collaborators is g c = 4.

In sections 3 and 4 we have introduced the assumption of non-degeneracy of the noise, in order to ensure that the mean field process was non degenerate. However, in the present example, there is no external noise in the evolution, so we can observe the effects of relaxing this hypothesis in a situation where the results of proposition 3.4 and of theorem 3.8 cannot be applied. First, we observed numerically that, without external noise, the process could become degenerate (namely some eigenvalues of the covariance matrix C α (t, s) become very small and even vanish.). This has also an incidence on the convergence of the method which presents numerical instabilities, though the iterations leads to a curve which is well fitted by the theoretical results of Sompolinsky et al. (see Fig. 3) . The instability essentially disappears if one adds a small noise. But, note that in this case, the solution does not match with Sompolinsky et al. theoretical calculation (see Fig. 3).

Modulo this remark, we have first considered the trivial case corresponding to small g values. We took g = 0.5 and T = 5. We choose as initial process the stationary Ornstein-Uhlenbeck process corresponding to the uncoupled system with Λ = 0.1. We drew µ α (0) randomly from the uniform distribution in [-1, 1] and v α (0) randomly from the uniform distribution in [0, 1].

Starting from this initial stationary process, we iterated the function F 1 . Then, during the iterations, we set s α = 0 in order to match the conditions imposed by Sompolinsky and coworkers. We observe that the method converges towards the expected solution: the mean function converges to zero, while the variance v(t) decreases exponentially fast in time towards a constant value corresponding to the stationary regime. This asymptotic value decreases between two consecutive iterations, which is consistent with the theoretical expectation that v(t) = 0 in the stationary regime of INRIA Figure 3: Numerical solution of the mean field equation after 14 iterations in the chaotic case (g = 5). We clearly see the numerical instabilities in the no-noise case, which do not exist in the low-noise case.

the trivial case. Finally, we observe that the covariance C(ts, s) stabilizes to a curve that does not depend on s and the stationary value (large ts) converges to zero.

We applied the same procedure for g = 5 corresponding to the "chaotic" regime. The behavior was the same for µ(t) but was quite different for the covariance function C(t, s). Indeed, while in the first case the stationary value of v(t) tends to zero with the number of iterations, in the chaotic case it stabilizes to a finite value. In the same way, the covariance C(ts, s) stabilizes to a curve that does not depend on s. The shape of this curve can be extrapolated thanks to Sompolinsky et al. results. We observe a very good agreement with the theoretical predictions with a fit f 4 (x) = a cosh(b(x-δ)) , corresponding to the fourth expansion of V q . Using a 6-th order expansion of

V q (x) = a 2 x 2 + b 4 x 4 + c 6 x 2 gives a fit f 6 (x) = ρ cosh(λ(x-δ)) 1 q 1+K 2 - 1 cosh 2 (λ(x-δ))
, where ρ, K, λ are explicit functions of a, b, c, we obtain a slightly better approximation.

Mean field equations for two populations with a negative feedback loop.

Let us now present a case where the fluctuations of the Gaussian field act on the dynamics of µ α (t) in a non trivial way, with a behavior strongly departing from the naive mean field picture. We consider two interacting populations where the connectivity weights are Gaussian random variables J αβ ≡ N ( Jαβ , σ αβ = 1) for (α, β) ∈ {1, 2} 2 . We set S αβ (x) = tanh(gx) and I α = 0, s α = 0, α = 1, 2.

This opens up the possibility to have an instability of the fixed point (µ = 0) leading to a regime where the average value of the membrane potential oscillates. This occurs if J11 + J22 > 0 and if g is larger than:

g c = 2 τ ( J11 + J22 )
.

The corresponding bifurcation is a Hopf bifurcation. The situation is different if one takes into account the fluctuations of the Gaussian field. Indeed, in this case the stability of the fixed point µ = 0 depends on v(t). More precisely, the real and imaginary part of the eigenvalues of DG(0, v(t)) depend on v(t). Therefore, the variations of v(t) act on the stability and oscillations period of v(t). Though the evolution of µ(t), v(t) are coupled we cannot consider this evolution as a coupled dynamical system, since v(t) = C(t, t) is determined by the mean field equation for C(t, s) which cannot be written as an ordinary differential equation. Note that we cannot assume stationarity here, as in the previous case, since µ(t) depends on time for sufficiently large g. This opens up the possibility of having complex dynamical regimes when g is large.

Numerical experiments

We have considered the case J11 = J22 = 5,τ = 0.1 giving a Hopf bifurcation for g c = 2 when J = 0 (fig. 4). The trajectory of µ 1 (t) and v 1 (t) is represented in Figure 4 in the case g = 3. When J = 0, µ 1 (t) presents regular oscillations (with non linear effects since g = 3 is larger than the critical value for the Hopf bifurcation, g c = 2). In this case, the solution v 1 (t) = 0 is stable as seen on the figure. When J = 0 the Gaussian field has (small) fluctuations which nevertheless strongly interact with the dynamics of µ 1 (t), leading to a regime where µ 1 (t) and v 1 (t) oscillate periodically

Conclusion

The problem of bridging scales is overwhelming in general when studying complex systems and in particular in neuroscience. After many others we look at this difficult problem from the theoretical and numerical viewpoints, hoping to get closer to its solution from relatively simple and physically/biologically plausible first principles and assumptions. One of our motivations is to better understand such phenomenological neural mass models as that of [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF].

We consider several populations of neurons and start from a microscopic, i.e. individual, description of the dynamics of the membrane potential of each neuron that contains four terms. The first one controls the intrinsic dynamics of the neuron. It is linear in this article but this assumption is not essential and could probably be safely removed if necessary. The second term is a stochastic input current, correlated or uncorrelated. The third one is a deterministic input current, and the fourth one describes the interaction between the neurons through random connectivity coefficients that weigh the contributions of other neurons through a set of functions that are applied to their membranes potentials. The only hypothesis on these functions is that they are smooth and bounded. The obvious Figure 4: Evolution of the mean µ 1 (t) and variance v 1 (t) for the mean field of population 1, for J = 0 and J = 2, over a time window [0,20]. n is the number of iterations of F 1 defined in section 3. This corresponds to a number of iterations for which the method has essentially converged (up to some precision). Note that v 1 (t) has been magnified by a factor of 100. Though Gaussian fluctuations are small, they have a strong influence on µ 1 (t). choice of sigmoids is motivated by standard rate models ideas. Another appealing choice is a smooth approximation to a Dirac delta function thereby opening a window on the world of spiking neurons.

We then derive the mean field equations and provide a constructive and new proof, under some mild assumptions, of the existence and uniqueness of a solution of these equations over finite and infinite time intervals. The key idea is to look at this mean field description as a global problem on the probability distribution of the membranes potentials, unlike previous studies. Our proof provides an efficient way of computing this solution and our numerical experiments show a good agreement with previous studies.

In the case where the nonlinearities are chosen to be sigmoidal our results shed a new and fascinating light on existing neural mass models. Indeed these appear as approximations of the mean field equations where the intricate but fundamental coupling between the time variations of the mean membrane potentials and their fluctuations, as represented by the covariance functions, is neglected. This article is just a small step toward answering from the theoretical and numerical standpoints the questions raised by this coupling but we are convinced that a host of interesting results can be found there.

INRIA

A Identification of the mean field equations

Ben-Arous and Guionnet studied from a mathematical point of view the problem of finding a meanfield description of large networks of spin glasses. They obtained using different methods of stochastic analysis a weak limit of the law of a given spin and proved their independence.

Our equations do not directly fit in their study: indeed, the spin intrinsic dynamics is nonlinear while the interaction is linear, and everything in done in dimension one. Nevertheless, their proof extends to our case which is somehow more simple. For instance in the case of the Simple Model with one population, we can readily adapt their proof in our case. More precisely, let P = 1, the equation of the network reads:

τ dV j t = (-V j t + N i=1 J ij S(V i t )) dt + σdW j t
In this case, we define for X ∈ M + 1 (C([t 0 , T ], Ê) the effective interaction term (U X t ) which is the effective interaction process defined in 2.1, i.e. the Gaussian process of mean Jαβ E [S(X t )] and of covariance: Cov U X t , U X s =: σ 2 αβ E [S(X t )S(X s )]. Let us note P the law of the membrane potential when there is no interaction (it is an Ornstein-Ulhenbeck process), and the empirical measure V N = 1 N N i=1 δ V i . We can prove that under the probability distribution averaged over the connectivities, see below, the empirical measure satisfies a large deviation principle with good rate function H defined as in [START_REF] Guionnet | Averaged and quenched propagation of chaos for spin glass dynamics[END_REF]. Using this large deviation result, we can prove annealed and quenched tightness of the empirical measure, and finally its convergence towards the unique process where the good rate function H achieves its unique minimum, which is defined by the property of having a density with respect to P and whose density satisfies the implicit equation:

Q ≪ P dQ dP = E exp T 0 U Q t dW t - 1 2 T 0 (U Q t ) 2 dt ( 40
)
where E denotes the expectation over the effective interaction process U Q .

Proof. Pathwise (strong) uniqueness of solution directly comes from the results on the SDE with Lipschitz coefficients (see e.g. [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF], Theorem 2.5 of Chapter 5)). It is clear that X t0 = X 0 . We use Itô's formula for the product of two stochastic processes to prove that the process ( 49) is solution of equation ( 48):

dX t = L(t)Φ L (t,

C Matrix norms

In this section we recall some definitions on matrix and vector norms. Let M n×n be the set of n × n real matrices. It is a vector space of dimension n 2 and the usual L p norms 1 ≤ p ≤ ∞ can be defined. Given L ∈ M n×n , we note L v p the corresponding norm. Given a vector norm, noted , on Ê n the induced norm, noted , on M n×n is defined as

L = sup x∈Ê n , x ≤1 Lx x
Since M n×n is finite dimensional all norms are equivalent. In this article we use the following norms

(i). L ∞ = max i n j=1 |L ij |. (ii). L v ∞ = max i, j |L ij | (iii). L 2 = sup x∈Ê n , x 2 ≤1 Lx 2 x 2
. This so-called spectral norm is equal to the square root of the largest singular value of L which is the largest eigenvalue of the positive matrix L T L. If L is positive definite this is its largest eigenvalue which is also called its spectral radius, noted ρ(L).

D Important quantities

Table 1 summarizes some notations which are introduced in the article and used in several places.

INRIA

Figure 1

 1 Figure 1: a. Neural mass model: a population of pyramidal cells interacts with two populations of inter-neurons: an excitatory one and an inhibitory one. b. Block representation of the model. The h boxes account for the synaptic integration between neuronal populations. S boxes simulate cell bodies of neurons by transforming the membrane potential of a population into an output firing rate. The coefficient J αβ is the random synaptic efficiency of population β on population α (P is the pyramidal population, E the excitatory and I the inhibitory ones), and the constants (C i ) model the mean strength of the synaptic connections between populations (it is the mean of the Js).

,

  Because of the properties of the stochastic integral, b a

Figure 2 :

 2 Figure 2: The change of coordinates.

  n=0 of Gaussian processes converges in law toward the Gaussian process Y with mean function µ and covariance function C. This process is clearly a fixed point of F k .

  t 0 )X 0 + Φ L (t, t)I(t) + L (s, t)Λ(s)dW s dt = L(t) Φ L (t, t 0 )X 0 +

		t
		L(t)Φ L (t, s)I(s) ds dt
		t0
	t
	+ Φ L (t, t)Λ(t)dW t +	L(t)Φ
	t0

t t0 Φ L (s, t)I(s) ds + t t0 Φ L (s, t)Λ(s)dW s + I(t) dt + Λ(t)dW t = (L(t)X(t) + I(t)) dt + Λ(t)dW t

Hence the theorem is proved.

RR n°6454

Note that in the cases we treat in this paper, the matrix J (N) is always almost surely invertible since it has non-degenerate Gaussian coefficients, and hence the equivalence in law will always be valid INRIA

In(Bogachev, 1998, Chapter

3.8), the property is stated whenever the mean and covariance are defined on a separable Hilbert space.RR n°6454

More precisely, this is the only minimum for the large deviation functional.RR n°6454

Aknowledgement

This research was partly supported by funding of the European Union under the grant no. 15879 (FACETS) and the Fondation d'Entreprise EADS.

Theoretical framework.

The dynamic mean field equation for µ α (t) is given, in differential form, by:

Let us denote by G α (µ, v(t)) the function in the righthand side of the equality. Since S is odd, ∞ -∞ S( v β (t)x)Dx = 0. Therefore, we have G α (0, v(t)) = 0 whatever v(t), and hence the point µ 1 = 0, µ 2 = 0 is always a fixed point of this equation.

Let us study the stability of this fixed point. To this purpose, we compute the partial derivatives of G α (µ, v(t)) with respect to µ β for (α, β) ∈ {1, 2} 2 . We have:

and hence at the point µ 1 = 0, µ 2 = 0, these derivatives read:

In the case v α (0) = 0, J = 0, s α = 0, implying v α (t) = 0, t ≥ 0, the equation for µ α reduces to:

which is the standard Amari-Cohen-Grossberg-Hopfield system. This corresponds to the naive mean field approach where Gaussian fluctuations are neglected. In this case the stability of the fixed point µ = 0 is given by the sign of the largest eigenvalue of the Jacobian matrix of the system that reads:

-

For the sake of simplicity we assume that the two time constants τ α are equal and we denote this value τ . The eigenvalues are in this case -1 τ + gλ, where λ are the eigenvalues of J and have the form:

Hence, they are complex whenever J12 J21 < -( J11 -J22 ) 2 /4, corresponding to a negative feedback loop between population 1 and 2. Moreover, they have a real part only if J11 + J22 is non zero (self interaction).

INRIA

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

We can also prove following the steps of Ben-Arous and Guionnet in (Ben-Arous and [START_REF] Guionnet | Averaged and quenched propagation of chaos for spin glass dynamics[END_REF] that there exists a unique solution to this equation, and that this solution satisfies the nonlinear nonmarkovian stochastic differential equation:

which can also be written as our mean field equation, averaged on the connectivities (see (Ben-Arous and [START_REF] Ben-Arous | Large deviations for Langevin spin glass dynamics[END_REF]). More precisely, let L V be the law of the solution of the equation:

which is exactly equation ( 13). They prove that V satisfies the nonlinear equation:

This result is likely extendable to the multi population case but the corresponding mathematical developments are out of the scope of this paper.

B The resolvent

In this appendix we introduce and give some useful properties of the resolvent Φ L of a homogeneous differential equation dx dt

where L : [t 0 , T ] → M P ×P (or (-∞, T ] → M P ×P ) is C 0 .

Definition B.1. The resolvent of ( 42) is defined as the unique solution of the linear equation:

where Id P is the P × P identity matrix.

Proposition B.1. The resolvent satisfies the following properties:

INRIA (iii). Let be a norm on M P ×P and assume that L(t) ≤ k L on [t 0 , T ]. Then we have:

Similarly, if L T (t) ≤ k L T on [t 0 , T ] we have:

(iv). We have

TrL(s) ds Proof. The properties (i) and (ii) are directly linked with the property of group of the flow of a reversible ODE. (iii) is an application of Gronwald's lemma. (iv) is obtained by a first order Taylor series expansion.

We also need in the article a lower bound on Φ L (t, t 0 ) for all t ∈ [t 0 , T ] in the general case where L is not constant. This can be achieved for example using Floquet's theory. Consider the interval [t 0 , 2Tt 0 ] and define the continuous periodic function L(t) of period 2(Tt 0 ) defined by

The corresponding resolvent Φ L(t, t 0 ) is equal to Φ L (t, t 0 ) for t 0 ≤ t ≤ T . Φ L(2Tt 0 , t 0 ) is invertible and hence there exists a ∈ Ê such that e 2a(T -t0) < |λ| for all eigenvalues λ of Φ L(2Tt 0 , t 0 ). One of Floquet's theorems states that there exists a norm on Ê P and γ > 0 such that γe a(t-t0) < Φ L(t, t 0 ) t ≥ t 0 , and in particular γe a(t-t0) < Φ L (t, t 0 ) t 0 ≤ t ≤ T (47)

Theorem B.2 (Solution of an inhomogeneous linear SDE). The solution of the inhomogeneous linear Stochastic Differential Equation:

dX t = (L(t)X(t) + I(t)) dt + Λ(s)dW s X t0 = X 0 (48) can be written using the resolvent: