How Gibbs distributions may naturally arise from synaptic adaptation mechanisms. A model-based argumentation

This paper addresses two questions in the context of neuronal networks dynamics, using methods from dynamical systems theory and statistical physics: (i) How to characterize the statistical properties of sequences of action potentials ("spike trains") produced by neuronal networks ? and; (ii) what are the effects of synaptic plasticity on these statistics ? We introduce a framework in which spike trains are associated to a coding of membrane potential trajectories, and actually, constitute a symbolic coding in important explicit examples (the so-called gIF models). On this basis, we use the thermodynamic formalism from ergodic theory to show how Gibbs distributions are natural probability measures to describe the statistics of spike trains, given the empirical averages of prescribed quantities. As a second result, we show that Gibbs distributions naturally arise when considering"slow"synaptic plasticity rules where the characteristic time for synapse adaptation is quite longer than the characteristic time for neurons dynamics.

Spike trains as a "neural code". Neurons activity results from complex and nonlinear mechanisms (Hodgkin & Huxley, 1952;Cronin, 1987;Dayan & Abbott, 2001;Gerstner & Kistler, 2002b), leading to a wide variety of dynamical behaviours (Cronin, 1987;Guckenheimer & Labouriau, 1993). This activity is revealed by the emission of action potentials or "spikes". While the shape of an action potential is essentially always the same for a given neuron, the succession of spikes emitted by this neuron can have a wide variety of patterns (isolated spikes, periodic spiking, bursting, tonic spiking, tonic bursting, etc . . . ) (Izhikevich, 2004;Brette & Gerstner, 2005;Touboul, 2008), depending on physiological parameters, but also on excitations coming either from other neurons or from external inputs. Thus, it seems natural to consider spikes as "information quanta" or "bits" and to seek the information exchanged by neurons in the structure of spike trains. Doing this, one switches from the description of neurons in terms of membrane potential dynamics, to a description in terms of spike trains. This point of view is used, in In this paper we suggest that these two questions are closely entangled and must both be addressed in the same setting. In order to have a good control on the mathematics we mainly consider the case of neural networks models where one has a full characterization of the generic dynamics (Cessac, 2008;Cessac & Viéville, 2008). Thus, our aim is not to provide general statements about biological neural networks. We simply want to have a good mathematical control of what is going on in specific models, with the hope that this analysis should shed some light on what happens (or does not happen) in "real world" neural systems. However, though part of the results used here are rigorous, this work relies also on "working assumptions" that we have not been able to check rigorously. These assumptions, which are essentially used to apply the standard theorems in ergodic theory and thermodynamic formalism, provide a logical chain which drives us to important conclusions that could be checked in experimental data.
The paper is organised as follows. In section 2 we introduce a framework in which spike trains are associated to a coding of membrane potential trajectories, and actually, constitute a symbolic coding in explicit examples. On this basis, we show how Gibbs distributions are natural probability measures to describe the statistics of spike trains, given the data of known empirical averages. Several authors have discussed the relevance of Gibbs distribution in the context of Hopfield model (see (Amit, 1989) and references therein) and more recently to spike trains (Wood, Roth, & Black, 2006;Kang & Amari, 2008). A breakthrough has been made in Tkacik et al., 2006). Actually, our approach has been greatly influenced by these two papers, though our results hold in a wider context. In section 3, we discuss the effect of synaptic adaptation when considering "slow" synaptic plasticity rules where the characteristic time for synapse adaptation is quite a bit longer that the characteristic time for neurons dynamics. These rules are formulated in the context of thermodynamic formalism, where we introduce a functional, closely related to thermodynamic potentials like free energy in statistical physics, and called "topological pressure" in ergodic theory. In this setting we show that the variations of synaptic weights leads to variation of the topological pressure that can be smooth ("regular periods"), or singular ("phase transitions"). Phase transitions are in particular associated to a change of "grammar" inducing modifications in the set of spike trains that the dynamics is able to produce. We exhibit a functional, closely related to the topological pressure, that decreases during regular periods. As a consequence, when the synaptic weights converge to a fixed value, this functional reaches a minimum. This minimum corresponds to a situation where spike trains statistics are characterized by a Gibbs distribution, whose potential can be explicitly written. An example with numerical computations is presented in section 4.

Neural state.
We consider a set of N neurons. Each neuron i is characterized by its state, X i , which belongs to some compact set I ∈ IR M . M is the number of variables characterizing the state of one neuron (we assume that all neurons are described by the same number of variables). A typical example is an integrate and fire model where M = 1 and X i = V i is the membrane potential of neuron i and I = [V min , V max ] (see section 2.5). Other examples are provided by conductances based models of Hodgkin-Huxley type 1 (Hodgkin & Huxley, 1952). Then X i = (V i , m i , n i , h i ) where m i , n i are respectively the activation variable for Sodium and Potassium channels and h i is the inactivation variable for the Sodium channel.
The evolution of these N neurons is given by a deterministic dynamical system of form: represents the dynamical state of a network of N neurons at time t, while time is discrete (for a discussion on time discretisation in spiking neural networks see (Cessac & Viéville, 2008)). Thus X ∈ M = I N where M is the phase space of (1), and F γ (M) ⊂ M. The map F γ : M → M depends on a set of parameters γ ∈ IR P . The typical case considered here is γ = W, I (ext) where W is the matrix of synaptic weights and I (ext) is some external current, assumed to be independent of time in the present paper (see section 2.5 for two explicit examples). Thus γ is a point in a P = N 2 + N dimensional space of control parameters.

Natural partition.
Neurons are excitable systems. Namely, neuron i "fires" (emits a spike or action potential), whenever its state X i belongs to some connected region P 1 of its phase space. Otherwise, it is quiescent (X ∈ P 0 = I \ P 1 ). In Integrate and Fire models neuron i fires whenever its membrane potential V i exceeds some threshold θ. In this case, the corresponding region is P 1 = [θ, V max ]. In Fitzhugh-Nagumo (FitzHugh, 1955, 1961Nagumo, Arimoto, & Yoshizawa, 1962) or Hodgkin-Huxley model (Hodgkin & Huxley, 1952), the firing corresponds to the crossing of a manifold called the threshold separatrix (Cronin, 1987) which separates the phase space into two connected regions P 0 and P 1 . For N identical neurons this leads to a "natural partition" P of the product phase space M. Call Λ = {0, 1} N , ω = (ω i ) N i=1 ∈ Λ. Then, P = {P ω } ω∈Λ , where P ω = P ω1 × P ω2 × · · · × P ωN . Equivalently, if X ∈ P ω , all neurons such that ω i = 1 are firing while neurons such that ω k = 0 are quiescent. We call therefore ω a "spiking pattern".

Raster plots.
To each initial condition X ∈ M we associate a "raster plot" ω = {ω(t)} +∞ t=0 such that X(t) ∈ P ω(t) , ∀t ≥ 0. We write X ⇀ ω. Thus, ω is the sequence of spiking patterns displayed by the neural network when prepared with the initial condition X. On the other way round, we say that an infinite sequence ω = {ω(t)} +∞ t=0 is an admissible raster plot if there exists X ∈ M such that X ⇀ ω. We call Σ γ the set of admissible raster plots for the set of parameters γ.
Pω i , ω being a spiking pattern. In this way, one associates naturally to an orbit of (1) a raster plot.

5
The dynamics (1) induce a dynamics on the set of raster plot in the following way: In this construction σ γ , called the "left shift", shifts the raster plot left-wise at each time step of the dynamics. Thus, in some sense, raster plots provide a code for the orbits of (1). But, the correspondence may not be one-to-one. That is why we use the notation ⇀ instead of →.
Let us introduce the following notation. If we are interested in a prescribed sequence of spiking patterns ω(s), . . . ω(t), from time s to time t, we denote by [ω] s,t , the set of raster plots whose spiking patterns from time s to time t match this sequence (cylinder set).

Asymptotic dynamics.
We are here mainly interested in the asymptotic behavior of (1) and set a few notations and notions.
The ω-limit set, Ω, is the set of accumulation points of F t γ (M). Since M is closed and invariant, we have Ω = ∞ t=0 F t γ (M). In dissipative systems (i.e. a volume element in the phase space is dynamically contracted), the ω-limit set typically contains the attractors of the system. A compact set A ∈ M is called an attractor for F γ if there exists a neighborhood U of A and a time n > 0 such that F n γ (U) ⊂ U In all examples considered in the present paper F γ is dissipative and the phase space is divided into finitely many attraction basins each of them containing an attractor (see Fig. 2). Simple examples of attractors are stable fixed points, or stable periodic orbits. More complex attractors such as chaotic attractors can be encountered as well. The attraction basins and attractors change when the parameters γ vary. These changes can be smooth (structural stability) or sharp (typically this arises at bifurcations points).
2.5 A generic example: generalized integrate and fire models 2.5.1 Quiescent stage of the neuron.
Among the various models of neural dynamics, generalized integrate and fire models play a central role, due to their (relative) simplicity, while the biological plausibility is well accepted (Rudolph & Destexhe, 2006;Gerstner & Kistler, 2002a). A representative class of such models is of form: defining the evolution of the membrane potential V k of neuron k. Here τ L = RC ≃ 10 − 20ms is the membrane time-constant related to the membrane resistance and its electric capacity, while E L ≃ −80mV is the related reversal potential. The term i (ext) k is an external "current" 2 assumed to be time constant in this paper. In the end of this section we shall however consider the case where some noise is superimposed is the n-th firing time 3 of neuron j and t (n) j t is the list of firing times of all neurons up to time t.
The synaptic currents reads: where E ± are reversal potential (typically E + ≃ 0mV and E − ≃ −75mV ). E and I refers respectively to excitatory and inhibitory neurons and the + (−) sign is relative to excitatory (inhibitory) synapses. Note that conductances are always positive thus the sign of the post-synaptic potential (PSP) is determined by the reversal potentials E ± . At rest (V k ∼ −70mV ) the + term leads to a positive PSP while − leads to a negative PSP. Conductances depend on past spikes via the relation: In this equation, M j (t, V) is the number of times neuron j has fired at time t and t (n) j t is the list of firing times of all neurons up to time t. G kj are positive constants, proportional to the synaptic efficacy: We use the convention W kj = 0 if there is no synapse from j to k. Finally, α represents the unweighted shape (called a α-shape) of the post-synaptic potentials. Known examples of α-shapes are where H is the Heaviside function.
Then, we may write (2) in the form (see (Cessac & Viéville, 2008) for more details): and: The previous equations hold whenever the neuron is quiescent, i.e. whenever membrane potential is smaller than a threshold θ > 0, usually depending on time (to account for characteristics such as refractory period of the neuron) and on the neuronal state. Here, however, we assume that θ is a constant (in fact without loss of generality (Gerstner & Kistler, 2002a)). When the membrane potential exceeds the threshold value, the neuron "fires" (emission of an action potential or "spike"). The spike shape depends on the model. In the present case, the membrane potential is reset instantaneously to a fixed value V reset ≃ E L , corresponding to the value of the membrane potential when the neuron is at rest. For simplicity we set V reset = 0 without loss of generality. This modeling of the spike introduces a natural notion of spike time, but the price to pay is to introduce a discontinuity in dynamics. Moreover, since spike is instantaneous the set of possible spike times occurring within time interval is uncountable in the continuous case. Introducing a minimal time discretisation at scale δ removes this pathology.

Time discretisation.
Assuming that spike times are only known within a precision δ > 0 (Cessac & Viéville, 2008) one shows that the dynamics of membrane potential, discretized at the time scale δ writes: where: is the number of spikes emitted by neuron j up to time t, in the raster plot ω. In the sequel we assume that α ± (u) = 0 whenever |u| ≥ τ M , i.e. g k (t, [ω] 0,t ) depends on past firing times over a finite time horizon τ M . We also set δ = 1.

Model I and model II.
A step further, two additional simplifications can be introduced: 1. Assuming that ρ k (t, t+ δ, [ω] 0,t ) = ρ is almost constant, which is equivalent to considering "current" synapses instead of "conductance" synapses, i.e. neglect the dependency of i (syn) k with respect to V k ; 2. Considering a much simpler form for J k (t, [ω] 0,t ), where each connected neuron simply increments the membrane potential during its firing state. This is equivalent to considering the post-synaptic profiles in the previous equations as instantaneous step-wise profiles.
Considering these assumptions leads the following dynamics of membrane potential: where F γ,i is the i-th component of F γ , Z(x) = 1 if x ≥ θ and 0 otherwise, both the integrate and firing regime being integrated in this equation. It turns out that this corresponds exactly to the time discretisation of the standard integrate and fire neuron model, which as discussed in e.g., (Izhikevich, 2003) provides a rough but realistic approximation of biological neurons behaviors. This model is called model I in the sequel whereas the model based on (4) is called model II.
2.5.5 Genericity results for models I and II.
The map F γ in models I and II is locally contracting (Cessac, 2008;Cessac & Viéville, 2008), but it is not smooth due to the sharp threshold in neurons firing definition. The singularity set where F γ is not continuous is: This is the set of membrane potential vectors such that at least one of the neurons has a membrane potential exactly equal to the threshold. Because of this, dynamics exhibit sensitivity to perturbations for orbits which approach too close the singularity set. Now, call Then, the following theorem holds (Cessac, 2008;Cessac & Viéville, 2008).
Theorem 1 For a generic (in a metric and topological sense) set of parameters γ, d(Ω, S) > 0. Consequently, 1. Ω is composed of finitely many periodic orbits with a finite period, Thus, attractors are generically stable period orbits.
2. There is a one-to-one correspondence between membrane potential trajectories and raster plots, except for a negligible set of M.
3. There is a finite Markov partition.

Remarks
1. Note however that, depending on parameters (synaptic weights, external current), some periods can be quite large (well beyond any accessible computational time) Also, the number of stable periodic orbits can grow exponentially with the number of neurons.
2. Result 2 means that spike trains provide a symbolic coding for the dynamics.
3. See section 2.6.1 for the implications of result 3.
These models constitute therefore nice examples where one has a good control on dynamics and where one can apply the machinery of thermodynamic formalism (see Appendix). Note that since dynamics is eventually periodic, one may figure out that there is nothing non trivial to say from the respect of spike train statistics. The key point is that period are practically so large that there is any hope to see the periodicity in real experiments. Thus, one may do "as if" the system were "chaotic". This is a fortiori the case if one superimposes upon the external current a small amount of noise to the external current I (ext) .
Also, though dynamics of model I and II may look rather trivial compared to what is expected from biological neural networks, it might be that such models are able to approximate trajectories of a real neural network, for suitable values of γ (e.g. after a suitable adaptation of the synaptic weights) and provided N , the number of neurons, is sufficiently large. This type of properties are currently sought by computational neuroscientist with interesting indications that 'IF models are good enough" to approximate biological neurons spike trains (Jolivet, Rauch, Lescher, & Gerstner, 2006). See (Rostro-Gonzalez, ) for a recent illustration of this.

Spikes dynamics and statistics.
2.6.1 Grammar.
Though dynamics (1) produces raster plots, it is important to remark that it is not able to produce any possible sequence of spiking patterns. This fundamental fact is most often neglected in computational neuroscience literature and leads, e.g. to severe overestimation of the system's entropy. Also, it plays a central role in determining spike train statistics.
There are therefore allowed and forbidden sequences depending on conditions involving the detailed form of F γ and on the parameters γ (synaptic weights and external current). Let ω, ω ′ be two spiking patterns. The transition ω → ω ′ is legal or admissible if there exists a neural state X such that neurons fire according to the firing pattern ω and, at the next time, according to the firing pattern ω ′ . Equivalently, F γ (P ω ) ∩ P ω ′ = ∅. An admissible transition must satisfy compatibility conditions depending on synaptic weights and currents.
As an example, for model I, compatibility conditions write, for all i = 1 . . . N : Note that, while the two first conditions only depend on the partition element P ω , the two last ones depend on the point X ∈ P ω . Basically, this means that the natural partition is not a Markov partition.
The idea is then to find a refinement of the natural partition such that allowed transitions depend only on the partition-elements. A nice situation occurs when the refined partition is finite. Equivalently, compatibility conditions can be obtained by a finite set of rules or a finite grammar. This is possible if there exists some finite integer r with the following property: Construct blocks of spiking patterns ω(0) . . . ω(r − 1) of width r. There are at most 2 N r such possible blocks (and in general quite a bit less admissible blocks). Label each of these blocks with a symbol α ∈ A, where A is called an alphabet. Each block corresponds to a connected domain ∩ r−1 s=0 F −s γ [P ω(s) ] ⊂ M and to a cylinder set [ω] 0,r−1 ⊂ Σ γ . Define a transition matrix G γ : A × A → {0, 1}, depending on γ, with entries g αβ , such that g αβ = 1 if the transition β → α is admissible, and 0 otherwise. To alleviate the notations in the next equation write then all admissible raster plots are obtained via the finite transition matrix G γ which provides the grammar of spiking patterns sequences.
In models I, II such a finite grammar exists for the (generic) set of parameters γ such that d(Ω, S) > 0. Moreover, there are open domains in the space of parameters where the grammar is fixed and is not affected by small variations of synaptic weights or current.
For more general models, it might be that there is no finite grammar to describe all admissible raster plots. When dealing with experiments, one has anyway only access to finite raster plots and the grammar extracted from empirical observations has a finite number of rules (see (Collet, Galves, & Lopez, 1995;Chazottes, Floriani, & Lima, 1998;Chazottes, 1999) for nice applications of this idea in the field of turbulence). This empirical grammar is compatible with the dynamics at least for the time of observation.
The grammar depends on γ. Assume now that we are continuously varying γ. This corresponds to following some path in the space of control parameters. Then, this path generically crosses domains where grammar is fixed, while at the boundary of these domains, there is a grammar modification, and thus a change in the set of admissible raster plots that dynamics is able to produce. Synaptic adaptation, as described in section 3, corresponds precisely to this situation. Thus, we expect that the neural network exhibits "regular" periods of adaptation where spike trains before and after synaptic changes, correspond to the same grammar. Between these regular periods, sharp changes in raster plots structure are expected, corresponding somehow to displaying new transitions and new rules in the spike code. Actually, the ability of neural networks to display, after adaptation, a specific subset of admissible sequences, provided by a specific grammar, is an important issue of this paper.

Spike responses of neurons.
Neurons respond to excitations or stimuli by finite sequences of spikes. In model I and II a stimulus is typically the external current but more general forms (spike trains coming from external neurons) could be considered as well. Thus, the dynamical response R of a neuronal network to a stimuli S (which can be applied to several neurons in the network), is a sequence ω(t) . . . ω(t + n) of spiking patterns. "Reading the neural code" means that one seeks a correspondence between responses and stimuli. However, the spike response does not only depend on the stimulus, but also on the network dynamics and therefore fluctuates randomly (note that this apparent randomness is provided by the dynamical evolution and does not require the invocation of an exogenous noise). Thus, the spike response is sought as a conditional probability P (R|S) (Rieke et al., 1996). Reading the code consists in inferring P (S|R) e.g. via Bayesian approaches, providing a loose dictionary where the observation of a fixed spikes sequences R does not provide a unique possible stimulus, but a set of stimuli, with different probabilities. Having models for conditional probabilities P (R|S) is therefore of central importance. For this, we need a good notion of statistics.

Performing statistics.
These statistics can be obtained in two different ways. Either one repeats a large number of experiments, submitting the system to the same stimulus S, and computes P (R|S) by an average over these experiments. This approach relies on the assumption that the system has the same statistical properties during the whole set of experiments (i.e. the system has not evolved, adapted or undergone bifurcations meanwhile). In our setting this amounts to assuming that the parameters γ have not changed.
Or, one performs a time average. For example, to compute P (R|S), one counts the number of times n(R, T, ω) when the finite sequence of spiking patterns R, appears in a spike train ω of length T , when the network is submitted to a stimulus S. Then, the probability P (R|S) is estimated by: T .
This approach implicitly assumes that the system is in a stationary state. The empirical approach is often "in-between". One fixes a time window of length T to compute the time average and then performs an average over a finite number N of experiments corresponding to selecting different initial conditions. In any case the implicit assumptions are essentially impossible to control in real (biological) experiments, and difficult to prove in models. So, they are basically used as "working" assumptions.
To summarize, one observes, from N repetitions of the same experiment, N raster plots ω m , m = 1 . . . N on a finite time horizon of length T . From this, one computes experimental averages allowing to estimate P (R|S) or, more generally, to estimate the average value, φ , of some prescribed observable φ(ω). These averages are estimated by : Typical examples of such observables are φ(ω) = ω i (0) in which case φ is the firing rate of neuron i; φ(ω) = ω i (0)ω j (0) then φ measures the probability of spike coincidence for neuron j and i; φ(ω) = ω i (τ )ω j (0) then φ measures the probability of the event "neuron j fires and neuron i fires τ time step later" (or sooner according to the sign of τ ). In the same way P (R|S) is the average of the indicatrix function χ R (ω) = 1 if ω ∈ R and 0 otherwise. Note that in (7) we have used the shift σ t γ for the time evolution of the raster plot. This notation is more compact and more adapted to the next developments than the classical formula, reading, e.g., for firing rates 1 . This estimation depends on T and N . However, one expects that, as N , T → ∞, the empirical averageφ (N ,T ) converges to the theoretical average φ , as stated e.g. from the law of large numbers. Unfortunately, one usually does not have access to these limits, and one is lead to extrapolate theoretical averages from empirical estimations. The main difficulty is that these observed raster plots are produced by an underlying dynamics which is usually not explicitly known (as it is the case in experiments) or impossible to fully characterize (as it is the case in most large dimensional neural networks models). Even for models I and II, where one has a full characterization of generic orbits, a numerical generation of raster plots can produce periodic orbits whose period is out of reach. Thus, one is constrained to propose ad hoc statistical models. This amounts to assuming that the underlying dynamics satisfies specific constraints.
2.6.4 Inferring statistics from an hidden dynamics.
We follow this track, proposing a generic method to construct statistical models from a prescribed set of observables. For this, we shall assume that the observed raster plots are generated by a uniformly hyperbolic dynamical system, where spike trains constitute a symbolic coding for dynamics (see the appendix for more details). This a technical assumption allowing us to develop the thermodynamic formalism on a safe mathematical ground. Basically, we assume that the observed raster plots can be generated by a dynamical system which is chaotic with exponential correlation decay, and that a finite Markov partition exists, obtained by a refinement of the natural partition. This is obviously a questionable assumption, but let us give a few arguments in its favor. First, it has been argued and numerically checked in (Rostro-Gonzalez et al., 2009) that finite sequence of spiking patterns, produced by an hidden neural network, including data coming from biological experiments, can be exactly reproduced by a gIF model of type II, by adding hidden neurons and noise to the dynamics. This suggests that gIF models are somehow "dense" in the space of spiking neural networks dynamical systems, meaning that any finite piece of trajectory from a spiking neural network can be approached by a gIF model trajectory. Second, adding noise to model I, II makes them uniformly hyperbolic 4 , though not continuous. So, we expect the theory developed below to apply to model I, II. Now, dealing with spike train coming from a neural network whose mathematical properties are unknown (this is especially the case with biological neural networks), and whose statistical properties are sought, the idea is to do as if these data were generated by a dynamical system like model I or II.
As a matter of fact, the choice of a statistical model always relies on assumptions. Here we make an attempt to formulate these assumptions in a compact way with the widest range of application. These assumptions are compatible with the statistical models commonly used in the literature like Poisson models or Ising like modelsà la Schneidman and collaborators , but lead also us to propose more general forms of statistics. Moreover, our approach incorporates additional elements such as the consideration of neurons dynamics, and the grammar. This last issue is, according to us, fundamental, and, to the best of our knowledge, has never been considered before in this field. Finally, this postulate allows us to propose algorithms that can be applied to real data with an posteriori check of the initial hypotheses ).
On this basis we propose the following definition. Fix a set φ l , l = 1 . . . K, of observables, i.e. functions Σ γ → IR which associate real numbers to sequences of spiking patterns. Assume that the empirical average (7) of these functions has been computed, for a finite T and N , and thatφ l (T,N ) = C l .
A statistical model is a probability distribution ν on the set of raster plots such that: 1. ν(Σ γ ) = 1, i.e. the set of non admissible raster plots has a zero ν-probability.
2. ν is ergodic for the left-shift σ γ (see section 6.1 in the appendix for a definition).
Note that item 2 amounts to assuming that statistics are invariant under time translation. On practical grounds, this hypothesis can be relaxed using sliding time windows. This issue is discussed in more details in . Note also that ν depends on the parameters γ.
Assuming that ν is ergodic has the advantage that one does not have to average both over experiments and time. It is sufficient to focus on time average for a single raster plot, via the time-empirical average: defined within more details in the appendix, section 6.2.
Item 3 can be generalized as follows.
where d(µ, ν) is the relative entropy or Kullback-Leibler divergence between two measures µ, ν (see eq. (54) in the appendix.). Dealing with real or numerical data, it is obviously not possible to extrapolate to T → ∞, but the main idea here is to minimize the "distance" d(π (T ) ω , ν) between the empirical measure, coming from the data, and the statistical model. Especially, if several models can be proposed then, the "best" one minimizes the Kullback-Leibler divergence (see th. 2 in the Appendix). The main advantage is that d(π (T ) ω , ν) can be numerically estimated using the thermodynamic formalism. Note however that this approach may fail at phase transition points where d(µ, ν) = 0 does not necessarily imply µ = ν (Chazottes, 1999). Phase transition can be numerically detected from empirical data (Comets, 1997).
2.6.5 Gibbs measures as statistical models.
Variational principle. Statistical physics naturally proposes a canonical way to construct a statistical model: "Maximizing the entropy under the constraints ν(φ l ) = C l , l = 1 . . . K", (see (Jaynes, 1957) for a beautiful and deep presentation of what became, since then, a "folklore" result). In the context of thermodynamic formalism this amounts to solving the following variational principle (see eq. (43) in the Appendix).
where m (inv) is the set of invariant measures for σ and h the Kolomogorov-Sinai entropy or entropy rate (see Appendix for more details). The "potential" ψ is given by ψ = K l=1 λ l φ l where the λ l 's are adjustable Lagrange multipliers. A measure ν ψ which realizes the supremum, i.e.
is called, in this context, an "equilibrium state". The function P [ψ] is called the "topological pressure".
Gibbs states. Uniformly hyperbolic dynamical systems have equilibrium states. Moreover, in this case, equilibrium states are Gibbs states (and vice-versa). A Gibbs state, or Gibbs measure, is a probability measure such that, one can find some constants c 1 , c 2 with 0 < c 1 ≤ 1 ≤ c 2 such that for all n ≥ 1 and for all ω ∈ Σ: Basically, this means that the probability that a raster plot starts with the bloc [ω] 0,n−1 behaves like exp(S (n) ψ(ω)) Zn . One recognizes the classical Gibbs form where space translation in lattice system is replaced by time translation (shift σ t γ ) and where the normalization factor Z n is the partition function. Note that P [ψ] = lim sup n→∞ 1 n log Z n , so that P [ψ] is the formal analog of a thermodynamic potential (like free energy).
Topological pressure. The topological pressure is a convex function. Moreover, this is the generating function for the cumulants of the Gibbs distribution. Especially, the Lagrange multipliers λ l can be tuned to a value λ * l such that ν ψ [φ l ] = C l , (item 3) using: It expresses that the Gibbs state is given the tangent of the pressure at λ * l (Keller, 1998). Note that we have here assumed that the topological pressure is differentiable, namely that we are away from a phase transition. Higher order moments are obtained in the same way, especially second order moments related to the central limit theorem obeyed by Gibbs distributions (Bowen, 1975(Bowen, , 2008Keller, 1998). It is also possible to obtain averages of more complex functions than moments (Ji, 1989). The topological pressure is obtained via the spectrum of the Ruelle-Perron-Frobenius operator and can be calculated numerically when ψ has a finite (and small enough) range (see appendix for more details).
Validating a Gibbs statistical model. The Kullback-Leibler divergence provides some notion of "distance" between two measures. For µ an invariant measure and ν ψ a Gibbs measure with a potential ψ, both defined on the same set of sequences Σ, one has (see eq. (55) in the appendix): is not invariant, this suggests to use this relation to compare different statistical models (i.e. when several choices of observables φ l are possible) by choosing the one which minimizes the quantity: (see th. 2 in the appendix). The advantage is that this quantity can be numerically estimated ).
Probability of spiking patterns blocs. In this context, the probability of a spiking pattern block R = [ω] 0,n−1 of length n corresponding to the response R to a stimuli S "behaves like" (in the sense of eq. (51)): where the conditioning is made explicit in the dependence of λ * in the stimulus S. (Typically, referring to model I and II, S is an external current and is incorporated in γ.) This equation holds as well when S = 0 (no stimulus). Obviously, for two different stimuli the probability P (R|S) may drastically change. Indeed, different stimuli imply different parameters γ thus a different dynamics and different raster plots, with different statistical weights on a specific bloc R. The grammar itself can also change, in which case R may be forbidden for a stimulus and allowed for another one.

Canonical examples.
Firing rates. If φ l (ω) = ω l (0), then π (T ) ω (φ l ) = r l is the average firing rate of neuron l within the time period T . Then, the corresponding statistical model is a Bernoulli distribution where neuron l has a probability r l to fire at a given time. The probability that neuron l fires k times within a time delay n is a binomial distribution and the inter-spike interval is Poisson distributed (Gerstner & Kistler, 2002a).
where, here, the index l is an enumeration for all (non-ordered) pairs (i, j), then the corresponding statistical models has the form of an Ising model, as discussed by Schneidman and collaborators in Tkacik et al., 2006). As shown by these authors in experiments on the salamander retina, the probability of spike blocs estimated from the "Ising" statistical model fits quite better to empirical date than the classical Poisson model.
Enlarged spikes coincidence. As a generalization one may consider the probability of co-occurrence of spikes from neuron i and j within some time interval τ . The corresponding functions are φ l (ω) = ω i (0)ω j (τ ) and the probability of a spike bloc R writes: Further generalizations are considered below.

Conclusion.
The main conclusion of this section is that Gibbs measures constitute optimal (in the sense of entropy maximization) statistical models whenever a set of observable has been prescribed. It also allows to select a model between several choices by minimizing the Kullback-Leibler divergence between the empirical measure and the statistical models. In all the examples presented above the statistical model is determined by the choice of an a priori form for the potential. In the next section, we explicitly compute the potential in neural networks with synaptic adaptation. 16 3 Synaptic plasticity.
3.1 Synapses update as an integration over spikes trains.
Synaptic plasticity corresponds to the evolution of synaptic efficacy (synaptic weights). More precisely, in our notations (eq. (3)), W ij essentially provides the maximal amplitude of the post-synaptic potential induced, at the synapse linking j to i, when neuron j fires a spike. Synaptic weights evolve in time according to the spikes emitted by the pre-and post-synaptic neuron. In other words, the variation of W ij at time t is a function of the spiking sequences of neurons i and j from time t − T s to time t, where T s is time scale characterizing the width of the spike trains influencing the synaptic change.
In most examples the synapse update writes: . Thus, typically, synaptic adaptation results from an integration of spikes over the time scale T s .

Spikes responses of synapses.
The synaptic variation δW ij is the integrated response of the synapse from neuron j to neuron i when neuron j sends a spike sequence [ω j ] t−Ts,t and neuron i fires according to [ω i ] t−Ts,t . This response is not a deterministic function, while g is. Thus (11) is an approximation. As a matter of fact, the explicit form of g is usually derived from phenomenological considerations as well as experimental results where synaptic changes can be induced by specific simulations conditions, defined through the firing frequency of pre-and post-synaptic neurons (Bliss & Gardner-Medwin, 1973;Dudek & Bear, 1993), the membrane potential of the post-synaptic neuron (Artola, Bröcher, & Singer, 1990), or spike timing (Levy & Stewart, 1983;Markram et al., 1997;Bi & Poo, 2001) (see (Malenka & Nicoll, 1999) for a review). Thus, these results are usually based on a repetition of experiments involving the excitation of pre-and post-synaptic neurons by specific spike trains. The phenomenological plasticity rules derived from these experiments are therefore of statistical nature. Namely, they do not tell us what will be the exact changes induced on synapses when this or this spike train is applied to pre-and post-synaptic neuron. Instead, they provide us the average synaptic change. Thus, the function g(W ij , [ω i ] t−Ts,t , [ω j ] t−Ts,t ) in (11) is typically a statistical average of the synaptic response when the spike train of neuron j (resp. i) is [ω j ] t−Ts,t (resp. [ω i ] t−Ts,t ), and the actual synaptic weight value is W ij . In this paper we investigate the effects of those synaptic plasticity rules when the characteristic time scale T s is quite a bit larger than the time scale of evolution of the neurons. Namely, we consider slow adaptation rules. In this situation the synaptic weights update can be written as a function of the empirical average (42). We therefore consider adaptation rules of form: where ǫ is a parameter that will be typically small. φ ij (W ij , ω) is a function that we now make explicit.
Following (Gerstner & Kistler, 2002a) canonical form of adaptation rules use an expansion in singlet, pairs, triplets etc of spikes. Formally, one may write a generic form for these rules. Fix T s < ∞ a positive integer. Let L be the finite set of ordered lists L in {−T s , −T s + 1, . . . , 0} with the form L = (t 1 , t 2 , . . . , t L ) with −T s ≤ t 1 < t 2 < · · · < t L ≤ 0. For L ∈ L, 1 ≤ i ≤ N , we call an i-monomial a function m i,L (ω) = ω i (t 1 ) . . . ω i (t L ). We call an (i, j) polynomial a function: where h ijL1L2 (W ij ) are smooth functions of W ij . They can be constant, linear functions of W ij ("multiplicative" rules) or nonlinear functions allowing for example to constrain W ij within bounded values (e.g. "hard bound" or "soft bounds" rules). The form (13) is the most general form of synaptic adaptation rules considered in this paper, while g has the explicit form: Though explicit, this formulation is heavy, and we use instead the form (12), (13).

Examples of adaptation rule instantiation.
Hebbian learning uses firing rates. A typical example corresponds to Spike-Time Dependent Plasticity as derived from Bi and Poo (Bi & Poo, 2001) provides the average amount of synaptic variation given the delay between the pre-and post-synaptic spike. Thus, "classical" STDP writes (Gerstner & Kistler, 2002a;Izhikevich & Desai, 2003): with: where the shape of f has been obtained from statistical extrapolations of experimental data. Hence STDP is based on a second order statistics (spikes correlations). There is, in this case, an evident time scale T s = max(τ − , τ + ), beyond which f is essentially zero.
Generalized STDP As a last example, (Gerstner & Kistler, 2002a) propose a rule which corresponds to: (18) We capture thus the main standard synaptic adaptation rules with (13).

Coupled dynamics.
We consider now the following coupled dynamics. Neurons are evolving according to (1). We focus here on slow synapses dynamics. Namely, synaptic weights are constant for T ≥ T s consecutive dynamics steps, where T is large. This defines an "adaptation epoch". At the end of the adaptation epoch, synaptic weights are updated according to (11) Recall that γ = (W, I (ext) ) (see section 2.1) and γ (τ ) is the set of parameters at adaptation epoch τ . In the present setting the external current I (ext) is kept fixed and only synaptic weights are evolving. Basically, I (ext) is used as an external stimulus.

Statistical effects of synaptic plasticity.
Synaptic plasticity has several prominent effects. First, modifying the synaptic weights has an action on the dynamics resulting, either in smooth changes, or in sharp changes. In the first case, corresponding to parameters variations in a domain where the dynamical system is structurally stable, small variations of the synaptic weights induce smooth variations on the ω-limit set structure (e.g. points are slightly moved) and in the statistics of orbits. The grammar is not changed. On the opposite, when crossing bifurcations manifolds, the slightest change in one synaptic weight value results typically in a drastic reorganization of the ω-limit set structure where attractors and their attraction basin are modified (see Fig. 2 and (Daucé et al., 1998;Siri et al., 2007) for an illustration of this in the case of frequency rates neural networks with Hebbian learning). But this can also change the grammar and the set of admissible raster plots. Some forbidden transitions become allowed, some allowed transitions become forbidden. Finally, the spikes train statistics are also modified. Typically, the time-empirical average of the raster plot changes (π (T ) ω (τ ) → π (T ) ω (τ +1) ) and the corresponding statistical model also evolves.
But synaptic adaptation is not a mere variation of synaptic weights. Indeed, synaptic weights are associated with the coupled dynamics (19), meaning that synaptic changes depends on neurons dynamics, itself depending on parameters. Thus, a synaptic adaptation process corresponds to following a path in the space of parameters γ; this path is not determined a priori but evolves according to neurons dynamics. Though this path can belong to a unique domain of structural stability, this is not the case in general. At some point during the adaptation, this path crosses a bifurcation manifold inducing sharp changes in the dynamics and in the grammar 5 . Hence, after several changes of this type one can end up with a system displaying raster plots with a structure rather different from the initial situation. These changes depend obviously upon the detailed form of neuron dynamics (1) and upon the synaptic update mechanism (11); they are also conditioned by parameters such as stimuli. We now Analise these effects in the light of thermodynamic formalism and Gibbs distributions.

Static synaptic weights.
Let us first consider the situation where δW = 0, corresponding to synaptic weights that do not evolve. Typically, this is the case if synaptic weights matrix converges to an asymptotic value W * . From eq. (12) this corresponds to : This imposes a condition on the average value of φ ij . Therefore, from section 2.6.5 this imposes that the statistical model is a Gibbs measure ν with a potential of form: where The potential Φ in (21) is such that Φ(ω) = 0 if ω is admissible and Φ(ω) = −∞ if it is forbidden, so that forbidden raster plots have zero probability. This is a way to include the grammar in the potential (see appendix). The statistical parameters λ * ij , are given by eq. (8) in section 2.6.5, and, making φ ij explicit (eq. (13)): Since m i,L1 m j,L2 are monomials, this equation thus imposes constraints on the probability of spikes n-uplets ω i (t 1 ) . . . ω i (t L )ω j (s 1 ) . . . ω j (s L ). Note that condition (22) corresponds to an extremum for the topological pressure as a function of λ. Note also that the variational formulation of Gibbs (equilibrium) state writes, in this case P [ψ * ] = sup ν∈m (inv) h [ν] since ν(ψ) = 0. Hence, the corresponding Gibbs measure has maximal entropy.
Let us emphasize what we have obtained. The statistical model that fits with the condition (20) in section 2.6.5 is a Gibbs distribution such that the probability of a spin block R of length n is given by : where ψ * depends on S via the statistical parameters λ * and via Φ (grammar). When the situation δW = 0 corresponds to the asymptotic state for a synaptic adaptation process, this potential provides us the form of the statistical model after adaptation, and integrates all past changes in the synaptic weights. We now discuss this process within details.
3.6.1 Synaptic weights update and related potentials.

Smooth variations.
Let us first assume that these variations are smooth, i.e. one stays inside a domain where dynamics is structurally stable and the grammar G γ (τ ) is not modified. We assume moreover that the topological pressure is differentiable with respect to the variation δγ (τ ) .
Then, using eq. (49) in appendix, the adaptation rule (23) can be written in the form: Since, in this section, pressure P is assumed to be smooth, one has, using (49),(50): where κ (τ ) is the tensor with entries: j2 ] is the correlation function of φ i1,j1 , φ i2,j2 for the measure ν ψ (τ ) . Using the explicit form (13) of φ ij one can see that κ (τ ) is a sum of time correlations between uplets of spikes. This is a version of the fluctuation-dissipation theorem where the response to a smooth variation of the potential ψ (τ ) is given in terms of a series involving the time correlations of the perturbation (Ruelle, 1969(Ruelle, , 1999. This series converges provided that dynamics is uniformly hyperbolic. For sufficiently small ǫ, the matrix I + ǫ 2 κ (τ ) is positive and:

It follows that the variation
This variation is therefore negative when ǫ is sufficiently small.
We come therefore to the following important conclusion. The adaptation rule (12) is a gradient system where the function F (τ ) φ decreases when iterating synaptic adaptation rules. Were the transition τ → τ + 1 to be smooth for all τ , would F (τ ) φ reach a minimum 6 at some W * as τ → ∞. Such a minimum corresponds to ∇W * F (τ ) φ = 0, thus to δW = 0 according to eq. (25). Hence, this minimum corresponds to a static distribution for the synaptic weights. Therefore, according to section 3.5.1 the potential ψ (τ ) would converge to the potential (21) as τ → +∞. However, we cannot expect these transitions to be smooth for all τ .

Singular variations.
As we saw, during the synaptic adaptation process, the corresponding path in the space of parameters usually crosses bifurcations manifold inducing sharp changes in the dynamics. At those points pressure may not be smooth corresponding to a phase transition. These changes are not necessarily easy to detect numerically, though algorithms allowing to detect phase transitions from finite samples exist, based on rigorous results (Comets, 1997). Indeed, searching "at blind" for all possible kind of phase transitions in so large dimensional dynamical systems, without any idea of what can happen, seems to be desperate, even for model I and II, especially when thinking of what can already happen in one dimensional systems (Beck & Schloegl, 1995). So we shall not discuss within more details this aspect, focusing on grammar changes 7 .
Indeed, from the point of view of neural network analysis, grammar changes implies changes in the type of raster plots that the network is able to produce. An interesting situation occurs when the set of admissible raster plots obtained after adaptation belongs to Σ γ (τ ) ∩ Σ γ (τ +1) . In this case, adaptation 6 Additional constraints are required ensuring that F (τ ) φ does not tend to −∞. Typically, such constraints amount to bound the synaptic weights variation (soft or hard-bounds rules) by a suitable choice of functions h ijL 1 L 2 in eq. (13). 7 We actually conjecture that the only possible phase transitions in model I and II are grammar changes occurring when crossing the set where d(Ω, S) = 0 (eq. (6)).
plays the role of a selective mechanism where the set of admissible raster plots, viewed as a neural code, is gradually reducing, producing after n steps of adaptation a set ∩ n m=1 Σ γ (m) which can be rather small. This has been observed by Soula et al. in (Soula et al., 2006) for model I with Spike Time Dependent Plasticity, though not analyzed in those terms. If we consider the situation where (1) is a neural network submitted to some stimulus, where a raster plot ω encodes the spike response to the stimulus, then Σ γ is the set of all possible raster plots encoding this stimulus. Adaptation results in a reduction of the possible coding, thus reducing the variability in the possible responses.

Evolution of the potential.
Iterating the adaptation process, one expects to have periods of smooth variations, punctuated by sharp transitions where potential and grammar change. We qualify them under the generic name of "phase transitions" without further specifications. We write τ ≡ (l, n) where l indexes the phase transition and n the number of epoch since the last phase transition. The adaptation process corresponds now to a sequence ψ (τ ) ≡ ψ (l,n) of Gibbs potentials. Thus, ψ (l,0) is the potential arising just after the l-th transition, with the convention that ψ (l,0) (ω) = −∞ if ω is forbidden, such that ψ (l,0) characterizes the grammar after the l-th phase transition. We call a regular period the succession of adaptation epochs between two phase transitions.
By definition, during a regular period the synaptic update W (τ +1) = W (τ ) + δW (τ ) induces a smooth variation of the potential ψ (τ +1) = ψ (τ ) + δψ (τ ) , and the variational principle (43) selects a new Gibbs measure ν ψ (τ +1) with a potential: It follows that: where τ l is the epoch where the l-th phase transition aroused. The explicit form of δψ (τ ) (as well as ψ (τ ) ) depends on the detailed form of the dynamics, and there is little hope to determine it in general, except when the adaptation process converges to a solution with δW = 0 (see section 3.5.1). According to section 3.6.2, the function F (τ ) φ (eq. (24)) decreases during regular periods.
When a phase transition occurs, the relation (29) does not hold anymore. We now write: where n l is the last epoch before the phase transition l + 1. δψ (l+1,0) reg contains the regular variations of the potential while δψ (l+1,0) sing contains the singular variations corresponding to a change of grammar. We end up with the following picture. During regular periods the grammar does not evolve, the potential changes according to (30), (31), and the function F (τ ) φ decreases. Then, there is a sharp change in the grammar and the potential. The function F (τ ) φ may also have singularities and may sharply increase. If the adaptation rule converges, the potential (30) converges to ψ * = Φ + λ * .φ (eq. (21)). This potential Φ characterizes the grammar and the statistical weight of allowed raster plots is given by λ * .φ. The potential ψ * contains all changes in the grammar and statistics arising during the evolution (30,31). Thus it contains the history of the system and the evolution of the grammar. As an example we consider an adaptation rule inspired from (15) with an additional term r d W where f (x) is given by (16) and with: Set : , and: where φ ij is a finite range potential with range 2T s . Then (32) has the form (12), δW

Effects of the adaptation rule on the synaptic weights distribution.
The term S i (ω) (eq. 33) can be either negative, inducing Long Term Depression, or positive inducing Long Term Potentiation. In particular, its average with respect to the empirical measure π (T ) ω (τ ) reads: where: and where r i (τ ) = π (T ) ω (τ ) (ω i ) is the frequency rate of neuron i in the τ -th adaptation epoch. The term η neither depend on ω nor on τ , but only on the adaptation rule parameters A − , A + , τ − , τ + , T s . Equation (34) makes explicit 3 regimes.
• Cooperative regime. If η > 0 then π (T ) ω (τ ) (S i ) > 0. Then synaptic weights have a tendency to become more positive. This corresponds to a cooperative system (Hirsch, 1989). When iterating adaptation, dynamics become trivial with neurons firing at each time step or remaining quiescent forever.
• Competitive regime. On the opposite if η < 0 synaptic weights become negative. This corresponds to a competitive system (Hirsch, 1989).
• Intermediate regime. The intermediate regime corresponds to η ∼ 0. Here no clear cut tendency can be distinguished from the average value of S i and spikes correlations have to be considered as well.

Static weights.
Thanks to the soft bound term r d W ij the synaptic adaptation rule admits a static solution given by: Note that this equation can have several solutions. Using the same decomposition as the one leading to (34) we obtain: Note that ω j (0) ω i (u) ≥ 0, thus the first term is negative and the second one is positive (see eq. (16)). The sign of W ij depend on the parameters A − , A + , T s , but also on the relative strength of the terms π (T )

Convergence to the static weights solution.
The synaptic adaptation rule (32) defines a mapping Q on the set of synaptic weights: where δ ij,kl = 1 if i = k, j = l and 0 otherwise. The second term is a linear response characterizing the variation of π with respect to small variations of W kl . Approximating π (T ) ω (τ ) by a Gibbs distribution with a potential ψ (τ ) , it writes C ψ (τ ) ,Hij (0) + 2 +∞ t=0 C ψ (τ ) ,Hij (t) where C ψ (τ ) ,Hij (t) is the time-t correlation function between ψ (τ ) and H ij . When ν ψ (τ ) is smooth with respect to W, the derivative is dominated by the first term and Q is contracting for sufficiently small ǫ. This is however not a sufficient condition for the convergence of (32) because Q is not continuous everywhere. It has jumps whenever γ crosses the boundary of a structural stability domain.
If the static solution is contained in a open ball where Q is continuous the contraction property ensures the convergence to a static solution for any initial condition in this ball (Brouwer theorem).
In this paper we postulate the convergence and check it numerically.
4.1.5 Spike train statistics in a static weights regime.
As emphasized in section 2.6.5 and 3.5.1, when the synaptic adaptation rule converges to a fixed point, the corresponding statistical model is a Gibbs measure with a potential: where Φ contains the grammar and λ are free statistical parameters. The value λ * of these parameters in the potential ψ * is determined by the relation: where the pressure is given by: This procedure provides us the explicit form of the raster plot probability distribution when the adaptation rule converges. But the price to pay is that we have to determine simultaneously the N 2 parameters λ * ij on which the Gibbs measure depends. Focusing on the joint probability of a small set of neurons (pairs, triplets) this constraint can be relaxed in order to be numerically tractable.

Numerical checks.
The main goal of section 4 is to provide an example of the theoretical concepts developed in this paper. The emphasis is however not put on numerical results which will be the main topic of a forthcoming paper. Therefore, we focus here on numerical simulations for model I only, since the simulations for model II are quite more computer-time consuming, and we consider only one case of η value corresponding to the intermediate regime defined above.

Implementation.
Neurons dynamics. Previous numerical explorations have shown that a Model I-network of N neurons, with synapses taken randomly from a distribution N (0, C 2 N ), where C is a control parameter, exhibits a dynamics with very large periods in specific regions of values of the space (ρ, C) (Cessac, 2008;Cessac & Viéville, 2008). On this basis, we choose N = 100, ρ = 0.95, C = 4.0, N = 100. The external current I (ext) in eq. (5) is given by I ext i = 0.06 + 0.01 N (0, 1). Note that fixing a sufficiently large average value for this current avoids a situation where neurons stops firing after a certain time ("neural death").
STDP implementation. An efficient implementation of a STDP rule does not only depend on the analytic form of the rule but also on the respective time scales characterizing neurons and synapses evolutions. In "on-line-protocols", where the synaptic changes occur at the same time scale as neurons dynamics the so-called recursive approach appears to be the more appropriated (see (Zou, 2006) and references therein). In our case, we use instead an offline protocol, where we register the dynamics of the system on a long time windows, then compute the STDP modification, then let the system evolve again to its attractor (see section 3). The offline protocols are very expensive in machine-time, especially when using long spike trains to get a reliable time average (π (T ) ω (τ ) (φ ij )). For this reason we need to add some words about our implementation.
We register spike trains in a binary code. Indeed, this is the cheapest way in memory requirements, though it might be expensive for accessing specific spike times. Also, bit-wise operations are faster than their equivalents on other types of data. Finally, there exist very fast methods for computing the number of bits on any specific variable of type of integer (The faster for large number of iterations is a look-up table of precomputed values, but in-line methods -using parallel methods based on masksare not so far in terms of performances. For details and speed comparison see the G. Manku website http://infolab.stanford.edu/ manku/bitcount/bitcount.html). Numerical comparison of this method with the direct one that records the dynamics in Boolean arrays and compute STDP spike by spike shows enormous performance difference growing exponentially as the length of trains increases.
Fix a potential ψ. If T is the length of the experimental raster plot, one divides it into k blocs of length n, such that T = k × n. Call: with S (n) j (ω) = jn+n−1 t=jn ψ(σ t ω). Then, the following result can be proved (Chazottes, 1999): for ν ψ almost-every ω. This computation requires the knowledge of ψ.
The result in (Chazottes, 1999) can be extended straightforwardly to the following case. If ψ 1 , ψ 2 are two potentials for the same grammar, and δψ = ψ 2 − ψ 1 then: Expanding this equation in power series of ǫ one recovers the expansion (26). In particular, subtracting to (41) the leading term (of order δW (τ ) .δW (τ ) ) one gets an estimation of the "linear response" term κ (τ ) i1,j1,i2,j2 in eq. (27). However, on practical grounds, the computation of (41) requires very long time series to reduce significantly the finite fluctuations of the empirical average. Similar problems are encountered in the computation of linear response from other methods (Cessac, 2007).
This result holds for any adaptation rule (any polynomial φ). In the case of the adaptation rule (32), this reduces to computing:

Simulation results.
We have run simulations for STDP parameters values: r d = 0.99, ǫ = 0.001, τ + = 16, τ − = 32, A + = 1.0, and T s = 64, corresponding to standard values (Izhikevich & Desai, 2003). We fix η = −0.01 corresponding to the intermediate regime discussed in section 4.1.2. This fixes the value of A − via eq. (35). The length of spike trains is T = 3941 = 4096 − 128 = 128 × 32 − 2 × T s , where 32 is the number of bits in long integer, in our (machine-dependent) implementation. An extended description will be published elsewhere. The main results are summarized in fig. 3. In this figure (Top) we have represented the evolution of the distribution of synaptic weights W (τ ) (see captions). For this value of η, the histogram evolves to a unimodal distribution with some nodes having strong synaptic weights. In fig. 3 (center) we depict the raster plots at the beginning and at the end of synaptic adaptation process. In fig. 3 (bottom-left) we have shown the evolution of the Frobenius norm for δW (τ ) (i.e, N i,j |δW (τ ) ij | 2 ), which converges to 0. This shows the convergence of the rule in this case. We have also plotted the average weights modification 1 N 2 i,j δW ij . Finally in fig. 3 (bottom -right) is represented δ (τ ) F φ computed with 512 block of size 512. We see clearly that this quantity is negative as expected from eq. (28) and tends to 0. We have also represented the term − 1 ǫ δW (τ ) .δW (τ ) corresponding to the first order expansion in eq. (28).

Conclusion.
In this paper, we have considered the questions of characterizing the statistics of spike trains generated by neuronal networks, and how these statistics evolve during synaptic adaptation, from the angle of dynamical systems theory and ergodic theory. We have introduced a framework where notions such as raster plots, empirical averages of observables, and synaptic plasticity rules, coming from neuroscience, can be formalized in the language of dynamical systems. In this context, we propose a strategy to produce optimal statistical models of spike trains, based on the maximization of statistical entropy where expectation of observables must be compatible with their empirical averages. These models are Gibbs measures. Considering adaptation rules with slow dynamics, we also suggest that such Gibbs measures may result from related adaptation rules whenever synaptic weights converge to a fix value. In this case, synaptic weights resume the whole history of the neuronal network contained in the structure of the generating function of cumulants (the topological pressure) and in the structure of allowed/forbidden raster plots (grammar).
ij (red) and the Frobenius norm for STDP matrix on log scale i.e, log 10 "

29
To our opinion this theoretical paper is a beginning. Our hope is now to apply this strategy to data coming from biological neuronal networks, while only formal neural networks were considered here. At the present stage there is indeed a crucial issue to be able to characterize spike train statistics and to discriminate statistical models. For example, it is still an open problem to decide, for a given spike train in a given experiment, whether firing rate are sufficient for the statistical characterization (leading to uncorrelated Poisson models) or if higher order orders statistics such as correlations are also necessary. This is more than an academic question. Beyond the choice of such or such statistical model, there are hypotheses on how neuronal network convert stimuli from the external word into an action. Discriminating statistical models is the way to select the best hypothesis.
Our approach opens up the possibility of producing numerical methods for such a discrimination, e.g. based on Kullback-Leibler divergence minimization obtained from eq. (9). Some of these methods are already available as open source code on http://enas.gforge.inria.fr/. They are based on the theory developed in the present paper, relying itself on ergodic theory and thermodynamic formalism. At the present stage, the overwhelming richness of thermodynamic formalism has not really been considered in neuroscience. We hope that this paper will attract the attention of the neuroscience community on this point.

Appendix
In this appendix we present the thermodynamic formalism material used in the paper. Though all quantities defined below depend on γ, we have dropped this dependence, which is not central here, in order to alleviate notations. We follow (Keller, 1998;Parry & Pollicott, 1990;Chazottes et al., 1998;Chazottes, 1999;Chazottes & Keller, 2009) in the presentation of Gibbs measures in ergodic theory. In accordance with the body of the text Σ is the set of admissible raster plots and we denote by B the related Borel sigma-field.

Ergodic measures.
Fix φ a continuous function. The time average of φ along the orbit of X is given by If µ is an invariant measure 8 , then, according to Birkhoff theorem, this limit exists for µ-almost every X. Standard theorems in ergodic theory ensure that (1) has at least one invariant measure (Katok & Hasselblatt, 1998) but there are typically many. Among them, ergodic measures play a distinguished role. An invariant measure µ is ergodic if any real measurable invariant function is µ-almost surely constant. As a corollary of Birkhoff's theorem π X (φ) = µ(φ), for µ-almost every X, where µ(φ) = φdµ is the expectation of φ with respect to µ. Namely, the empirical average π X (φ) is equal to the "ensemble" average µ(φ) whatever the initial condition, provided it is chosen in the support of µ. Any invariant measure can be written as a convex decomposition of ergodic measures.

Raster plots statistics.
Fix n > 0, a set of times t 1 , . . . t n , and a set of prescribed spiking patterns ω(t 1 ) . . . ω(t n ). The set of raster plots ω ′ ∈ Σ such that ω ′ (t k ) = ω(t k ), k = 1 . . . n contains therefore raster plots where the spiking patterns at prescribed times t 1 . . . t n are imposed (cylinder set). By (countable) intersection and union of cylinder sets one can generate all possible events in Σ, such as "neuron i is firing at time t 1 ", or "neuron i 1 is firing at time t 1 , and neuron i 2 is firing at time t 2 , . . . and neuron i n is firing at time t n ", etc . . . .
The statistical properties of raster plots are inherited from the statistical properties of orbits of (1), via the correspondence ν[C] = µ [{X ⇀ ω, ω ∈ C}], where C is a cylinder set. On practical grounds, these statistical properties are obtained by empirical average. The probability of a cylinder set C is given by: where χ is the indicatrix function, while the time average of some function φ : Σ → IR is given by: 8 Namely µ(Fγ −1 A) = Fγ (A) for any measurable set A ⊂ B.

31
When ν is ergodic the following holds.
for ν-almost every ω, where π ω (.) is called the empirical measure for the raster plot ω.
In this paper, we have adopted a pragmatic approach based on the idea that statistical models obtained via the applications of principles in thermodynamic formalism (which are basically the same as those used in statistical physics) could constitute good prototypes for the analysis of real spike trains statistics. In this way, our goal is not to prove theorems but instead to use already proved theorems. However, theorems have hypotheses that we now explicit and whose plausible validity is discussed. Basically, to apply the machinery described below we need first to assume that raster plots constitute a symbolic coding for the orbits of (1), i.e. there is a one-to-one correspondence between membrane potential trajectories and raster plots, except for a negligible subset of initial conditions. We also assume that symbolic trajectories (raster plots) are generated by a finite grammar. This hypothesis is discussed in section 2.6.1, and it can actually be rigorously established for models I and II without noise. The case with noise is under current investigations. The corresponding transition matrix G defines a sub-shift of finite type or topological Markov chain. Actually, part of the results presented below are extensions of standard results on Markov chains. We also assumed that G is primitive i.e. ∃n 0 > 0 such that ∀i, j, ∀n > n 0 (G n ) ij > 0. This amounts to assuming that the dynamics on each attractor is topologically mixing 9 .
All these hypotheses hold if the dynamical system generating the raster plots is continuous and uniformly hyperbolic. Following the definition of (Parry & Pollicott, 1990), a map F : M → M is uniformly hyperbolic on a set Λ if: 1. There exists a splitting of the tangent space in unstable (E u ) and stable (E s ) subspaces such that there exists C, µ > 0 with DF t E s , DF −t E u ≤ Cµ t , t > 0.
3. The periodic orbits in Λ are dense (and Λ consists of more than a single closed orbit).

There exists an open set
Adding noise to model I, II renders them uniformly hyperbolic (see footnote 4 in the text) but they are not continuous due to the singularity of the firing threshold. Hence, additional work is necessary to prove that there is a finite Markov partition (Cessac & Fernandez, in preparation). Throughout this Appendix we assume that all these hypotheses hold.
A specific example of potential, used in this paper, corresponds to integrating the grammar into a potential Φ such that Φ(ω) = 0 is ω is allowed, and Φ(ω) = −∞ otherwise. In this case, however, the potential is not continuous anymore but only upper semi-continuous.

Equilibrium states
Let ν be a σ-invariant measure and Σ (n) the set of admissible cylinders of lenght n. Call: n is the entropy of ν. Let m (inv) be the set of invariant measures for σ, then the pressure of a potential ψ is defined by: This supremum is attained-not necessarily at a unique measure-where: ν ψ is called an equilibrium state, as it maximises some version of the (neg) free energy. Equivalently, a measure ν is an equilibrium state for a potential ψ ∈ C(Σ) if and only if, for all potential η ∈ C(Σ): This means that ν is a tangent functional for P at ψ. In particular ψ has a unique equilibrium state if and only if P is differentiable at ψ, i.e. lim t→0 6.3.4 Properties of the topological pressure.
We review some important properties of the topological pressure used in this paper. Basically, the topological pressure has the same properties of thermodynamic potentials like free energy. First, where Z n (ψ) is a partition function: and S (n) ψ(ω) = n−1 t=0 ψ(σ t ω).
Also, the topological pressure is a convex functional of the potentials i.e.

Gibbs measure.
Assume that the transition matrix of the shift σ is primitive and that ψ is regular. Then, there is a unique ergodic measure ν ψ , called a Gibbs measure, for which one can find some constants c 1 , c 2 with 0 < c 1 ≤ 1 ≤ c 2 such that for all n ≥ 1 and for all ω ∈ Σ: where S (n) ψ(ω) is given by (47). ν ψ is also the unique equilibrium state.
The Ruelle-Perron-Frobenius (RPF) operator for the potential ψ, denoted by L ψ , acts on functions g ∈ C(Σ) as L ψ g(ω) = ω ′ , ω=σω ′ e ψ(ω ′ ) g(ω ′ ). This operator has a unique maximal eigenvalue s = e P [ψ] associated to a right eigenfunction b ψ and a left eigenfunction ρ ψ (probability measure) such that L ψ b ψ (ω) = sb ψ (ω) and L ψ vdρ ψ = s vdρ ψ , for all v ∈ C(Σ). The remaining part of the spectrum is located in a disk in the complex plane, of radius strictly lower than s. Finally, for all v ∈ C(Σ) The Gibbs measure is ν ψ = b ψ ρ ψ . Given ψ and knowing its right eigenfunction b ψ and pressure P [ψ] one can construct a new potential: such that L Ψ 1 = 1, where 1 is the constant function 1(ω) = 1. Such a potential is called "normalised". Its pressure is zero. If ψ has a finite range-r then the RPF operator reduces to the transition matrix of a (r + 1)-step Markov chain. Thus, b ψ and ρ ψ can be easily determined. The Gibbs measure is none other than the invariant mesure of this chain. This can be used to generate a raster plot distributed according to a Gibbs distribution with a potential ψ. Moreover, for ν ψ -almost-every raster plot ω: This allows to obtain of ψ numerically (Chazottes et al., 1998).
There is a last important property. Let µ be an invariant measure and ν ψ a Gibbs measure with a potential ψ, both defined on the same set of sequences Σ. Let be the relative entropy (or Kullack-Leibler divergence) between µ and ν. Then, If µ = ν, d(µ, ν) = 0. The converse may not be true if the potential is not regular. If a raster ω is typical for the Gibbs measure ν ψ then one expects that π (T ) ω becomes closer to ν ψ than any other Gibbs measure (for another potential ψ ′ ) as T growths. This provides a criterion to compare two Gibbs measures (and to discrimate between several statistical models). Indeed, the following theorem holds (Chazottes et al., 1998).
This says that d π (T ) ω , ν ψ becomes, for sufficiently large T , smaller than the Kullback-Leibler divergence between π (T ) ω and any other Gibbs measure.