Bruno Cessac
email: bruno.cessac@inria.fr

Adrian G Palacios
email: adrian.palacios@uv.clvii

Spike Train Statistics from Empirical Facts to Theory: The Case of the Retina

This chapter focuses on methods from statistical physics and probability theory allowing the analysis of spike trains in neural networks. Taking as an example the retina we present recent works attempting to understand how retina ganglion cells encode the information transmitted to the visual cortex via the optical nerve, by analyzing their spike train statistics. We compare the maximal entropy models used in the literature of retina spike train analysis to rigorous results establishing the exact form of spike train statistics in conductance-based Integrate-and-Fire neural networks.

Introduction

Given a stimulus from the external world (e.g., visual scene, sound or smell) biological sensors at the periphery of the nervous system are able to transduce the physical manifestations of this stimulus (light emission, air pressure variations, chemical concentrations) into sequences of action potentials (spike trains), which propagate through the nervous system. Then, the brain is able to analyze those spike trains and infer crucial information on the nature of the stimulus. Critical -yet unsolved -questions in neuroscience are How is the physical signal encoded by the nervous system? How does the brain analyze the spike trains? What are the underlying computational coding principles? At the current stage of scientific knowledge, answering those questions is still a challenge for biology and computational neuroscience.

Among sensory systems the retina provides functionality such as detection of movement, orientation, temporal and spatial prediction, response to flash omissions and contrast, that were up to recently viewed as the exclusive duty of higher brain centers [START_REF] Gollisch | Eye smarter than scientists believed: neural computations in circuits of the retina[END_REF]. The retina is an accessible part of the brain [START_REF] Dowling | The retina: an approachable part of the brain[END_REF] and a prominent system to study the neurobiology and the underlying computational capacity of the neural coding. As a matter of fact, there is currently a wide research activity in understanding how the retina encodes visual information. However, basic questions are still open, such as: Are the ganglion cells (which send spikes from the eyes to the brain via the optical nerve), independent signal-encoders or are neural correlations important for coding a visual scene, and how to interpret them?

Outline

In this chapter, we present a state of the art about neural coding in the retina considered from the point of view of statistical physics and probability theory. As a consequence, this chapter contains both recent biological results and mathematical developments. The chapter is organized as follows. In Sect. 1.2 we introduce the current challenge of unraveling the neural code via spike trains statistics analysis. Such an analysis requires elaborated mathematical tools introduced in Sect. 1.3. We mainly focus on the so-called Gibbs distributions. This concept come from statistical physics but our presentation departs from the classical physics courses since it is based on transition probabilities of Markov process. This way, as we show, allows to handle non-stationary dynamics, and is adapted to statistical analysis, of data as well as neural networks models. As an illustration, we present, in Sect. 1.4, two "success stories" where spike train statistics analysis has allowed to make a step further in our understanding of information encoding by the retina. In the same section, we also present an example of a rigorous spike train analysis in a neural network and compare the spike trains probability distribution to the models currently used on the experimental side. The vertebrate retina is a tightly packed neural tissue, exhibiting a rich diversity of neurons. It is structured in three cells nuclei layers and two plexiform synaptic layers [START_REF] Masland | The fundamental plan of the retina[END_REF][START_REF] Wohrer | Virtual retina: a biological retina model and simulator, with contrast gain control[END_REF] (Fig. 1.1(a)). The outer nuclear layer (ONL) contains the rods and cones photoreceptors (P) somata; the inner nuclear layer (INL) contains bipolar (B), horizontal (H) and amacrine cells (A). Finally, the most internal nuclear layer is composed with ganglion cells (G) and displaced amacrine cells. The outer plexiform layer (OPL) corresponds to synaptic contacts between P, B and H cells. The inner plexiform layer (IPL) corresponds to synaptic contacts between B, A and G cells.

The retina is about 300 -500µm thick, depending on species, and has about 100 -130 millions of photoreceptors, 10-12 millions of bipolar, horizontal and amacrines cells and 0.4 to 1.6 millions of G cells.

Together with this high and compact number of cells there is a very large number of synapses present in dendrites and axons terminal, that has been roughly estimated to 1 billion of synapses [START_REF] Sarpeshkar | Ultra Low Power Bioelectronics: Fundamentals, Biomedical Applications, and Bio-Inspired Systems[END_REF]. The retina is also rich in terms of the variability of neurotransmitters, where rods, cones, and bipolar cells liberate glutamate, horizontal and amacrines cells can liberate gaba, glycine, serotonin, acetylcholine, dopamine among others. Together with the richness in chemical slow synapses circuits the retina has a variety of electrical ("gap-junctions"), fast synapses endowing the retina with specific functional circuits.

Single photons are converted by photoreceptors into a graded change in the resting potential, resulting in a neurotransmitter liberation (glutamate) into the synaptic region connecting photoreceptors with B and H cells. Those cells make synapses with G and A cells. Therefore, photons fluxes generate a chain of changes in the resting potential of B,H,A, and G cells with consequence the emission of action potentials ("spikes") by G cells. They are the principal neural encoders through the integration of neural signals.

The retina output, formed by spike train sequence, is carried by different types of G cells through the optical nerve to the brain higher visual structures: e.g.lateral geniculate nucleus (LGN) or visual cortex layers (Fig. 1.2).

Retina Circuits and Receptive Fields

As a result of its stratified, horizontal and vertical structure, and of the various type of synaptic connections (electrical fast synapses ∼ 0.1 ms for short distance; chemical slow synapses ∼ 10 ms for long distances) between the different type of neurons (P,H,B,A,G) a large number of "circuits" are present in LGN) in the thalamus and then for a second synapses to the main cortical visual area V1, from where many other projections target secondary cortical areas (V2, etc). Reproduced from [START_REF] Logothetis | Vision: A window on consciousness[END_REF].

the retina. The main connectivity structure of the retina is shown in Fig. 1.1(b). This circuitry results in the capacity of specific G cells to respond to specific stimuli in the visual field.

The receptive field (RF) of a sensory neuron is a region of space where the presence of a stimulus modifies the activity of that neuron. In the retina this change of activity is precisely the result of the transduction chain, from photoreceptor to G cells, converting photons into spike trains. As a consequence, one also defines the RF of a G cell as the input from all of the photoreceptors which synapse with it (via B, H, A) cells.

The RF of a cell can have different forms, depending on the network of neurons connected to this cell.

A prominent example, is the antagonist ON-OFF center-surround arrangement. First, photoreceptors make synapses with ON (excitatory) B cells and OFF (inhibitory) B cells according to their response to light. The physiological properties of G cells are determined at the center and the surround of their RF by the input of ON or OFF B cell. As consequence of dynamical and complex interaction (spatial and temporal) opposite functions for, e.g., color, contrast, intensity are likewise found for a single G cell, depending where the stimulus is present in their RF. It has been long believed that retina was mainly acting as an image transducer, absorbing photons and producing electrical signals or acting as a temporal and spatial linear filter. It was also believed that the retina doesn't perform any pre-processing of the image before sending spike trains to the brain. More recently, researchers pointed out that retina, in some species, is "smarter" than previously believed and is able to detect salient features or properties in a image such as approaching motion, motion detection and discrimination, texture and object motion, creating predictive or anticipatory coding thanks to "specialized" G cells (see, e.g., [START_REF] Gollisch | Eye smarter than scientists believed: neural computations in circuits of the retina[END_REF] for a review). The specificity of these population of cells for the detection of differential motion results largely from the circuit they belong to. An example is shown in Fig. 1.4 (detection of differential motion) where A cells play a prominent role.

Multi-Electrodes Array Acquisition

The pioneering work of Hubel and Wiesel based on anatomy and single cell recording on brain visual areas was very useful. However at that time, little was known about the properties of the reti- MEA devices constitute an excellent tool to track the physiological properties of G cells [START_REF] Petrusca | Identification and characterization of a y-like primate retinal ganglion cell type[END_REF] as well as their coding capacity [START_REF] Koch | Efficiency of information transmission by retinal ganglion cells[END_REF][START_REF] Koch | How much the eye tells the brain[END_REF]. Before the introduction of MEA devices, the neural coding properties of single G cells was study using intra or extra cellular electrodes, giving a limited sense of their collective role. In that respect, the work of Markus Meister et al.. [START_REF] Meister | Multi-neuronal signals from the retina: acquisition and analysis[END_REF] using MEA devices was pioneer. With simple stimulus, like checkerboard random white noise, and spike sorting algorithms these authors were able to determinate the number of spiking cells and their respective RF. They have shown that concerted G cells are critical, not only for retina development, but for the neural coding processing.

Encoding a Visual Scene

The distribution and fluctuations of visual signals in the environment can be treated as a statistical problem [START_REF] Olshausen | Natural image statistics and efficient coding[END_REF][START_REF] Simoncelli | Natural image statistics and neural representation[END_REF]. Natural scenes (a digital image or movie from a natural scenario) differ in their particular statistical structure and therefore the encoding capacity of a visual system should be able to match the properties and distribution of visual signals in the environment where the organism lives [START_REF] Atick | Could information theory provide an ecological theory of sensory processing?[END_REF][START_REF] Geisler | Visual perception and the statistical properties of natural scenes[END_REF][START_REF] Geisler | Natural systems analysis[END_REF][START_REF] Thompson | Ways of coloring: Comparative color vision as case study for cognitive science[END_REF][START_REF] Van Hemmen | 23 problems in systems neuroscience[END_REF].

The anatomical and physiological segregation of different aspects of a visual scene in separate spatial, temporal and chromatic channels start at the retina and rely on local "circuits" [START_REF] Atick | Could information theory provide an ecological theory of sensory processing?[END_REF]. However, how the precise articulation of this neural network contributes to local solutions and global perception is still largely a mystery.

G cells, as well as most neurons in the nervous system respond to excitations, coming from other neurons or from external stimuli, by emitting spike trains. In the contemporary view [START_REF] Rieke | Spikes: Exploring the Neural Code[END_REF], spikes are quantum or bits of information and their spatial (neuron-dependent) and temporal (spike times) structure carry "information": This is called "the" "neural code". Although this view is strongly based on a contemporary analogy with computers, spike trains are not computer-like binary codes. Indeed, an experiment reproduced several times (e.g., an image presented several times to the retina) does not reproduce exactly the same spike train, although some regularity is observed. As a consequence, current attempts to deciphering the neural code are based on statistical models.

The Ganglion Cells Diversity

The recent use of MEA in retina has lead to the description of a diversity of G cells type and to the question about their actual early visual capacity. The vertebrate retina has in fact 15-22 anatomically different class of G cells making it a much more complex functional neural network than expected [START_REF] Masland | The fundamental plan of the retina[END_REF][START_REF] Rockhill | The diversity of ganglion cells in a mammalian retina[END_REF][START_REF] Gollisch | Eye smarter than scientists believed: neural computations in circuits of the retina[END_REF].

The three most frequent G cells in the retina can be classified from their morphology in: parasol (primates but α or Y in cats and rabbits) corresponding to 3-8% of the total number of G cells; midget (β or X in cats and rabbits) corresponding to 45-70%; and bi-stratified G cells. In physiological terms parasol (Y) cells can be classified as brisk-transient, and midget (X) as brisk-sustained. They can have an ON or OFF function.

Although only a reduced fraction of the existing G cells [START_REF] Masland | The fundamental plan of the retina[END_REF][START_REF] Masland | The unsolved mystery of vision[END_REF][START_REF] Schwartz | Sophisticated temporal pattern recognition in retinal ganglion cells[END_REF] has been studied in detail [START_REF] Gollisch | Eye smarter than scientists believed: neural computations in circuits of the retina[END_REF],

their diversity raises questions such as: How do G cells encode an image? Which features from a natural visual scene are they coding? Are G cells independent or collective encoders?

An interesting approach has been advanced by [START_REF] Balasubramanian | Receptive fields and functional architecture in the retina[END_REF]. The authors propose that the retina organization should use simple coding principles to carry the maximum of information at low energetic cost. However, as the authors point out, the statistic distribution (e.g., color, contrast) for natural images is not Gaussian. Therefore, the classical Gaussian estimator for Shannon information:

I = 1 2 log 2 (1 + SN R) , 1 Neural networks dynamics. xiii
where SNR is the signal to noise ratio is not appropriate. Instead, "pixels" in natural images are highly correlated and the general form of statistical entropy (see Eq. ( 2) in [START_REF] Balasubramanian | Receptive fields and functional architecture in the retina[END_REF]) is required to calculate the spike capacity of G cells to carry information. In that respect, the coding capacity for different G cells has been estimated (see, e.g., Eq. ( 5) in [START_REF] Balasubramanian | Receptive fields and functional architecture in the retina[END_REF] ). The larger capacity for information transmission comes from, e.g., "sluggish" G cells (32%); local-edge (16%), brisk-transient (9%).

Population Code

This term refers to the computational capacity of a neural assembly (or circuit) to solve specific questions [START_REF] Averbeck | Neural correlations, population coding and computation[END_REF][START_REF] Pouget | Information processing with population codes[END_REF]. Assuming that living systems have evolved to optimize the population code, how is this optimum reached in the retina? Are G cells sensors independent encoders or, on the opposite, are neural correlations important for coding? In an influential article, Nirenberg et al.. [START_REF] Nirenberg | Retinal ganglion cells act largely as independent encoders[END_REF] suggest that G cells act as independent encoder. However, orchestrated spikes train from G cells were reported by pioneer work of Rodieck [START_REF] Rodieck | Maintained activity of cat retinal ganglion cells[END_REF] and Mastronarde [START_REF] Mastronarde | Correlated firing of cat retinal ganglion cells. I. Spontaneously active inputs to X-and Y-cells[END_REF]. Mastronarde shows that G cells responses tend to fire together and dynamically adapt to light or dark background [START_REF] Mastronarde | Correlated firing of cat retinal ganglion cells. I. Spontaneously active inputs to X-and Y-cells[END_REF]. This suggests that they act in a correlated way. However, this approach is by itself incomplete, since different sources of correlation were not clearly considered [START_REF] Panzeri | A unified approach to the study of temporal, correlational, and rate coding[END_REF][START_REF] Schneidman | Synergy, redundancy, and independence in population codes[END_REF].

On the other hand, MEA can now record many G cells from small pieces of retina (< 500µm) [START_REF] Devries | Correlated firing in rabbit retinal ganglion cells[END_REF][START_REF] Meister | Multi-neuronal signals from the retina: acquisition and analysis[END_REF][START_REF] Field | Information processing in the primate retina: circuitry and coding[END_REF] and help us to asses the importance, and origin, of neural synchrony for the neural coding. For example, in darkness, salamander G cells shows 3 types of synchrony depending on the time laps: (i) a common photoreceptor source through B cells (broadcorrelation : 40 -100ms) (ii) A cells and G cells connected through gap junctions (medium : 10 -50ms) (iii) gap junction between G cells (narrow :< 1ms) [START_REF] Brivanlou | Mechanisms of concerted firing among retinal ganglion cells[END_REF].

At present and although a large bunch of experimental facts enlighten our knowledge about the retina structure as well as its activity, basic questions on the way how a visual scene is encoded by spike trains remain still open. This is largely due to (i) the complex structure of the retina; (ii) its large number of cells; (iii) the lack of sufficiently accurate statistical models and methods to discriminate competing hypotheses. Apparently elementary questions such as determining whether correlations are significant from the analysis of MEA recordings requires in fact the use of smart statistical analysis techniques, based on "statistical models" defined by a set of a priori hypothesis as we see in the next section.

Spike Train Statistics from a Theoretical Perspective

In this section we develop the mathematical framework to analyze spike train statistics. The collective neuron dynamics, which is generally submitted to noise, produces spike trains with randomness though some statistical regularity can be observed. Spike trains statistics is assumed to be summarized by an hidden probability µ characterizing the probability of spiking patterns. One current goal in experimental analysis of spike trains is to approximate µ from data. We describe here several theoretical tools allowing to handle this question. Our presentation is based on the notion of transition probabilities. In this context we introduce Gibbs distributions, which is one of the main theoretical concept of this chapter.

Gibbs distributions are usually considered in the stationary case where they are obtained from the maximal entropy principle. Their definition via transition probabilities, adopted in this chapter, affords the consideration of Gibbs distribution in the more general context of non-stationary dynamics with possibly infinite memory.

Spike Statistics

Raster Plots

We consider a network of N neurons. We assume that there is a minimal time scale δ > 0 corresponding to the minimal resolution of the spike time, constrained by biophysics and by measurements methods (typically δ ∼ 1 ms) [START_REF] Cessac | On dynamics of integrate-and-fire neural networks with adaptive conductances[END_REF][START_REF] Cessac | Overview of facts and issues about neural coding by spikes[END_REF]. Without loss of generality (change of time units) we set δ = 1, so that spikes are recorded at integer times. One then associates to each neuron k and each integer time n a variable ω k (n) = 1 if neuron k fires at time n and ω k (n) = 0 otherwise. A spiking pattern is a vector

ω(n) def = [ ω k (n) ] N k=1
which tells us which neurons are firing at time n. We note A = { 0, 1 } N the set of spiking patterns. A spike block is a finite ordered list of spiking patterns, written:

ω n2 n1 = { ω(n) } {n1≤n≤n2} ,
where spike times have been prescribed between the times n 1 to n 2 (i.e., n 2n 1 + 1 time steps). The depth of the block is the number of time steps where time has been prescribed (in the example this is

n 2 -n 1 + 1)
. The set of such blocks is A n2-n1+1 . Thus, there are 2 N n possible blocks with N neurons and depth n. For example, N = 3 neurons and n = 2 time steps the possible blocks are:

  0 0 0 0 0 0   ;   0 1 0 0 0 0   ;   1 0 0 0 0 0   ;   1 1 0 0 0 0   ; . . .   1 1 1 1 1 0   ;   1 1 1 1 1 1   .
We call a raster plot a bi-infinite sequence

ω def = {ω(n)} +∞ n=-∞
, of spiking patterns. This notion corresponds to its biological counterpart (Sect. 1.2.2) with the obvious difference that experimental raster plots are finite. The consideration of infinite sequences is more convenient on the mathematical side but, at several places, we discuss the effects of having finite experimental rasters on spike statistics estimation. The set of raster plots is denoted X = A ) .

Transition Probabilities

The probability that a neuron emits a spike at some time n depends on the history of the neural network.

However, it is impossible to know explicitly its form in the general case since it depends on the past evolution of all variables determining the neural network state. A possible simplification is to consider that this probability depends only on the spikes emitted in the past by the network. In this way, we are seeking a family of transition probabilities of the form P ω(n) ω n-1 n-D , the probability that the firing pattern ω(n) occurs at time n, given a past spiking sequence ω n-1 n-D . Here D is the memory depth of the probability, i.e., how far in the past does the transition probability depend on the past spike sequence.

We use here the convention that P ω(n

) ω n-1 n-D = P [ ω(n) ] if D = 0 (memory-less case).
Transition probabilities depend on the neural network characteristics such as neurons conductances, synaptic responses or external currents. They give information on the dynamics that takes place in the observed neural networks. Especially, they have a causal structure where the probability of an event depends on the past. This reflects underlying biophysical mechanisms in the neural networks which are also causal. The explicit computation of transition probabilities can be done in some model-examples (Sect. 1.4.4). From them, one is able to characterize statistical properties of rasters generated by the network, as we now develop.

Markov Chains

Transition probabilities with a finite memory depth D define a "Markov chain", i.e., a random process where the probability to be in some state at time n (here a spiking pattern ω(n)) depends only upon a finite past (here on ω(n -1), . . . , ω(n -D)). Markov chains have the following property. Assume that we know the probability of occurrence of the block ω m+D-1 m ,

P ω m+D-1 m = P [ ω(m + D -1), ω(m + D -2), . . . , ω(m) ] . (1.1)
Note that, mathematically, the order of the spiking patterns does not matter in the right-hand side since we are dealing with a joint probability, but choosing this specific order is useful for subsequent explanations. Then, by definition, the probability of the block ω m+D m is:

P ω m+D m =P [ ω(m + D), ω(m + D -1), . . . , ω(m) ] =P [ ω(m + D) | ω(m + D -1), . . . , ω(m) ] P [ ω(m + D -1), . . . , ω(m) ] .
Thus:

P ω m+D m = P ω(m + D) ω m+D-1 m P ω m+D-1 m , 1 Neural networks dynamics.
xv and, by induction, the probability of a block ω n m , ∀n, nm ≥ D is given by:

P [ ω n m ] = n l=m+D P ω(l) ω l-1 l-D P ω m+D-1 m . (1.2)
Thus, knowing the probability of occurrence of the block ω m+D-1 m one can infer the probability of forthcoming blocks by the mere multiplication of transition probabilities.

Given the (joint) probability P [ ω n m ] the (marginal) probability of sub-blocks can be easily obtained,

since for m ≤ n 1 ≤ n 2 ≤ n, P ω n2 n1 = * (n1,n2) m,n P [ ω n m ] , (1.3)
where

* (n1,n2) m,n
means that we sum up over all possible spiking patterns in the interval { m, n } excluding the interval { n 1 , n 2 } (i.e., we sum up over all possible values of ω(n), . . . , ω(n 1 -1), ω(n 2 + 1), . . . , ω(n)).

As a consequence, from (1.2), (1.3), the probability of the spike block ω n n-D , of depth D, is:

P ω n n-D = * (n-D,n) m,n n l=m+D P ω(l) ω l-1 l-D P ω m+D-1 m . (1.4)
Knowing the probability of an initial block of depth D (here ω m+D-1 m

) one infers from this equation the probability of subsequent blocks of depth D. Equation (1.4) can also be expressed in terms of vectormatrices multiplication, and the main properties of the Markov chain can be deduced from linear algebra and matrices spectra theorems [START_REF] Seneta | Non-negative Matrices and Markov Chains[END_REF]. For compactness we shall not use this possibility here though, (see [START_REF] Vasquez | Entropy-based parametric estimation of spike train statistics[END_REF] for further details).

However, this equation shows that the "future" of the Markov chain (the probability of occurrence of blocks) depends on an initial condition (here

P ω m+D-1 m
), which is a priori undetermined. Moreover, there are a priori infinitely many possible choices for the initial probability.

Asymptotic of the Markov Chain

Assume now that nm → +∞ in Eq. (1.4) and more precisely that m → -∞. Practically, this limit corresponds to considering that the system began to exist in a distant past (defined by the initial condition of the Markov chain) and that it has evolved long enough, i.e., over a time larger than relaxation times in the system, so that it has reached sort of an "adult" age where its structure is essentially fixed. Note that this does not exclude adaptation processes, e.g., if the transition probabilities depend explicitly on time.

Mathematically, the limit nm → +∞ corresponds to studying the asymptotic of the Markov chain and related questions are: Is there a limit probability P ω n n-D ? Does it depend on the initial condition

P ω m+D-1 m , m → -∞?
Let us first consider the easiest case where transition probabilities are invariant under time translation. This means that for each possible spiking pattern α ∈ A , for all possible "memory" blocks

α -1 -D ∈ A D and ∀n, P ω(n) = α ω n-1 n-D = α -1 -D = P ω(0) = α ω -1 -D = α -1 -D .
We call this property stationarity referring rather to the physics literature than to the Markov chains literature (where this property is called homogeneity). If, additionally, all transition probabilities are strictly positive then there is a unique probability µ, called the asymptotic probability of the chain, such that, whatever the initial choice of a probability P ω m+D-1 m in (1.4) the probability of a block ω n n-D converges to µ ω n n-D as m tends to -∞. One says that the chain is ergodic (Note that positivity of all transition probabilities is a sufficient but not necessary condition for ergodicity [START_REF] Seneta | Non-negative Matrices and Markov Chains[END_REF]). In this sense, dynamics somewhat "selects" the probability µ, since, whatever the initial condition P ω m+D-1 m , it provides the statistics of spikes observed after a sufficiently long time. Additionally, µ has the following property: for any time

n 1 , n 2 , n 2 -n 1 ≥ D, µ ω n2 n1 = n2 l=n1+D P ω(l) ω l-1 l-D µ ω n1+D-1 n1 . (1.5)
Let us return to the problem of choosing the initial probability in Eq. (1.4). If one wants to determine the evolution of the Markov chain after a initial observation time n 1 one has to fix the initial probability

P ω n1+D-1 n1
and to use (1.4) (where m is replaced by n 1 ) and there is an indeterminacy in the choice

of P ω n1+D-1 n1
. This indeterminacy is released, though, if the system has started to exist in the infinite past. Then, P ω n1+D-1 n1 has to be replaced by µ ω n1+D-1 n1 and Eq. (1.4) becomes:

µ ω n2 n2-D = * (n2-D,n2) n1,n2 n2 l=n1+D P ω(l) ω l-1 l-D µ ω n1+D-1 n1 . (1.6)
In this way, taking the limit m → -∞ for an ergodic Markov chain, resolves the indeterminacy in the initial condition.

Positivity and stationary assumptions may not hold. If positivity is violated then several situations can arise: several asymptotic probability distributions,can exist, depending on the choice of the initial probability P ω m+D-1 m

; it can also be that no asymptotic probability exist at all. If stationarity does not hold, as it is the case e.g.for a neural network with a time-dependent stimulus, then one can nevertheless define a probability µ selected by dynamics. In short, this is a probability µ on the set of raster plots A )

which still obeys (1.5) but without the conditions of stationarity (transition probabilities are not timetranslation invariant [START_REF] Fernandez | Chains with complete connections : General theory, uniqueness, loss of memory and mixing properties[END_REF]). In this case, which is realistic when dealing with living systems submitted, e.g., to time-dependent stimuli, the statistics of spikes is time-dependent. For example, the probability that a neuron emits a spike at time n depends on n, while it is not the case when dynamics is stationary.

1.3.1.5 Gibbs Distributions Assume that P ω(n) ω n-1 n-D
> 0 for all n ∈ ). Then, a probability distribution µ that obeys (1.5) is called a Gibbs distribution, and the function has to be modified [START_REF] Fernandez | Chains with complete connections : General theory, uniqueness, loss of memory and mixing properties[END_REF].

φ n ω n n-D def = log P ω(n) ω n-1 n-D , (1.7
The Gibbs potential depends on the block ω n-1 n-D and on the spiking pattern ω(n), thus, finally, this is a function of the block ω n n-D of depth D + 1. The term "normalized" refers to the fact that the potential in (1.7) is the logarithm of a transition probability. Below, we give example of Gibbs distributions where the potential is not normalized: this is an arbitrary function of the block ω n n-D . We call R = D + 1 the range of the potential. A Gibbs potential can have an infinite range (D → -∞ in our setting).

The condition P ω(n) ω n-1 n-D
> 0 for all n ∈ ) ensures that there is a one-to-one correspondence between a Gibbs potential and a Gibbs distribution. If this condition is relaxed, i.e., some transitions are forbidden, then several Gibbs distribution can be associated with a Gibbs potential. This corresponds to a first-order phase transition in statistical physics [START_REF] Georgii | Gibbs measures and phase transitions[END_REF]. In the infinite range case, the existence and uniqueness of a Gibbs distribution associated with this potential requires additional assumptions to the positivity of transition probabilities [START_REF] Fernandez | Chains with complete connections : General theory, uniqueness, loss of memory and mixing properties[END_REF].

From (1.2), we have ∀nm ≥ D:

µ ω n m | ω m+D-1 m = exp n l=m+D φ l ( l, ω ) . (1.8)
This form reminds the Gibbs distribution on spin lattices in statistical physics where one looks for lattice translation-invariant probability distributions given specific boundary conditions. Given a potential of range D the probability of a spin block depends on the states of spins in a neighborhood of size D of that block. Thus, the conditional probability of this block given a fixed neighborhood is the exponential of the energy characterizing physical interactions within the block as well as with the boundaries. Here, spins are replaced by spiking patterns; space is replaced with time which is mono-dimensional and oriented:

there is no dependence in the future. Boundary conditions are replaced by the dependence in the past.

(ω) = ω k (n) ω k ′ (n) is '1' if

Probabilities and Averages

Let µ be a probability on the set of rasters (typically the Gibbs distribution introduced above). Mathematically, the knowledge of µ is equivalent to knowing the probability µ [ ω n m ] of any possible spike block.

For an observable f we denote µ [ f ] def = f dµ the average of f with respect to µ. If f is only a function of finite blocks ω n m then:

µ [ f ] = ω n m f (ω n m ) µ [ ω n m ] , (1.9)
where the sum holds on all possible (2 n-m+1 ) values of ω n m . For example the average value of

f (ω) = ω k (n) is given by µ [ ω k (n) ] = ω k (n) ω k (n) µ [ ω k (n) ]
where the sum holds on all possible values of ω k (n) (0 or 1). Thus, finally

µ [ ω k (n) ] = µ [ ω k (n) = 1 ] , (1.10)
which is the probability of firing of neuron k at time n. This quantity is called the instantaneous firing rate. Likewise, the average value of ω k1 (n) ω k2 (n) is the probability that neuron k 1 and k 2 fire at the same time n: this is a measure of pairwise synchronization at time n.

Empirical Averages

In experiments, raster plots have a finite duration T and one has only access to a finite number N of those rasters, denoted ω (1) , . . . , ω (N ) . From these data one computes empirical averages of observables.

Depending on the hypotheses made on the underlying system there are several ways of computing those averages.

A classical (though questionable assumption as far as experiments are concerned) is stationarity: the statistics of spike is time-translation invariant. In this case the empirical average reduces to a time average.

We denote π (T )

ω [ f ] the time average of the function f computed for the raster ω of T . For example, when

f (ω) = ω k (n), π (T ) ω [ f ] = 1 T T -1 n=0 ω k (n)
, which provides an estimation of the firing rate of neuron k (it is independent of time from the stationarity assumption). If f is a monomial ω k1 (n 1 ) . . . ω km (n m ),

1 ≤ n 1 ≤ n 2 ≤ n m < T , then π (T ) ω [ f ] = 1 T -nm T -nm n=0 ω k1 (n 1 + n) . . . ω km (n m + n),
ω (1) [ f ] = π (T ) ω (2) [ f ] for distinct rasters
ω (1) , ω (2) . Moreover, those fluctuations depend on T .

Assume now that all empirical rasters have all been generated by an hidden Markov chain and additionally that this chain is ergodic with a Gibbs distribution µ. Then, all those rasters obey π (T )

ω (r) [ f ] → µ [ f ],
r = 1, . . . , N , as T → +∞, whatever f : the time average converges to the average with respect to the hidden probability µ (this is one of the definitions of ergodicity). As a consequence the fluctuations of the time-average about the exact mean µ [ f ] tends to 0, typically like

K f √
T , where K f is a constant depending on f . This is the celebrated central limit theorem stating moreover that fluctuations about the mean are Gaussian [START_REF] Gikhman | The Theory of Stochastic Processes[END_REF]. We come back to this point in Sect. 1.3.2.4.

The remarkable consequence of ergodicity (which implies stationarity) is that the empirical average can be estimated from one raster only. Now, if we have N rasters available we can use them to enlarge artificially the sample size, e.g.computing empirical average by 1

N N r=1 π (T ) ω (r) [ f ].
This also allows the computation of error bars as well as more elaborated statistical estimation techniques [START_REF] Rieke | Spikes: Exploring the Neural Code[END_REF].

What if the stationarity assumption is violated? Then, the average of f depends on time and one computes the empirical average from the N rasters. We denote π (N ) [ f (n) ] the average of f at time n, performed over N rasters. For example when

f (ω) = ω k (n), π (N ) [ f (n) ] = 1 N N r=1 ω (r)
k (n) is the sample-averaged probability that neuron k fires at time n . If all rasters are described by the same probability (the Gibbs distribution which is also defined in the non-stationary case), then

π (N ) [ f (n) ] → µ [ f (n) ] as N → +∞.

Example of Empirical Average: Estimating Instantaneous Pairwise Correlations

Assume that spikes are distributed according to an hidden probability µ supposed to be stationary for simplicity. The instantaneous pairwise correlations of neurons k, j with respect to µ is:

C ( k, j ) = µ [ ω k (0) ω j (0) ] -µ [ ω k (0) ] µ [ ω j (0) ] . (1.11)
Since µ is stationary the index 0 can be replaced by any time index (time-translation invariance of statistics).

Assume now that we have a raster ω distributed according to µ. An estimator of C ( k, j ) is:

C (T ) ω ( k, j ) = π (T ) ω [ ω k (0)ω j (0) ] -π (T ) ω [ ω k (0) ] π (T ) ω [ ω j (0) ] . (1.12)
It converges to C ( k, j ) as T → +∞.

The events "neuron k fires at time 0" (ω k (n) = 1) and "neuron j fires at time 0" (ω j (n) = 1) are

independent if µ [ ω k (0) ω j (0) ] = µ [ ω k (0) ] µ [ ω j (0) ],
thus C ( k, j ) = 0. (Note that independence implies vanishing correlation but the reverse is not true in general. Here the two properties are equivalent thanks to the binary 0, 1 form of the random variables ω k (0), ω j (0)).

Assume now that the observed raster has been drawn from a probability where these events are independent, but the experimentalist who analyzes this raster does not know it. To check independence she computes C (T ) ω ( k, j ) from the experimental raster ω. However, since T is finite, C

ω ( k, j ) will not be exactly 0. More precisely, from the central limit theorem the following holds. The probability that the random variable

C (T ) ω ( k, j
) is larger than ǫ, is well approximated (for large T and small ǫ) by e -ǫ 2 T 2 K .

K can be exactly computed (Sect. 1.3.1.5). In the simplest case where spikes are drawn independently with a probability p of having a spike, K is equal to p 2 (1p 2 ). Thus, fluctuations are Gaussian and their mean-square deviation decay with T as K T . As a consequence, even if neuron j and k are independent, the quantity C (T ) ω ( k, j ) will never be 0: it has fluctuations around 0.

This can be seen by a short computer program drawing at random 0's and 1's independently, with the probability p to have a '1', and plotting C As a consequence, it is stricto-sensu not possible to determine whether random variables are uncorrelated, by only computing the empirical correlation from samples of size T , since even if these variables are uncorrelated, the empirical correlation will never be zero. There exist statistical tests of independence from empirical data, beyond the scope of this chapter. A simple test consists of plotting the empirical correlation versus T and check whether it tends to zero as K T . Now, experiments affords only sample of limited size, where T rarely exceeds 10 6 . So, fluctuations are of order √ K ×10 -3 and it makes a difference whether K is small or big. It is therefore difficult to interpret weak empirical correlations. Are they sample fluctuations of a system where neurons are indeed independent, or are they really significant, although weak? This issue is further addressed in Sect. 1.4.2.

Matching Experimental Averages

Assume that an experimentalist observes N rasters, and assume that all those rasters are distributed according to an hidden probability distribution µ. Is to possible to determine or, at least, to approach µ from those rasters? One possibility relies on the maximal entropy principle described in the next sections.

We assume for the moment that statistics is stationary.

ω [ O k ] = C k , k = 1, .
. . , K. Now, to estimate the hidden probability µ, by some approximated probability µ ap , she has to assume, as a minimal requirement, that:

π (T ) ω [ O k ] = C k = µ ap [ O k ] , k = 1, . . . , K, (1.13)
i.e., the expected average of each observable, computed with respect to µ ap is equal to the average found in the experiment. This fixes a set of constraints to approach µ. We call µ ap a statistical model.

Unfortunately, this set of conditions does not fix a unique solution but infinitely many ! As an example if we have only one neuron whose firing rate is 12 , then a straightforward choice for µ ap is the probability where successive spikes are independent (P [

ω k (n)ω k (n -1) ] = P [ ω k (n) ] P [ ω k (n -1) ])
and where the probability of a spike is 1 2 . However, one can also take a one-step memory model where transition probabilities obey

P [ ω k (n) = 0 | ω k (n -1) = 0 ] = P [ ω k (n) = 1 | ω k (n -1) = 1 ] = p, P [ ω k (n) = 0 | ω k (n -1) = 1 ] = P [ ω k (n) = 1 | ω k (n -1) = 0 ] = 1 -p, p ∈ [0, 1]. In this case, indeed the invariant probability of the corresponding Markov chain is µ ap [ ω k (n) = 0, 1 ] = 1
2 , since from Eq. (1.5),

µ ap [ ω k (n) = 0 ] = ω k (n-1)=0,1 P [ ω k (n) = 0 | ω k (n -1) ] µ ap [ ω k (n -1) ] p 2 + 1 -p 2 = 1 2 .
The same holds for µ ap [ ω k (n) = 1 ]. In this case, we match the constraint too but with a model where successive spikes are not independent. Now, since p takes values in the interval [0, 1] there are uncountably many Markov chains with memory depth 1 matching the constraint. One could also likewise consider memory depth D = 2, 3 and so on.

Since transition probabilities reflect the underlying (causal) mechanisms taking place in the observed of the neural network, the choice of the statistical model defined by those transition probabilities is not anecdotal. In the example above, that can be easily generalized, one model considers that spikes are emitted like a coin tossing, without memory, while other models involve a causal structure with a memory of the past. Even worse, there are infinitely many choices for µ ap since (i) the memory depth can be arbitrary; (ii) for a given memory depth there are (infinitely) many Markov chains whose Gibbs distribution matches the constraints (1.13). Is there a way to selecting, in fine, only one model from constraints (1.13), by adding some additional requirement? The answer is "yes".

Entropy

The entropy rate or Kolmogorov-Sinai entropy of a stationary probability distribution µ is:

h [ µ ] = -lim n→∞ 1 n ω n 1 µ [ ω n 1 ] log µ [ ω n 1 ] , (1.14)
where the sum holds over all possible blocks ω n 1 . This definition holds for systems with finite or infinite memory. In the case of a Markov chain with memory depth D > 0, we have [START_REF] Cornfeld | Ergodic Theory[END_REF] h

[ µ ] = - ω D+1 1 µ ω D 1 P ω(D + 1) ω D 1 log P ω(D + 1) ω D 1 , (1.15)
Note that, from time-translation invariance the block ω D 1 can be replaced by ω D+n-1 n , for any integer n.

When D = 0, the entropy reduces to the classical definition:

h [ µ ] = - ω(0) µ [ ω(0) ] log µ [ ω(0) ] .
(1.16)

Gibbs Distributions in the Stationary Case

In the stationary case Gibbs distributions obey the following variational principle [START_REF] Ruelle | Thermodynamic formalism[END_REF][START_REF] Keller | Equilibrium States in Ergodic Theory[END_REF][START_REF] Chazottes | Pressure and Equilibrium States in Ergodic Theory, chapter Ergodic Theory[END_REF]:

P(φ) = sup ν∈Minv ( h [ ν ] + ν [ φ) ] ) = h [ µ ] + µ [ φ ] , (1.17)
where M inv is the set of all possible stationary probabilities ν on the set of rasters with

N neurons; h [ ν ]
is the entropy of ν and ν [ φ ] is the average value of φ with respect to the probability ν. Looking at the second equality, the variational principle (1.17) selects, among all possible probability ν, one probability which realizes the supremum, the Gibbs distribution µ.

The quantity P(φ) is called the topological pressure. For a normalized potential it is equal to 0. However, the variational principle (1.17) holds for non-normalized potentials as well i.e., functions which are not the logarithm of a probability [START_REF] Ruelle | Thermodynamic formalism[END_REF][START_REF] Keller | Equilibrium States in Ergodic Theory[END_REF][START_REF] Chazottes | Pressure and Equilibrium States in Ergodic Theory, chapter Ergodic Theory[END_REF].

In particular, consider a function of the form:

1 Neural networks dynamics.

xxi

H β (ω 0 -D ) = K k=1 β k O k (ω), (1.18)
where O k are observables, β k real numbers and β denotes the vector of β k 's, k = 1, . . . , K. We assume that each observable depends on spikes in a time interval { -D, . . . , 0 }.

To the non-normalized potential H β (ω 0 -D ) one can associate a normalized potential φ of the form:

φ(ω 0 -D ) = H β (ω 0 -D ) -log ζ β (ω 0 -D ), (1.19)
where ζ(ω 0 -D ) is a function that can explicitly computed. In short, one can associate to the potential

H β (ω 0 -D ) a matrix with positive coefficient; ζ β (ω 0 -D
) is is a function of the (real positive) largest eigenvalue of this matrix as well as of the corresponding right eigenvector (see [START_REF] Vasquez | Gibbs distribution analysis of temporal correlation structure on multicell spike trains from retina ganglion cells[END_REF] for details). This function depends on the model-parameters β. The topological pressure is the logarithm of the largest eigenvalue.

In this way, H β defines a stationary Markov chain with memory depth D, with transition probabilities:

P ω(0) ω -1 -D = e H β (ω 0 -D ) ζ β (ω -1 -D )
.

(1.20)

Denote µ β the Gibbs distribution of this Markov chain. The topological pressure P(φ β ) obeys:

∂P(φ β ) ∂β k = µ β [ O k ] , (1.21)
while its second derivative controls the covariance of the Gaussian matrix characterizing the fluctuations of empirical averages of observables about their mean. Note that those fluctuations are Gaussian if the second derivative of P is defined. This holds if all transitions probabilities are positive.

In the memory-less case D = 0 where only the statistics of instantaneous spiking patterns is considered, the Gibbs distribution reads:

µ β (ω(0)) = e H β (β,ω(0))
ω(0) e H β (β,ω(0)) .

(1.22)

In this case,

ζ β = ω(0)
e H β (β,ω(0)) .

(1.23) This is a constant (it does not depend on the raster). It is called partition function in statistical physics.

The Maximal Entropy Principle

Assume now that we want to approximate the exact (unknown) probability µ by an approximated probability µ ap that matches the constraints (1.13). The idea is to take as a statistical model µ ap the Gibbs distribution of a function of the form (1.18), corresponding to a set of constraints attached to observables O k , where the β k 's are free parameters of the model. Thus, the statistical model is fixed by the set of observables and by the value of β. We write then, from now on, µ β instead of µ ap .

Looking at the variational principle (1.17), we have to take the supremum over all probabilities ν that matches (1.13), i.e., µ

β [ O k ] = C k so that µ β [ H β ]
is a constant for fixed β. Therefore, in this case (1.13) reduces to maximizing the entropy rate given the constraints (1.13). This the classical way of introducing Gibbs distributions in physics courses. Then, the β k 's appear as Lagrange multipliers, that have to be tuned to match (1.13). This can been done thanks to (1.21). Note that the topological pressure is convex so that the solution of (1.21) is unique.

The important point is that procedure provides a unique statistical model defined by the transition probabilities (1.20). Thus, we have solved the degeneracy problem of Sect. 1.3.2.5 in the stationary case.

Range-1 Potentials

Let us now present a few examples used in the context of spike train analysis of MEA data, among others.

The easiest examples are potentials with a zero memory depth, in the stationary case, where therefore the spiking pattern ω(0) is independent of ω(-1). This corresponds to range-1 potentials.

Among them, the simplest potential has the form:

φ β (ω(0)) = N k=1 β k ω k (0) -log ( ζ β ) . (1.24)
It corresponds to impose constraints only on firing rates of neurons. We have

ζ β = N k=1 (1 + e β k ) and
the corresponding Gibbs distribution is easy to compute:

µ [ ω n m ] = n l=m N k=1 e β k ω k (l) 1 + e β k . (1.25)
Thus, the corresponding statistical model is such that spikes are independent. We call it a Bernoulli model . The parameter β k is directly related to the firing rate r k since r k = µ ( ω k (0) = 1 ) = e β k 1+e β k , so that we may rewrite (1.25) as:

µ [ ω n m ] = n l=m N k=1 r ω k (l) k (1 -r k ) 1-ω k (l) ,
the classical probability of coin tossing with independent probabilities.

Another prominent example of range-1 potential is inspired from statistical physics of magnetic systems and has been used by Schneidman and collaborators in [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF] for the analysis of retina data (Sect. 1.4). It is called Ising potential and reads, with our notations:

φ(ω(0)) = N k=1 β k ω k (0) + 1≤j<k≤N β kj ω k (0)ω j (0) -log ζ β . (1.26)
The corresponding Gibbs distribution provides a statistical model where synchronous pairwise synchronizations ω k (0)ω j (0) between neurons are taken into account, but neither higher order spatial correlations nor other time correlations are considered. The function ζ β is the classical partition function (1.23).

The Ising model is well known in statistical physics and the analysis of spike statistics with this type of potential benefits from a diversity of methods leading to really efficient algorithms to obtain the parameters β from data ( [START_REF] Tkačik | Spin glass models for a network of real neurons[END_REF][START_REF] Roudi | Statistical physics of pairwise probability models[END_REF][START_REF] Roudi | Ising model for neural data: Model quality and approximate methods for extracting functional connectivity[END_REF][START_REF] Cocco | Neuronal couplings between retinal ganglion cells inferred by efficient inverse statistical physics methods[END_REF]).

Markovian Potentials

Let us now consider potentials of the form (1.7) allowing to consider spatial dependence as well as time dependence upon a past of depth D.

Consider first a stationary Markov chain with memory depth 1. The potential has the form:

φ(ω 0 -1 ) = N -1 k=0 β k ω k (0) + N -1 k=0 k-1 j=0 0 τ =-1 β kjτ ω k (0) ω j (τ ) -log ζ β (ω(-1)). (1.27)
This case has been investigated in [START_REF] Marre | Prediction of spatiotemporal patterns of neural activity from pairwise correlations[END_REF] for spike train analysis of spike trains in the parietal cat cortex, assuming stationarity.

All examples treated above concerns stationary situations. In the non-stationary case, the entropy rate is not defined and the Gibbs distributions cannot be defined via the maximal entropy principle, while it is still possible to define them as done in Sect. 1.3.1.5. Now, the most general form for non-stationary Markov chain with memory depth D corresponds to potentials of the form:

β i1,n1,...,i l ,n l (n) ω i1 (n + n 1 ) . . . ω i l (n + n l ), (1.28)
where the sum P(N,D) holds on the set of non repeated pairs of integers (i, n) with i ∈ {1, . . . , N } and n ∈ {-D, . . . , 0}. Indeed, it can be shown that any function of blocks with depth D, f (ω n n-D ) can be written as a linear combination of all possible monomials, a polynomial, constructed on blocks of depth D [START_REF] Hammersley | Markov fields on finite graphs and lattices[END_REF]. In the non-stationary case the coefficients β i1,n1,...,i l ,n l (n) depend explicitly on n. They can be chosen so that (1.28) is normalized.

Non-Markovian Potentials

One can conceptually extend the definition of Markovian potentials to the case when D → -∞. This corresponds to a process with an infinite memory called "chain with complete connections" which is widely studied in the mathematical and mathematical physics literature (see [START_REF] Maillard | Introduction to chains with complete connections[END_REF] for a review). In this case the potential is a "function" φ n ω n -∞ depending on a infinite past ω n-1 -∞ . Although this case seems rather abstract, it turns out that the only known examples where spike train statistics can be exactly characterized in neural networks models are of this form. An example is given below. Moreover, this potential form allows to recover all the examples above.

How Good is the Estimation?

Once we have chosen a set of constraints and once we have found the parameters β to match (1.13), how can we check the goodness of the model? Additionally, changing the set of constraints provides another model. How can we compare two statistical models?

There is a wide literature in statistics dealing with the subject. In the realm of spike train analysis an important reference is [START_REF] Rieke | Spikes: Exploring the Neural Code[END_REF] and references therein. Here, we point out two criteria for model comparison, used in this chapter as an illustration.

A first and straightforward criterion consists of computing the empirical probability of blocks of depth 1, 2, . . . and to compare it to the probability predicted by the model. Of course, the number of blocks of depth R increases like 2 N R ; moreover the probability of large blocks is expected to decrease fast with the block depth. So, practically, one considers a subset of possible blocks. The challenge is in fact to predict the probability of events which have not been included as constraints in the Gibbs potential of the model.

For example, does an Ising model well predict the probability of occurrence of triplets, quadruplets, of non simultaneous spikes?

The typical representation of this criterion is a graph, with, on abscissa, the observed probability of blocks and, on ordinate, the predicted probability. Thus, to each block corresponds a point in this two-dimensional graph. A "good" model is such that all points spread around the diagonal y = x. The distance to the diagonal is expected to increase as the probability of the block decreases thanks to the central limit theorem. Indeed, if the exact probability of a block is P , then the empirical estimation of this probability is a random variable with a Gaussian distribution, of mean P , and a variance that can be computed from the topological pressure. A usual approximation of this variance is P (1 -P ). Thus, in a similar way as in Fig. 1.6, the set of points in the graph spreads around the diagonal in a region delimited by the curves ±3 P (1-P ) T called "confidence bounds". An example is given in Fig. 1.7.

Another criterion is provided by the Kullback-Leibler divergence (KL) which provides some notion of asymmetric "distance" between two probabilities. Its computation is numerically delicate but, in the present context of Gibbs distributions, the following holds. If µ is the hidden (time-translation invariant probability) and µ β a Gibbs distribution with a potential φ β , one has, [START_REF] Keller | Equilibrium States in Ergodic Theory[END_REF][START_REF] Chazottes | Pressure and Equilibrium States in Ergodic Theory, chapter Ergodic Theory[END_REF]: [START_REF] Vasquez | Gibbs distribution analysis of temporal correlation structure on multicell spike trains from retina ganglion cells[END_REF]. The estimated block probability versus the observed block probability for all blocks from range 1 to 4 (coded by colors), for N = 4 neurons with a model of range R = 3 for pairs and triplets. We include the equality line y = x and the confidence bounds (black lines) for each model, corresponding to π (T ) (w) ± 3σw with σw being the standard deviation for each estimated probability given the total sample length T ∼ 3 • 10 5 . In the figure h t corresponds to h, Eq. (1.30).

d (µ, µ β ) = P(φ β ) -µ [ φ β ] -h(µ). ( 1
which can be computed from the raster. Additionally, the entropy h(µ) is unknown and its estimation by numerical algorithms for a large number of neurons is delicate [START_REF] Strong | Entropy and information in neural spike trains[END_REF]. However, when considering two statistical models µ β1 , µ β2 with potentials φ β1 , φ β2 to analyze the same data, h(µ) is a constant (it only depends on data). Thus, comparing these two models amounts to comparing

P [φ β1 ] -π (T )
ω [ φ β1 ] and

P [φ β2 ] -π (T ) ω [ φ β2 ]. Thus, the quantity h [φ] = P [φ] -π (T ) ω [ φ ] , (1.30)
provides a relative criterion to compare models, i.e., determining if model φ β2 is significantly "better" that model φ β1 , reduces to the condition:

h [φ β2 ] ≪ h [φ β1 ] . (1.31)
Its computation is detailed in [START_REF] Vasquez | Entropy-based parametric estimation of spike train statistics[END_REF][START_REF] Vasquez | Gibbs distribution analysis of temporal correlation structure on multicell spike trains from retina ganglion cells[END_REF].

Using Gibbs Distributions to Analysis Spike Trains Statistics

In this section we show how the statistical tools presented in this chapter can be used to analyze spike trains statistics. In the "challenge" section we mention the current controversy about the question: Are G cells sensors independent encoders or, on the opposite, are neural correlations important for coding?

We present here recent works where Gibbs distributions has been be used to address this question with important implications on neural coding. However, as we show, those examples also raise additional and fundamental questions. Some of them that can be addressed on theoretical grounds, by studying neural networks models. A third section presents an example of such a model where spike trains is known to have a Gibbs statistics and where the potential is explicitly known. We compare those results to the current state of the art in spike train analysis with Ising distributions.

Are Ganglion Cells Independent Encoders?

This question can be now reformulated in the context of Gibbs distributions. Independence between neurons means that spike statistics is described by a potential, possibly non-stationary, of the form:

1 Neural networks dynamics. xxv φ n ( ω ) = N k=1 φ n,k ( ω k ) . (1.32)
This assumption can be stated independently on the memory depth, so we write here φ n (ω) instead of φ n (ω n-D n ) to alleviate notations and to be as generic as possible. In (1.32) ω k is the spike train { ω k (l) } l≤n produced by neuron k only. In this way the transition probabilities of the global network are products of transition probability for each neuron and the Gibbs distribution (1.5) is a product of marginal distributions for one neuron. On the opposite, if one believes that spike correlations play an important role in statistics one has to include them in the Gibbs potential. Typically, spike correlations are characterized by monomials and the potential takes the generic form (1.28). Obviously, there are many possible choices for this potential, depending on the set of observables assumed to be relevant.

To compare different models one can use the criteria described in Sect. 1.3.2.12. Does an independent model predicts correctly the spike blocks occurring in the observations? If not, which correlations has to be included? How evolves the Kullback-Leibler divergence as the type of correlations (monomials) taken into account growth?

However, the application of those criteria is delicate in experimental data, taking into account the large number of cells, their different types and the relatively small size of spike train samples obtained from experiments. For these reasons analysis of retina data has been performed either for memory-less models where the number of neurons can be up to 100 neurons, or to models with memory with small range and a small number of neurons. Let us present some of those works.

Weak-Pairwise Correlations Imply Strongly Correlated Network

States in a Neural Population.

A breakthrough in spike train analysis has been made by the seminal paper of E. Schneidman and collaborator [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF]. They have considered carefully instantaneous pairwise spike correlations in experiments on the vertebrate retina. It appears that, mostly, those correlations are weak. But are they significant? In Sect. 1.3.2.4 we have shown that weak correlations can be hardly interpretable without further analysis.

The maximal entropy principle provides a way to quantify the role of those correlations. Take a statistical model (a Gibbs potential) without constraint on instantaneous pairwise spike correlations, and fix only constraints on rates. This provides a Bernoulli model with a potential (1.24). Take another model where pairwise correlations are taken into account (but no higher-order interactions). This is an Ising model with potential (1.26). The authors compare those two models using the comparison between estimated block probability versus the observed block probability for spatial blocks, but also some other criteria such as the probability of having n spiking G cells in a time window of 20 ms.

Their work convincingly shows that although pairwise correlations are weak, they are nevertheless necessary to characterize spike trains statistics. The Ising model clearly does quite better than a Bernoulli model and the authors claim to "predict" about 90% of the multi spiking structure of a large G cells population in salamander and guinea pig. Note however that they focus on spatial pattern. This shows that weak correlations between pairs of neurons coexist with strongly collective behavior in the responses of ten or more neurons. Similarly, Shlens et al. [START_REF] Shlens | The structure of multi-neuron firing patterns in primate retina[END_REF] predict 99% of a complete ON and OFF parasol G cells population in primates with a Ising model . This suggests that the neural code is dominated by correlation effects.

On the other hand, this work does not directly unravel the importance of neural correlation for carrying information. Moreover, in [START_REF] Roudi | Pairwise maximum entropy models for studying large biological systems: when they can work and when they can't[END_REF], it has been shown that although Ising model is good for small populations, this is an artifact of the way data is binned and of the small size of the system. Additionally, it might be questionable whether more general forms of Gibbs distributions (e.g.involving more general monomials) could improve the estimation and account for deviations to Ising-model ( [START_REF] Shlens | The structure of large-scale synchronized firing in primate retina[END_REF][START_REF] Tkačik | Spin glass models for a network of real neurons[END_REF][START_REF] Ohiorhenuan | Sparse coding and high-order correlations in fine-scale cortical networks[END_REF]) and provide a better understanding of the neural code from the point of view of the maximal entropy principle [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]. Very recently, Ganmor and collaborators [START_REF] Ganmor | Sparse low-order interaction network underlies a highly correlated and learnable neural population code[END_REF] have extended the maximal entropy principle introducing higher order instantaneous spikes correlations. Triplets and so on are considered although all spikes arise at the same time. This therefore still corresponds to a memory-less model. These authors have convincingly shown that such model describes more accurately retina responses to natural images than Ising model.

In particular spatio temporal patterns where well predicted: binary words of 10 retinal G cells over 10 time steps. Note however that those data are binned (Sect. 1.5.1 for a definition and a discussion).

As a matter of fact, back to 1995, [START_REF] Martignon | Detecting higher-order interactions among the spiking events in a group of neurons[END_REF] already considered multi-unit synchronizations and proposed several tests to understand the statistical significance of spike synchronizations. A few years later, [START_REF] Martignon | Neural coding: Higher-order temporal patterns in the neurostatistics of cell assemblies[END_REF] generalized this approach to arbitrary spatio-temporal spike patterns and compared this method to other existing estimators of high-order correlations or Bayesian approaches. More recently, several papers have pointed out the importance of temporal patterns of activity at the network level [START_REF] Lindsey | Repeated patterns of distributed synchrony in neuronal assemblies[END_REF][START_REF] Villa | Spatiotemporal activity patterns of rat cortical neurons predict responses in a conditioned task[END_REF][START_REF] Segev | Hidden neuronal correlations in cultured networks[END_REF]. As a consequence, a few authors ( [START_REF] Marre | Prediction of spatiotemporal patterns of neural activity from pairwise correlations[END_REF][START_REF] Amari | Information geometry of multiple spike trains[END_REF][START_REF] Roudi | Mean field theory for non-equilibrium network reconstruction[END_REF]) have attempted to define time-dependent Gibbs distributions on the base of a Markovian approach, but with one time step memory only, and under assumptions such as detailed balance, [START_REF] Marre | Prediction of spatiotemporal patterns of neural activity from pairwise correlations[END_REF] or conditional independence between neurons, see Eq. ( 1) in [START_REF] Roudi | Mean field theory for non-equilibrium network reconstruction[END_REF]. A more general method relying on Gibbs distributions has been proposed in [START_REF] Vasquez | Entropy-based parametric estimation of spike train statistics[END_REF], and applied to the same data as [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF],

in [START_REF] Vasquez | Gibbs distribution analysis of temporal correlation structure on multicell spike trains from retina ganglion cells[END_REF]. The results shows a clear increase in the model accuracy (measured with the KL divergence 1.29) when adding spatio-temporal constraints.

In some sense, it is clear that the more parameters (or constraints) the more the model fits the data, up to some limit, fixed by the data and especially the raster size, where models become indistinguishable. Is there a minimal statistical model? A related question is: what does a statistical model teach us about the underlying neural network (e.g., the retina) and about neural coding? Let us first present an important work addressing the second question, before addressing the first one.

The Architecture of Functional Interaction Networks in the Retina

A critical question to elucidating the neural code, at least from a theoretical point of view, is to confront multi-electrodes real data against sophisticated statistical models testing for the underlying neural structure involved. Since the connectivity of any neural network can in principle growth exponentially as a function of the number of neurons, classical numerical methods can rapidly become inefficient. The work of Schneidman and collaborators [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF] has nevertheless suggested that although the number of possible activity patterns and underlying interactions is exponentially large in a neural network, a pairwise-based

Ising model gives a surprisingly accurate description of neural population activity patterns. So, an economical assumption to reducing a putative network dimensionality, is to use a pairwise correlation model able to recover the main network structure involved in the encoding of a natural movie [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF].

More recently, searching for further reduction on the dimensionality of the network structure, Ganmor et al. [START_REF] Ganmor | The architecture of functional interaction networks in the retina[END_REF] have shown the presence of small groups of neurons having strong correlated activity. As an outcome of their MEA spike train analysis the authors introduce the notion of functional connectivity.

This corresponds to associating with the parameters in the Gibbs potential a network of "effective" interactions (β ij for spike pairs, β ijk for triplets, . . . ).

The performance of their model is clearly higher than when assuming independence between neurons. It is able to predict the neural activity, including synchrony, for larger neural networks as in [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF] and the contribution of small groups of neurons to the complete population activity. To further reduce the structure of the network to most critical neural interactions, the authors apply a nearest-neighbors paradigm. An interesting result is that the reduction of close to 50% of the original pairwise interactions does not change the accuracy of predictions, and models most functional groups of nearest neurons with an accuracy > 95%. Additionally, small functional overlapping units (10-20 neurons) seem to be a critical structure for the encoding of natural movies stimulus.

This work, together with [START_REF] Ganmor | Sparse low-order interaction network underlies a highly correlated and learnable neural population code[END_REF], shows that a Gibbs potential with a relatively small number of parameters, corresponding to effective interactions between neurons, is able to reproduce spike trains of populations of neurons in the retina submitted to natural images. The network of these effective interactions is organized in a hierarchical and modular manner so that large network models can be constructed from smaller sub-networks in a modular fashion. Thus, by a suitable scaling of parameters, one could be able to extrapolate the Gibbs potential of a small population of neurons to large populations. Moreover, in some sense, this effective network "underlies the code", from the terminology of the authors. This means that the spike generation, as a response to a stimulus (an image), results from a dynamical process which can be summarized by the Gibbs potential of the model.

This work raises however several questions. First, the potential considered by the authors is memoryless. No time dependent process takes place in the potential. In some sense, the time-causality expected in a neural network is hidden by the effective potential. Another critical aspect is the interpretation of the effective interaction. It is stated in [START_REF] Ganmor | The architecture of functional interaction networks in the retina[END_REF] that "although the pairwise interactions in the model do not necessarily reflect a physical interaction between cells, they give a unique functional interaction map between the neurons in the network and represent statistical dependence between pairs of units." But, if xxvii they do not represent physical interactions (synapses), what do these functional interactions correspond to? To our knowledge this question has not been yet resolved.

Spike Train Analysis in a Neural Network Model

The maximal entropy principle relies on the assumption of stationarity as well as an a priori and somewhat ad hoc choice of observables. This choice severely constrains the form of the statistical model and the information that can be extracted about the underlying neuronal network producing the spike. In particular, the choice of observables determines the transition probabilities and implicitly fixes a causal structure to analyze spike events. Especially, memory-less models focuses on synchronous events, hiding somewhat temporal causality.

Obviously, it is extremely difficult to obtain a clear cut characterization of spike trains statistics from experimental data, taking into account the experimental set up, spike acquisition, spike sorting, but also the relatively small size of spike trains (typically, in retina experiments T < 10 5-6 ). So, a natural question is: "Can one have a reasonable idea of what spike statistics could be in neural network model "?

The answer is "yes".

In neural networks spikes result from the collective and non linear dynamics of neurons coupled by synapses (electrical or chemical) and submitted to "external" stimuli. As a consequence statistics of spike train is closely related to this network structure (neurons and synapses characteristics) and to the stimulus. The idea is to show here the relationships between the neural network structure and the form of transition probabilities in an explicit example, a conductance-based neural network model of Integrate-and-Fire (IF) type called "generalized Integrate and Fire" (gIF) and introduced by Rudolph and Destexhe [START_REF] Rudolph | Analytical integrate and fire neuron models with conductance-based dynamics for event driven simulation strategies[END_REF]. This section summarizes the paper [START_REF] Cessac | Statistics of spike trains in conductance-based neural networks: Rigorous results[END_REF].

The Model

We consider the evolution of a set of N neurons. Here, neurons are considered as "points" instead of spatially extended and structured objects. As a consequence, we define, for each neuron k ∈ {1 . . . N }, a variable V k (t) called the "membrane potential of neuron k at time t" without specification of which part of a real neuron (axon, soma, dendritic spine, ...) it corresponds to.

Fix a firing threshold θ. The sub-threshold dynamics of the model is: for some r. Then,

C k dV k dt + g k ( t, ω ) V k = i k (t, ω). ( 1
g k ( t, ω ) = g L,k + N j=1 G kj α kj ( t, ω ) , (1.34)
where α kj ( t, ω ) = t (r) j (ω)<t α kj tt Returning to Eq. (1.33) the term i k (t, ω) is a current given by i

k (t, ω) = g L,k E L + N j=1 W kj α kj ( t, ω ) + i (ext) k (t) + σ B dB k dt
, where E L is the Nernst leak potential, W kj is the synaptic weight from j to k, i As in all IF models, when V k reaches the firing threshold θ, it is reset. Here, however, it is not instantaneously reset to a constant value, but to a Gaussian random variable with mean zero and variance σ 2 R , after a time delay δ > 0 including the depolarization-repolarization and refractory period. We call τ k ( t, ω ) the last time before t where neuron k has been reset.

Membrane Potential Decomposition

Given the spike history of the network it is easy to integrate (1.33) and to find an explicit expression for the membrane potential at time t. It depends on the past spike history ω. Denote V k (t, ω) the membrane potential at time t given the spike history before t. Set:

Γ k (t 1 , t 2 , ω) = e -1 C k R t 2 t 1 g k ( u,ω ) du , (1.35)
corresponding to the flow of (1.33). We have

V k (t, ω) = V (det) k (t, ω) + V (noise) k (t, ω), V (det) k (t, ω) = V (syn) k (t, ω) + V (ext) k (t, ω) (1.36)
where

V (syn) k (t, ω) = 1 C k N j=1 W kj t τ k ( t,ω ) Γ k (t 1 , t, ω) α kj ( t 1 , ω ) dt 1 , (1.37)
is the synaptic interaction term which integrates the pre-synaptic spikes from the last time where neuron k has been reset. Additionally,

V (ext) k (t, ω) = 1 C k g L,k E L + t τ k ( t,ω ) i (ext) k (t 1 ) Γ k (t 1 , t, ω)dt 1 , (1.38)
contains the stimulus i (t, ω) is the stochastic part of the membrane potential. This is a Gaussian variable with mean zero and variance

σ 2 k ( t, ω ) = Γ 2 k (τ k ( t, ω ) , t, ω) σ 2 R + σ B C k 2 t τ k ( t,ω ) Γ 2 k (t 1 , t, ω) dt 1 . (1.39)
The first term in the right-hand side comes from the reset of the membrane potential to a random value after resetting. The second one is due the integration of synaptic noise from τ k ( t, ω ) to t.

Statistics of Raster Plots

It has been shown in [START_REF] Cessac | Statistics of spike trains in conductance-based neural networks: Rigorous results[END_REF] that the gIF model (1.33) has a unique Gibbs distribution with a non-stationary potential:

φ n ( ω ) = N k=1 φ n,k ( ω ) , (1.40)
with,

φ n,k ( ω ) def = ω k (n) log π ( X k (n -1, ω) ) + ( 1 -ω k (n) ) log ( 1 -π ( X k (n -1, ω) ) ) . (1.41)
Here As announced in the beginning of this section, the Gibbs potential that constraints spike statistics summarizes several contributions. The effect of the synaptic network structure, integrated over the past spike history as well as over the contribution of the stimulus (external current), also integrated over the past spike history, appears in the term V (det) k

X k (n -1, ω) = θ -V (det) k (n -1, ω) σ k (n -1, ω) , (1.42
(n -1, ω). The synaptic noise and the reset to a random value, integrated over the history, depends on the term σ k (n -1, ω). We insist that the result holds for a time-dependent external current, i.e., a non-stationary dynamics.

The potential exhibits therefore clearly the causal structure of spikes generation. The probability to have a spike at time n depends explicitly on synaptic weights, defining the neural network structure, and on the stimulus, via the external current. Moreover, it introduces a clear history dependence. Now, there are several important remarks.

1. The number of parameters on which spike train statistics depend is relatively small. It increases like a polynomial in the number on neurons, e.g., at most N 2 synaptic weights.

2. Those parameters are physical parameters: the synaptic weights between neurons, the stimulus, plus physical characteristics of the neuron such as leak Nernst potential. Thus, the potential outlines a network which is not effective, as in Sect. 1.4.3, but is the real underlying network. Additionally, this potential is non-linear function of the synaptic weights and stimulus while Ising-like models are linear in the interactions.

3. The memory is mathematically infinite. In other words, to handle the spike statistics exactly, one has to consider a non-Markovian potential with an infinite history (D → -∞).

From Non-Markovian to Markovian Potentials

Since the dependence on the past decays exponentially fast, thanks to the exponential decay of synaptic response, it is possible to provide Markovian approximations of the potential (1.40). Basically, one truncates the synaptic response after a characteristic time D and performs a series expansion of the function (1.43), using the fact that the power of a monomial is the same monomial. So, the series becomes a polynomial, which provides a Markovian potential of the form (1.28). Here, the coefficients β i1,n1,...,i l ,n l (n)'s depend explicitly on the synaptic weights (network structure) as well as on the stimulus. Now, for N neurons and a memory depth D, the truncated potential contains 2 N D coefficients β i1,n1,...,i l ,n l (n), while the exact (D → -∞) potential depends only on a polynomial number of parameters. This shows that, in this model, a potential of the form (1.28) induces a strong, and somehow pathological, redundancy.

Additionally, the truncated potential is far from Ising, or more elaborated potentials used in experimental spike train analysis. As we have seen most of these models are memory-less and non causal. This model is not sufficient to describe the retina since it neglects the specific types of bipolar, horizontal and some amacrine cells that do not "fire". Additionally, it neglects electric synapses (gap junctions)

playing an important role in the connectivity as shown in Fig. 1.1(b). Recent investigations show that the conditional factorization property (1.40) disappear in the presence of gap junctions, so that statistics is expected to be even more complex, with no independence at all (Cessac, Cofre, in preparation).

Conclusion

In this chapter we have attempted to give a short overview of recent questions related with the concept of spike train analysis. corresponds to an integration of physical interactions over the binning time scale. For example, in the Ising model, the pairwise coefficient β kj integrates the effect of several circuits connecting neuron j to neuron k, as well as dynamic-dependent effects. As a matter of fact its interpretation is rather delicate. This is however certainly not the end of the story and this aspect has to be still investigated on the theoretical and experimental side.

Linear Potentials versus Combinatorial Explosion

Experimental attempts to go "beyond Ising" [START_REF] Marre | Prediction of spatiotemporal patterns of neural activity from pairwise correlations[END_REF][START_REF] Vasquez | Parametric estimation of spike train statistics by gibbs distributions : an application to bio-inspired and experimental data[END_REF] suggest that Markovian models with increasing range should describe better and better the spike statistics. This is also expected from the theoretical analysis summarized in Sect. 1.4.4. However, this raises several remarks and questions. First, it is evident that the more parameters, the best is the model, but what do we learn from this plethora of parameters? Second,

1 Neural networks dynamics.

xxxi one has to face the critical issue of an exponential increase of parameters, with the potential range and with the number of neurons, so that numerical methods can rapidly become inefficient. Moreover, the sample size required to determine those coefficients is expected to increase also exponentially, ruining any hope to extract reliable coefficients from empirical data. Finally, as outlined in the previous section, the interpretation of the coefficients is difficult even for the simplest pairwise interactions β kj .

Our point of view is that the linear potential approach, based on the maximal entropy principle, is maybe inappropriate. On the opposite, non linear potentials of the form (1.40), truncated to a finite memory depend on a number of parameters, the physical parameters of the network, which increases only polynomially with the number of neurons. Although, the number of blocks determining the potential increases exponentially fast with the memory depth, it could well be that only a small proportion of blocks are sufficient to extract most of the information about the hidden parameters. Finally, the interpretation of parameters is here straightforward and such a model allows to treat the non-stationary case. This may provide an alternative to Ising-like models to study spike train statistics in experiments.

Outlook

In this section we would like to point out a few challenges for the next future, on the theoretical and experimental sides.

Gibbs Distributions and the Neural Code

Although interesting results have come out from the Gibbs analysis of retina spike trains, the link between spike statistics and their modeling with Gibbs distribution on one hand, and the way how a visual scene is encoded by G cells spikes emission on the other hand, remains rather tight. Defining a priori a Gibbs potential from a set of constraints superimposes upon spike trains a causal structure, purely spatial or spatio-temporal, associated with "effective interactions", e.g.the coefficients β i1,n1,...,i l ,n l (n) in the polynomial expansion (1.28). What do these interactions teach us about the neural code? How are they related to a visual scene? Given a Gibbs distribution that fits well retina spike statistics is it possible to infer information about the visual scene perceived by this retina? Is it possible to build retina "decoders" from Gibbs statistics? If yes, does a "correlated decoder", with correlated G cells perform better than a "rate decoder" where G cells outputs are independent? Although interesting advances has been done on these questions (see, e.g., [START_REF] Schaub | The ising decoder: reading out the activity of large neural ensembles[END_REF]) we believe that they are far from being solved, and that they constitute a challenge for the next years.

A related question concerns the concept of receptive field. We have presented in Sect. rates without need of considering higher order statistics. It is a current approach to model RF as filters, typically linear, followed by a non-linear correction [START_REF] Pillow | Prediction and decoding of retinal ganglion cell responses with a probabilistic spiking model[END_REF]. How does the paradigm of linear-non linear filter connects with the paradigm of Gibbs distribution? Can we infer the shape of the RF filter from the Gibbs potential? Classical RF filters are based on firing rates for individual G cells; on the opposite Gibbs statistics deals with spike events (monomials) for collective behaviors. Are these two visions coherent? Is there a link e.g.between effective interactions and RF?

To our best knowledge those questions have not been resolved. On the theoretical side they can be addressed in the context of Gibbs distributions. Given a Gibbs potential modeling retina response it could be possible to compute the linear response to a stimulus considered as a weak perturbation of dynamics.

This linear response is characterized by a convolution kernel which could be compared to the models of receptive field filters used in the literature. This work remains still to be done though.

Experimental Limits

To our opinion the current experimental set up is faced with the following limits.

• Natural stimulus must reproduce ecological environment, where animals lives, including the way how animals explore it, how they are in action, moving their head and eyes. As a consequence the captured images used in experiments need to be dynamically displayed to the retina, reproducing natural motions.

• Sophisticated experiment require MEA devices recording from a large numbers of cells (e.g.> 100).

For example a new MEA technology using 4096 electrodes matrix makes it possible to recording most of the neurons in a single small surface, but the recording data string would take several gigabytes of saving space and requires fast technologies to access data. Moreover, their numerical analysis, spike sorting + spike train analysis, in any computer or even a cluster will take a very long time. Actually,

• The MEA outputs need adequate spike sorting algorithms able to deal with larger and larger numbers of cells. Current algorithms allow to treat about 200 electrodes signals.

• Spike train analysis necessitates adequate statistical tools applying to a large number of interacting cells to evaluate different possible models for neural encoding (e.g.population coding). Current algorithms allow to treat less than 100 neurons for an Ising model.

• A validation of any neural coding model is required. This is done by contrasting its performance against real behavioral sensory results for the animals under study [START_REF] Jacobs | Ruling out and ruling in neural codes[END_REF]. Additionally, without the precise quantification of the animal performance for a particular behavioral task, responding to natural stimulus, it will not be possible to access the extended validity of any proposed model. Thus, both the animal capacity and the theoretical model need to be contrasted.

Clearly, those constraints constitute high level challenges for the scientific community. Probably, this is only the beginning of a long story.

1. 2

 2 Unraveling the Neural Code in the Retina via Spike Train Statistics Analysis 1.2.1 Retina Structure and Functions 1.2.1.1 Retina Structure

Fig. 1 . 1 Fig. 1 . 2

 1112 Fig. 1.1 Processing steps of the visual stream. (a) The cellular organization of the retina (from Expert Reviews in Molecular Medicine by Cambridge University Press 2004); (b) Main connectivity structure between retina cells types (from [24])

Fig. 1 .

 1 Fig.1.3(a) explains in a schematic way how this property results from the connectivity between P, B, and H cells. In the example, the illumination of the photoreceptors in the center of the RF results in a depolarization of ON B cells so in an increase of spikes rate in the respective connected G cells. On the opposite, the illumination of the photoreceptors at the periphery of the RF results in a hyper-polarization of OFF B cells so in a decrease of spikes rate in the respective connected G cells. In more general terms, as a consequence of this architecture (Fig.1.3(a)), a G cell connected to that B cell fires spikes at the maximal rate when the center of the RF is illuminated and when the surround is dark (Fig.1.3(b-2), case 3). On the opposite it fires no spike at all when the center of the RF is dark and the surround is illuminated (Fig.1.3(b-4) case 4). Fig. 1.3(b) summarizes the different patterns of illuminations -G cell response in terms of spike firing and the functional implication of RF organization. For example, a full, uniform, illumination of the RF

  Fig. 1.3 Center-Surround antagonism (a) Illumination of a piece of retina (from http://www.webexhibits.org/colorart/ganglion.html ). (b) ON-center and OFF-center RF, Figure from [1]. The first line shows center-surround architecture of the cell while lines 2-6 shows a typical response of the G cells and the illumination pattern leading to that response. (b-1) Center-surround architecture of an ON-center cell and illumination pattern. (b-2) Time duration of the stimulus and spike response of the cell. Time is in abscissa. (b-3) and (b-4) Same as columns (b-1) and (b-2) for an OFF-center cell. Case 1 left (right) is a ON-center (OFF-center) G cell where a light spot (yellow) in the center of the RF generates an increase (decrease) of spike firing. In case 2 a spot stimulus in the surround generates a decrease (increase) of the spike rate. In case 3,4 an increase in the size of the stimuli leads a sharper response. In case 5 a diffuse stimulus covering the center-periphery has no effect on the spike firing rate.

Fig. 1 . 4

 14 Fig. 1.4 Detection of differential motion. (a) An object-motion-sensitive G cell remains silent under global motion of the entire image but fires when the image patch in its RF moves differently from the background. (b) Scheme summarizing the circuitry behind this computation. Rectification (see[START_REF] Demb | Cellular basis for the response to second-order motion cues in y retinal ganglion cells[END_REF] for a description of rectification mechanism.) of B cell signals in the RF center creates sensitivity to motion. Polyaxonal A cells in the periphery are excited by the same motion-sensitive circuit and send inhibitory inputs to the center. If motion in the periphery is synchronous with that in the center, the excitatory transients will coincide with the inhibitory ones, and firing is suppressed (Fig. from[START_REF] Gollisch | Eye smarter than scientists believed: neural computations in circuits of the retina[END_REF]. The legend is adapted from this reference).

  ) is called a (normalized) Gibbs potential . The advantage of this definition of Gibbs distribution is that it holds for time-dependent transition probabilities contrarily to the classical definition from the maximal entropy principle (Sect. 1.3.2.8). Moreover, in the case (1.7) the Gibbs potential depends explicitly on time (index n). This definition also extends to system with infinite memory (Sect. 1.4.4) although Eq.(1.5)

  and so on. Why using the cumbersome notation π (T ) ω [ f ]? This is to remind the reader that such empirical averages are random variables. They fluctuate from one raster to another i.e., π (T )

ω

  ( k, j ) for different values of ω, while increasing T (Fig.1.6).

1

  Fig.1.6 Correlation (1.12) as a function of sample length T in a model where spikes are independent. For each T we have generated 1000 rasters of length T , with two independent neurons, drawn with a firing rate p = 1 2 . For each raster we have computed the pairwise correlation (1.12) and plotted it in log-scale for the abscissa (red point). In this way we have a view of the fluctuations of the empirical pairwise correlation about its (zero) expectation. The full lines represent respectively the curves 3 qp 2 (1-p 2 ) T(blue) and -3 q p 2 (1-p 2 ) T (green) accounting for the Gaussian fluctuations of C (T ) ω ( k, j ): 99% of the C (T ) ω ( k, j )'s values lie between these two curves.

Fix

  K observables O k , k = 1, . . . , K, and compute their empirical average π (T ) ω [ O k ]. The remarks of the previous sections hold: since all rasters are distributed according to µ, π (T ) ω [ O k ] is a random variable with mean µ [ O k ] and Gaussian 1 fluctuations about its mean, of order 1 √ T . If there are N > 1 rasters the experimentalist can estimate the order of magnitude of those fluctuations and also analyze the raster-length dependence. In fine, she obtains an empirical average value for each observable, π (T )

j

  (ω) sums up the spike responses of post-synaptic neuron k to spikes emitted by the pre-synaptic neuron j at times t (r) j (ω) < t. g L,k is the leak conductance of neuron k.

  -dependent) external stimulus received by neuron k, and dB k is a Brownian noise whose variance is controlled by σ B .

  ω) contains the deterministic part of the membrane potential. On the opposite, V

  1.2.1.2 the classical notion of RF which is associated with an isolated G cell, independently of what the neighboring G cells are perceiving. Additionally, e.g.fig. 1.3, describes well the response of a G cell in terms of firing

  Determining the "Best" Markov Chain to Describe an Experimental RasterWe now show how the formalism of the previous section can be used to analyze spike trains statistics in

	1 Neural networks dynamics.	xvii
	1.3.2 experimental rasters.
	1.3.2.1 Observables
	We call observable a function which associates to a raster plot a real number. Typical examples are
	f (ω) = ω

k (n) which is equal to '1' neuron k spikes at time n in the raster ω and is '0' otherwise; likewise the function f

  and only if neuron k and k ′ fire synchronously at time n in the raster ω. These two cases are example of what we call monomials in the chapter, namely functions of the form ω k1 (n 1 ) ω k2 (n 2 ) . . . ω km (n m ) which is equal to 1 if and only if neuron k 1 fires at time n 1 , . . . , neuron k m fires at time n m in the raster ω. Thus monomials attribute the value '1' to characteristic spike events. One can also consider more general forms of observables, e.g.non linear functions of spike events

(see for example Eq. (1.37), (1.41) below).

  .33) C k is the membrane capacity. g k ( t, ω ) is the integrated conductance at time t given the past spike activity encoded in the raster ω. It is defined in the following way. Denote α kj ( tτ ) the synaptic response of neuron k, at time t, to a pre-synaptic spike from neuron j that aroused at time τ . Classical examples of synaptic responses are α kj ( t ) = e where H the Heaviside function (that mimics causality) and τ kj is the characteristic decay times of the synaptic response. In gIF model the conductance g k ( t, ω ) integrates the synaptic responses of neuron k to spikes arising in the past. Call

	-t τ kj H(t) or α kj ( t ) = t τ kj e -t τ kj H(t), t (r) j (ω) the r-th spike-time emitted by neuron j in the raster ω, namely ω j (n) = 1 if and only if n = t	(r) j (ω)

  to replace this function by a constant, giving therefore a Bernoulli potential. So, Ising model as well as memory less models are rather poor in describing the statistics of model(1.33). But, then, how can we explain the success of Ising model to analyze retina data? We return to this point in the conclusion section. the model obeys a mean-field theory as the number of neurons tends to infinity. This requires, in general, strong constraints on the synaptic weights (vanishing correlations), not necessarily realistic.

	1.4.4.5 Are Neurons Independent? For this model we can answer the question of neurons independence. The potential (1.40) is a sum over neurons, similarly to (1.32), but it has not the same form as (1.32). The difference is subtle but essential. While in (1.32) the potential of neuron k, φ k depends upon the past via the spike-history ω 1.4.4.6 Links with the Retina

Now, the best approximation of the potential

(1.40)

by a memory-less potential is . . . Bernoulli. This is because of the specific form of φ: a term ω k (0) multiplied by a function of ω -1 -∞ . To have a memory-less potential one has k of neuron k only, in (1.40) φ k depends upon the past via the spike-history ω of the entire network. The factorization in (1.40) reflects only a conditional independence: if the past spike history of the network is fixed then the only source of randomness is the noise, which is statistically independent by hypothesis, for all neurons. So, there is nothing deep in the factorization (1.40). On the opposite, a factorization like (1.32) would reflect a deeper property. Neurons would somewhat be able to produce responses which are well approximated by a function of their own history only, although each neuron receives inputs from many other neurons. Considering the form of the potential φ, given by (1.41) there is no chance to obtain the strong factorization property (1.32) unless neurons are disconnected. This property could however arise if

  Taking as an example the case of the retina we have presented a summary of recent experimental progresses from MEA recording to spike train analysis. On the theoretical side, we have introduced a general formalism connecting spike train statistics to Markov chains and Gibbs distributions. This formalism looks appropriate since, on one hand it allows to recover the Gibbs distributions forms used currently in the literature of spike train analysis, and on the other hand it affords analytical developments to characterize spike train probabilities in neural networks models. Finally, we have presented three examples of recent successes in spike trains analysis. These examples are encouraging but raise salient questions that we would like now to address. if and only if neuron k as fired at least once in this slice. In this way, one smooths out the dynamical interactions occurring at a time scale smaller than τ (especially synaptic interactions). So the coefficient β kj 's in a binned-Ising model with a binning of 10 -20ms somewhat integrate the synaptic transmission effects and neurons pairwise interactions appear as instantaneous. In this way, one looses an important part of the dynamics and of the network structure. The "functional interactions" evoked in Sect.1.4.3

	1.5.1 Ising or not Ising?

In Sect. 1.4.2, 1.4.3 we have outlined the relative success of Ising model to analyze retina data, while in Sect. 1.4.4.4 we have computed explicitly the potential and concluded that it is quite far from Ising. What is the reason of this discrepancy? A first explanation, exposed in Sect. 1.4.4.6, is that the model (1.33) is not a good model for the retina. Another possible reason is the difference of time scales considered in both approaches. While the theoretical results of Sect. 1.4.4 consider neurons dynamics at the time scale of a spike (about 1 ms), statistical analysis of experimental data use, in all the examples we know, data binning. From preliminary analyzes of spike train (correllograms), one extracts a characteristic time scale τ (about 10 -20ms) from which spike trains are binned. Recall that a binned spike train is a raster Ω, obtained by cutting the original raster ω into time-slices of size τ and setting Ω k (n) = 1 in the slice n

Fluctuations are not necessarily Gaussian, if the system undergoes a second order phase transition where the topological pressure introduced in Sect. 1.3.1.5 is not twice differentiable.

Acknowledgment

This work has been supported by ERC grant Nervi 227747 (BC), European grant BrainScales (BC), ANR-CONICYT grant (KEOPS), FONDECYT 1110292 (AP) and ICM-IC09-022-P (AP).

http://webvision.med.utah.edu/ This site summarizes recent advances in knowledge and understanding of the visual system through dedicated chapters and evolving discussion to serve as a clearing house for all things related to retina and vision science.

The brain from top to bottom: the retina http://thebrain.mcgill.ca/flash/d/d_02/d_02_cl/d_02_cl_vis/d_02_cl_vis.html This web site contains a series of topics dealing with the brain: "memory and the brain", "evolution and the brain" and so on. Each topic is developed at different levels: beginner, intermediate, advanced. This is a very useful and didactic reference.

Information Processing in the Retina

http://www.sumanasinc.com/webcontent/animations/content/receptivefields.html This page contains an animation illustrating the functioning of On-Off Receptive Fields. This is part of the web site http://www.sumanasinc.com/webcontent/animations.html which contains nice animations on different topics, including neuroscience. Multi Electrode arrays http://en.wikipedia.org/wiki/Multielectrode_array From the famous web site wikipedia. Paul Avery's Home Page. Image Gallery: Vision and the Eye

http://www.phys.ufl.edu/ ~avery/course/3400/gallery/gallery_vision.html This page contains a series of nice pictures illustrating the functioning of eyes and vision.