Denis Jouvet
email: denis.jouvet@loria.fr

David Langlois
email: david.langlois@loria.fr

Mohamed Amine Menacer
email: mohamed-amine.menacer@loria.fr

Dominique Fohr
email: dominique.fohr@loria.fr

Odile Mella
email: odile.mella@loria.fr

Kamel Smaïli
email: kamel.smaili@loria.fr

About vocabulary adaptation for automatic speech recognition of video data

Keywords: Speech recognition, vocabulary, vocabulary adaptation, vocabulary selection

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The vocabulary is one of the key components of an automatic speech recognition (ASR) system. It needs to be adequate with respect to the considered speech recognition task, and this is usually achieved through a training or adaptation process. That is the object of this paper, which discusses the adaptation of vocabularies for automatic speech transcription of videos in French, English and Arabic, for AMIS (Access Multilingual Information opinionS) project 1 . AMIS project aims at helping users to access information from videos that are in a foreign language, that is to understand the main ideas of the video. The best way to do that, is to summarize the video for having access to the essential information. Therefore, AMIS focuses on the most relevant information in videos by summarizing and translating it to the user. Obviously, the process starts by an automatic transcription of the audio channel.

Baseline ASR systems used at the beginning of the project have been developed from available corpora. For what concerns the linguistic part, that means that the vocabularies and the associated language models have been elaborated from quite old text data. Consequently, the vocabularies are somewhat outdated, and they are not relevant for a proper processing of person names and locations that have recently emerged in the news. Besides the fact that out-of-vocabulary (OOV) words affect speech recognition performance (in average, each out-ofvocabulary word produces 1.2 errors [START_REF] Rosenfeld | Optimizing lexical and ngram coverage via judicious use of linguistic data[END_REF]), names of persons and of locations convey a very important and useful information for understanding the content of the videos. One way to cope with this aspect is to collect large amounts of text data over the web, that correspond to about the same time period as that of the videos to be processed, and build new speech recognition vocabularies from this new text data.

Unknows words are also problematic in natural language processing, for example for syntactic parsing and for machine translation. Several papers have investigated the handling of unknown words [START_REF] Attia | Handling unknown words in statistical latentvariable parsing models for Arabic, English and French[END_REF], including the use of a probabilistic model for guessing base forms [START_REF] Lindén | A probabilistic model for guessing base forms of new words by analogy[END_REF] in English and Finish, and a morphological guesser for lemmatization in Arabic [START_REF] Attia | The Floating Arabic Dictionary: An Automatic Method for Updating a Lexical Database through the Detection and Lemmatization of Unknown Words[END_REF]. However, such approaches for dealing with written texts are not applicable to speech recognition.

With respect to speech recognition, several approaches have been developed in the past for elaborating vocabularies that are adequate for a given task. When a single text corpus is available, and when this corpus is homogeneous, the selection method is straightforward, it simply consists in selecting the most frequent words in the training corpus. However, since many years, the selection is done from numerous and heterogeneous corpora, which differ strongly in term of source or content (e.g., various radio or TV channels, journals, speech transcripts, …), time period, and size (from a few million words up to more than several hundred million words). In such case, it is not suitable to concatenate all the text corpora and just select the most frequent words. A frequent word in a small corpus, thus interesting to select, may end up with a small frequency in the concatenated data, and would thus not be selected.

When dealing with a heterogeneous set of text corpora, various selection methods have been proposed that rely on the unigram distribution of the words in each sub-corpus. A conventional approach consists in finding the linear combination of the unigrams associated to each sub-corpus, that matches the best with the unigram distribution of some development set [START_REF] Venkataraman | Techniques for effective vocabulary selection[END_REF],and [START_REF] Allauzen | Automatic building of the vocabulary of a speech transcription system[END_REF]; then the words having the largest unigram values (according to the combined unigram distribution) are selected. The combination parameters are obtained through an expectation-maximization process. Selection approaches based on neural networks have also been investigated [START_REF] Jouvet | A machine learning based approach for vocabulary selection for speech transcription[END_REF]. It should be noted that all these techniques require the availability of a development set, representative of the task, for optimizing the unigram combination weights.

This paper investigates the selection of speech recognition vocabularies in French, English and Arabic, for the automatic transcription of videos in AMIS project. It is organized as follows. Section 2 presents the baseline speech recognition systems. Section 3 describes the collection of the textual data over internet. Section 4 presents an analysis of the collected data, with a comparison to the GigaWord data sets. Finally, Section 5 details the selection of speech recognition vocabularies and discusses some evaluation results.

Baseline ASR systems

The speech recognition systems are based on the KALDI speech recognition toolkit [START_REF] Povey | The kaldi speech recognition toolkit[END_REF].

Acoustic modeling is based on Deep Neural Networks (DNN), as such modeling provides the best performance [START_REF] Hinton | Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups[END_REF]. The DNN has an input layer of 440 neurons (11 frames of 40 coefficients each), 6 hidden layers of 2048 neurons each, and the output layer has about 4000 neurons, corresponding to the number of shared densities of the initial GMM-based speech recognition system. The classical n-gram approach is used for language modeling. 1 presents some characteristics of the baseline ASR systems, with respect to some linguistic aspects. Vocabulary sizes vary from about 100 k words (for French and Arabic) to 150 k words (for English). The average number of pronunciations variants vary from 1.1 for English, to 2.1 for French, and 5.1 for Arabic. In French, most of the pronunciation variants are due to the optional mute-e at the end of many words, and to possible liaison consonants with following words starting by a vowel. In Arabic, the larger number of pronunciation variants is due to the absence of diacritic marks, which indicate short vowels, in the spelling of the vocabulary words.

Web textual data

As the vocabularies in the baseline ASR systems haven been defined according to available text corpora, that are rather old, the vocabularies are somewhat outdated, and they do not properly reflect the names of persons and locations observed in the recently collected videos of AMIS project. To update the vocabularies, new text data has been collected over the internet, in a period matching the period of the videos. This section also describes the elaboration of the test and development sets.

Training corpus

A few newspaper, radio and TV web sites in French, English and Arabic have been selected for collecting text data. A script was used to crawl web pages from the given sites over several months. The period over which text data was collected, was the same for the three languages.

A preprocessing has been applied on the raw text data collected from the various web sites. It mainly consists in removing useless data (e.g., date tags, hour tags, some keywords such as "view image", "download", ...), long non-Arabic text in Arabic web pages, ... Moreover, all duplicated sentences were also removed. About 80% of the amount of collected data is thus ignored. The amount of word occurrences available per language, after this preprocessing, is reported in

Test corpus

The videos processed in AMIS project have been collected from Youtube. They correspond to various channels such as Alarabiya, Alquds, BBC, EnnaharTV, Euronews, France24, RT, SkynewsArabia… For most of the videos, Youtube provides short descriptions which correspond to the content of the video, and thus contain names of persons and locations occurring in the videos. Such data has been collected for all AMIS videos (a few thousand videos per language). This data is used as a test data set for evaluating the percentage of occurrences of out-of-vocabulary words. Table 3 indicates the amount of word occurrences in the development sets, for each language. An example of YouTube description is available in the top part of Figure 1.

Development corpus

On Euronews web site, one can find descriptions of Euronews videos. Such descriptions generally provide detailed information on the content of the video, which in some cases, is rather similar to a transcription of its content. Independently of the collection of videos for AMIS project, another set of about 8000 videos, in Arabic language, were collected from Euronews web site. Cross language links available in the Euronews descriptions make possible to collect also the descriptions in French and in English for those videos. This led to about 8000 text descriptions in French, in English and in Arabic. This data set, which is not associated to AMIS videos, but comes from a similar period was used as a development set for the selection of the new vocabularies.

Among the videos collected in AMIS project, a part of them corresponds to Euronews. Hence, you can find in the top part of Figure 1 an example of a Euronews description (long description), along with the YouTube description (which is much shorter, four lines only).

Analysis of the collected web data

An analysis of the data was carried out. For each data set collected from internet (i.e., French, English and Arabic), the frequency of occurrences of the words has been analyzed. The same analysis was applied on the GigaWord corpora available from the LDC (French [START_REF] Graff | French Gigaword Third Edition LDC2011T10[END_REF], English [START_REF] Graff | English Gigaword LDC2003T05[END_REF], and Arabic [START_REF] Parker | Arabic Gigaword Fifth Edition LDC2011T11[END_REF]).

As a result, Figure 2 displays, for each corpus, the coverage of the word occurrences with respect to the most frequent word tokens of the corresponding corpus. For example, for English with data collected over internet the 100,000 most frequent words cover about 96.6% of the 2,880 million word occurrences of the English data; whereas for Arabic, the 100,000 most frequent words cover only 92.7% of the 690 million word occurrences of the Arabic data. Solid lines correspond to the data collected on internet. Dotted lines correspond to the GigaWord corpora.

Figure 2: Coverage of text data with respect to the most frequent words (of each data set).

On the figure, the 4 curves corresponding to "French (gigaword corpus)", "French (internet data)", "English (gigaword corpus)", and "English (internet data)" are very similar. For each language, internet data and GigaWord data leads to very similar results. The figure also shows that to reach a given coverage, much more words are needed in Arabic than in French and English languages, probably due to the morphological richness of Arabic.

Updated vocabularies

The new text data collected over internet is used here as text material from which new ASR vocabularies are selected for transcribing AMIS videos. Results are then analyzed mainly in terms of percentage of out-of-vocabulary words in the test sets for the different languages and different vocabulary sizes.

Selection of vocabulary words

The selection process relies on the conventional approach. First a unigram is estimated on each subset of the training corpus, corresponding to a radio channel, a TV channel, a journal, etc. For example, for the French language, the various subsets correspond to Euronews, France 24, France Inter, Le Monde, Le Figaro, L'Humanité, and so on. Overall, there are about 30 subsets for the French language. A similar splitting, according to web site, is done also for English and Arabic data, leading to 22 subsets for English and 29 subsets for Arabic.

Once unigrams models are trained on each subset, they are linearly combined to make a global unigram. The weights of the linear combination are estimated with an Estimation-Maximization (EM) algorithm to match as best as possible the unigram estimated on the development data. The objective function that the E.M. algorithm optimizes is the Kullback-Leibler distance between the unigram distribution corresponding to the linear combination of the unigrams estimated on each sub-corpus, and the unigram distribution estimated on the development corpus. On the French data, the two largest combination weights and the associated sub-corpus are the following: 0.876 for Euronews, and 0.106 for France24. All the other weights are below 0.01. A similar behavior is observed for the other languages. The large weight obtained for the Euronews channel may be due to the fact that the development set is made of descriptions of Euronews videos.

Finally, the selected vocabulary corresponds to the words that have the largest probability in the combined unigram. For each language, four vocabularies have been extracted corresponding respectively to the 100 k, 200 k, 400 k and 800 k most probable words.

Analysis of results

The best analysis that could be carried out requires a manual transcription of a large subset of AMIS videos. As such transcription is not available for AMIS videos, we have used the text data corresponding to the Youtube descriptions as test sets. On the test sets, we evaluated the amount of out-of-vocabulary words for the various vocabularies: baseline ASR vocabulary, and new vocabularies (100 k, 200 k, 400 k, and 800 k words). Results for the 3 languages are reported in Table 5. For comparison purpose, Table 4 reports the percentages of out-ofvocabulary words on the development sets. As can be seen on these tables, the percentage of out-ofvocabulary words is much lower with the new vocabularies than with the old ones. The same behavior is observed on the development and on the test sets. In all cases, increasing the size of the vocabularies significantly reduces the percentage of outof-vocabulary words. For example, for the English data, on the test set, the OOV rate was reduced from 5.5% with the baseline vocabulary (150 k words) to 3.3% with the new 100 k word vocabulary, and then to 2.7%, 1.9% and 1.5% respectively with the 200 k, 400 k and 800 k vocabularies.

Comparing the languages, the OOV rates are smaller for the French data than for the English data. The largest OOV rates are observed for the Arabic language. The large OOV rate on Arabic data was also observed in other studies related to statistical modeling of Arabic [START_REF] Meftouh | Comparative Study of Arabic and French Statistical Language Models[END_REF] and [START_REF] Meftouh | Modeling Arabic Language using statistical methods[END_REF].

To check the benefit of the new vocabularies, they have been used for a new transcription of AMIS videos. The two speech recognition results obtained with the old vocabulary, and with the new vocabulary (100 k words), are displayed as simultaneous subtitles. Figure 1 provides a typical example of the recovery of names of persons thanks to the new vocabularies. "John Kerry" was not present in the old French vocabulary, and thus the corresponding occurrence was replaced by a sequence of short words which are acoustically close. As the person name is missing in the old vocabualry, the corresponding transcription (cf. line ASR-V01 in Figure 1) gets difficult to understand, and an important information (the name "John Kerry") is missing; such behavior will also impact the machine translation process. With the new vocabularies, this problem is overcome.

Conclusion

This paper has investigated the problem of out-of-vocabulary word in the transcription of videos in French, English and Arabic. A large part of out-of-vocabulary words concerns names of persons and locations, which convey an important information for understanding the content of videos. To elaborate speech recognition vocabularies that are adequate for the transcription of the videos, large amount of data has been collected over internet in a period matching the period of the videos. This data collected over internet has been compared to the well-known GigaWord corpora, available from LDC. The behavior (coverage) of the frequent words of each corpus, is similar between the data collected over the web and the GigaWord data. Nevertheless, the comparison shows that much more (frequent) words are needed in Arabic than in French or English to achieve a similar coverage of the word occurrences.

The collected data has been used to elaborate updated vocabularies in French, English and Arabic. Different sizes have been considered from 100 k words up to 800 k words. Noticeable reductions in the OOV rates are observed when the vocabulary size increases. The smallest OOV rates are observed on French data, and the largest ones on Arabic data.

Figure 1 :

 1 Figure 1: Display of speech recognition results achieved with the old and new vocabularies. Both speech recognition results (ASR-V01 corresponding to the baseline ASR, and ASR-V02 corresponding to ASR with the updated vocabulary) are displayed as synchronized subtitles (bottom-left) and in separate frames (bottom-right). For helping checking recognition performance, when available, the YouTube and Euronews description are also displayed (middle part).

Table 1 .

 1 Some characteristics related to linguistics aspects of the baseline ASR systems.

		French English Arabic
	Text training data (number of word occurrences)	1,620 M	155 M 1,000 M
	Vocabulary size (number of words)	97 k	150 k	95 k
	Number of pronunciation variants per word.	2.1	1.1	5.1
	Table

Table 2 .

 2 Note that during this preprocessing, all capital letters have been kept.

Table 2 .

 2 Amount of word occurrences per language for the web training data, and for the GigaWord data.

	Language	Web data	GigaWord
	French	1.9 G	0.8 G
	English	2.9 G	4.1 G
	Arabic	0.7 G	1.1 G

Table 3 .

 3 Amount of word occurrences per language in development and test sets.

		French	English	Arabic
	Development set	1500 k	1720 k	1240 k
	Test set	250 k	280 k	70 k

Table 4 .

 4 Percentage of out-of-vocabulary words in the development sets for each language and vocabulary. Sizes of baseline vocabularies are specified in Table1.

		French	English	Arabic
	Nb. words	51 k	64 k	129 k
	Nb. occurrences	1500 k	1720 k	1240 k
	Baseline (95 to 150 k)	1.8%	7.2%	17.4%
	New 100 k	0.4%	1.1%	5.5%
	New 200 k	0.1%	0.4%	3.1%
	New 400 k	0.1%	0.3%	1.5%
	New 800 k	0.1%	0.3%	0.2%

Table 5 .

 5 Percentage of out-of-vocabulary words in the test sets for each language and vocabulary. . Sizes of baseline vocabularies are specified in Table1.

		French	English	Arabic
	Nb. words	20 k	21 k	20 k
	Nb. occurrences	250 k	280 k	70 k
	Baseline (95 to 150 k)	1.8%	5.5%	16.4%
	New 100 k	0.8%	3.3%	6.8%
	New 200 k	0.4%	2.7%	4.5%
	New 400 k	0.2%	1.9%	3.1%
	New 800 k	0.2%	1.5%	2.0%

http://deustotechlife.deusto.es/amis/project/

Acknowledgements

Part of this work was supported by the Chist-Era AMIS (Access Multilingual Information opinionS) project. Also, some experiments presented in this paper have been carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https: //www.grid5000.fr).


Rubén Herzog
email: rubenherzog@ug.uchile.cl

María-José Escobar

Adrián G Palacios

Bruno Cessac

Dimensionality Reduction on Maximum Entropy Models on Spiking Networks

Introduction

The nature and importance of correlations for the neural coding, has been debated since long time, including the role of weak pairwise correlations [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF][START_REF] Shlens | The Structure of Multi-Neuron Firing Patterns in Primate Retina[END_REF]; for error supression, enhancing channel capacity [?]; or to capture the statistical properties of the network activity [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF][START_REF] Shlens | The structure of multi-neuron firing patterns in primate retina[END_REF].

The introduction of maximum entropy models (MEM) to the analysis of neural networks started to gain popularity recently. MEM have been proposed to capture and predict the collective spatio-temporal pattern activity from a vertebrate retina network [4-6, 14, 15] or from a cortex cultures cells [START_REF] Marre | Prediction of spatiotemporal patterns of neural activity from pairwise correlations[END_REF]. However, most of the studies has been limited to capture spatial correlations by using an Ising model [START_REF] Ganmor | The architecture of functional interaction networks in the retina[END_REF][START_REF] Ganmor | Sparse low-order interaction network underlies a highly correlated and learnable neural population code[END_REF][START_REF] Tkacik | The simplest maximum entropy model for collective behavior in a neural network[END_REF], where events occurring at different times are independent. Thus, these approaches have neglected spatio-temporal correlations, memory and causality, which, presumably, play a central role on the nervous system. More recently, an extension of MEM that includes spatio-temporal correlations has been proposed at the theoretical [START_REF] Cofré | Dynamics and spike trains statistics in conductance-based Integrate-and-Fire neural networks with chemical and electric synapses[END_REF][START_REF] Cofre | Exact computation of the maximum-entropy potential of spiking neural-network models[END_REF]; numerical [START_REF] Nasser | Spatio-temporal spike train analysis for large scale networks using the maximum entropy principle and Montecarlo method[END_REF][START_REF] Nasser | Parameters estimation for spatio-temporal maximum entropy distributions: application to neural spike trains[END_REF]; and experimental level [START_REF] Vasquez | Gibbs distribution analysis of temporal correlation structure on multicell spike trains from retina ganglion cells[END_REF].

In all cases, a MEM is characterized by a function ("energy" or "Hamiltonian") containing functional and/or effective interactions between spiking neurons. In its simplest form this energy contains only self-interactions and the corresponding MEM is similar to a Bernoulli model, where spikes are independent, that well fits the dominant terms in the statistics. However adding interactions, e.g. pairwise or higher, allows in principle to handle higher and higher order statistics. In its simplest form we can PLOS 1/22 assume that a neural population is coding solely by means of firing rates, while including pairwise spatio-temporal interactions we assumes the neural coding include also the temporal co-activation of the neurons [?, ?]. However, in terms of a computational model, the overwhelming explosion of parameters has corrolary effects: over-fitting and bad reliability in parameters estimations. This proliferation of parameters (i.e. increase of model dimensionality), especially whith a growing number of neurons, is a major criticism of MEM. Would be it possible to lower the dimensionality of the space where a MEM can be mapped, given the redundancies of the code? [?, ?, ?, 3] If so, we are facing a compressible MEM that reflects the compressibility of the neural code. So, in analogy with signal compression, where the presence of statistical redundancy can be exploited to map the signal onto a subset of independent channels (i.e. non redundant). In this paper, we propose to exploit the redundancies of the neural code to find an optimal set of independent MEM dimensions that captures the information about the neural code.

Here, we propose a method based on information geometry and MEM that allow us to detect the optimal set of MEM independent dimensions capturing the information contained on neural spiking data. Specifically, we propose a method that is general enough to consider spatio-temporal interactions at several lags between neurons. The core of this approach is a matrix capturing the effect of the change of one parameter of the MEM on the second-order statistics of the network activity. This matrix is known under the name "Fisher metric" in statistics and information geometry [?] and Susceptibility Matrix (χ) in statistical physics. Even when this matrix is related to the effect of varying one parameter on the network activity, it can be analytically computed without actually fitting a MEM model being then possible to perform this computation in large set of neurons (∼100 neurons). Based on the spectral properties of this matrix, we can find the optimal number of independent channels that better represent the neural activity.

The method would be first validated in synthetic data, where data would be generated with a known underlying statistics. So, using χ to propose a MEM with more parameters than the underlying data statistics, our method would helps to detect the dimensionality of the underlying statistics, even in the presence of unobserved events. Nevertheless, if the activity is too dense, the underlying statistics hides under the noisy activity, hindering the estimation of the underlying dimensionality.

Additionally, the method was tested in a neural population of retinal ganglion cells (RGC) in vitro Octodon degus recorded with 252-MEA (multi electrode array), an experimental set-up that allow to stimulate and record a population of ∼ 200-300 of RGC. These cells are part of the retina processing and sent it output trough the optic nerve to the brain. Our method was applied to three different set of neural response obtained after applying three type of visual stimuli: i) spontaneous photopic activity; ii) white-noise checkerboard and iii) a short natural movie. Based on our experimental observations, under i) we shows the activity is sparser: small correlations and silent cells. During ii) and iii) responses are denser, exhibiting highly correlated activity mainly driven by common inputs coming from overlapped receptive fields. As expected, we found that RGC activity is highly compressible (∼ 50% of the imposed model dimensionality) and that the stimuli spatio-temporal modulation increases the number of independent dimensions required to optimally represent the neural activity. This suggest that RGC population activity is able to adapts to a stimuli conditis changing the number of coding channels according to stimuli correlations.

I PLOS 2/22

Results

Context. We give here the main elements to understand the results described in this section. Technical details are given in the section "Methods". Spike trains are denoted by ω and the state of neuron i in the time bin t is denoted by ω i (t) ∈ {0, 1} (= 1 when spiking). MEMs are probability distributions assigning an exponential weight to spike trains. This weight (also called "energy") is a linear combination of specific terms called in the sequel "interactions" or "monomials".

A paradigmatic example is the Ising model where the energy at time t reads

φ(t) = i b i ω i (t) + i,j J ij ω i (t)ω j (t).
The terms b i and J ij are parameters tuning the probability, whereas the terms ω i (t) and ω i (t)ω j (t) depend on the spike configuration and are called "interactions" (self-interaction for the term ω i (t), and pairwise interactions ω i (t)ω j (t)). In Ising model, there is no coupling of spikes at different times so that the probability of observing a raster of length T factorizes into probabilities of spike states at each time: successive time events are therefore independent. More generally, we can extend energy to a general form:

φ = L l=1 h l m l . (1)
The terms h l are parameters tuning the probability. The form of the energy, i.,e. the arbitrary choice of interactions, defines a MEM with which one is attempting to "explain" data, by a suitable tuning of parameters h l . This is actually achieved by maximizing the statistical entropy under the constraints that, for each l, the average value of m l predicted by the model (noted µ [m l ]) is equal to the average observed in the experimental raster, (π(m l )). Hence the terminology Maximum Entropy Model (MEM). Considering each parameter h l as spanning a direction in a L dimensional space, a MEM is a point in this space which can be very huge. Indeed, the number of parameters increases with the number of neurons leading rapidly to a plethora.

Our method proposes an analytic way to find the optimal set of linearly independent dimensions that represents more accurately the inner structure of the neural code. It is based on measuring the interdependence between the MEM parameters and re-mapping those parameters on a different coordinate systems where each dimension is a linear combinations of them and different dimensions are linearly independent. Once we are

PLOS

3/22 on this new coordinate system, we can find an optimal set of dimensions based on the geometric properties of the MEM (see eq. 13). Thus, finding the optimal number of dimensions is equivalent to finding the optimal number of independent coding channels of the neural code under this MEM assumptions: coding is based on firing rates, spatial and temporal correlations.

We have applied the method to two types of data: synthetic and retinal, where the synthetic was defined using different MEM:

1. Independent model: Neurons are independent, i.e. there are only self-interactions (firing rates) of a neuron to itself and the parameter associated to this interaction controls the firing rate of the neuron. A rate model with N neurons has therefore L = N parameters. We call this model Indep., because neurons are independent. and 270 neurons for each recording. We also generated shuffled versions of these rasters, where the neurons firing rates were maintained, but the dependence between them were disrupted by randomizing the individual spike trains.

General Method Description and Applications to Synthetic Data

Our analysis relies on a matrix χ, called Susceptibility matrix in the sequel, using here the statistical physics terminology. Its detailed definition and computation is presented in the Methods. Here it is sufficient to know that it is numerically computed from the empirical time correlations between monomials (see eq. 12) without needing to fit a MEM on the observable data. Susceptibility matrix has the following properties: (ii) χ is symmetric and positive, as a consequence, it has real positive eigenvalues1 (iv) From the spectrum of χ one can define an ellipsoid with principal axes the eigenvectors v k and with extension 1

λ 1 ≥ λ 2 ≥ . . . λ k ≥ λ L > 0.
√ λ k
on the principal axis k (Fig 1C). Clearly, the smaller λ k the longer the ellipsoid in the direction k. As elaborated in the methods section, this ellipsoid delimitates a volume, called confidence volume.

Points inside the confidence volume correspond to indistinguishable models reproducing equally well the empirical observations (Fig 1C) within an accuracy (see eq. 13 for more details). As a consequence, the inverse of the eigenvalues tells us how much a small variation in the estimation of the parameters affects the statistics. For a large eigenvalue λ k , a tiny variation in the direction v k has a dramatic impact on statistics. On the opposite, small eigenvalues correspond to sloppy dimension where a big change in the corresponding direction has small impact. The notion of stiff and sloppy dimension in statistics is not new and has been used by several authors, including for the analysis of spike trains [?], but the treatment we propose, for MEM with spatio-temporal interactions (in contrast to previous papers dealing with Ising model) is, to our best knowledge, a novelty. k c : a parameter representing the structure of the spike train interactions Independent on the nature of the synthetic data, we computed the χ matrix for a PWR2 statistics, i.e., overfitting the data generated from an independent model. From the χ matrix we obtained their eigenvalues spectrum, which is shown in Fig 1B .     The eigenvalue spectrum presents a particular shape which analysis is the core of our method. Moving along the spectrum from left to right (increasing index k, decreasing eigenvalue λ k ) we observe a first sharp decrease (cut-off ) at k = N , for the Indep. raster (Fig 1B ,black) and a corresponding minima on the volume (Fig ??A, black square). We also observe a second cut-off further in the spectrum. The k first eigenvalues and associated eigenvectors define a k-dimensional ellipsoid with volume V(k, ) depending on k and (eq. [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF]. The relation between the volume of the ellipsoid and the number of the MEM dimensions given by k, for a given accuracy , is shown in Fig 1D . The presence of two cut-offs means that there are two values of k where the confidence volume is minimal (Fig 1D ),i.e. the model is more accurately determined. The first cut-off occurs at k = N and it is a trivial solution representing only the statistics given by the neuron firing rates and not capturing the spatio-temporal interactions. The second cut-off, from now on called k c represents a second minimal point of interest, capturing the spatio-temporal interactions, and which value is proposed in this article.

According to our results, the volume of indistinguishable models depends on the number of dimensions considered in the χ matrix (k) and the imposed accuracy (Eq. 13). A large volume indicates sloppy directions that could be neglected, while a small volume is a straight-forward dimensionality reduction framework. To illustrate this idea Fig 1C shows a system that lives in a 3D space represented as an ellipsoid, whose principal axes (dimensions) are linear combinations of the original axes (x, y, z ) and the

PLOS 5/22
amplitude of these dimensions is given by the inverse of the corresponding eigenvalues, representing the degrees of freedom that the system has on each dimension. k c depending on the the network size and recording length As we just described, k c is a global parameter measuring the number of dimensions at a given accuracy that minimizes the confidence volume, representing the minimal set of models that are equally good explaining the observed data. Each of this dimensions, i.e. eigenvectors, are a linear combination of model parameters, as expected from the eigendecomposition process of a symmetric and positive matrix (recall that the dimensionality of the matrix is given by the number of model parameters). So, if the underlying statistics is random (no linear dependences between model parameters) and has the same dimensionality than the imposed model, k c approximates L for an infinite raster (Fig S2A -B). Otherwise, if there is too much noise or if there are linear relationships on the underlying statistics, we expect to see k c < L. Formally, we can decompose the full MEM dimensionality as However, increasing the network size, as stated before, increases the number of possible patterns rapidly, requiring longer and longer rasters in order to observe all the possible patterns. So, for a finite raster, k c and L will diverge as the network size grows, given that some patterns will never occur in that time. Nevertheless, despite the unavoidable finiteness of the rasters, at N ≤ 80 and T = 10 k c values on the independent and spatio-temporal correlated cases Using the MEM described above (Indep. and PWR2), we generated 100 different rasters (N = 20, T = 10 6 ) for each MEM, looked for the minimal volume and corresponding k c (Fig ? ?A). For Indep. we see two minima, where the global one is the first, showing that we can capture the second-order statistics in a raster without this kind of underlying statistics, but would not correspond to the global minima (i.e. overfitting). In the PWR2 case we see just one minima, beyond which the volume explodes. Thus, we found for Indep. case k c = 19 ± 0 and for PWR2 k c = 447 ± 7.00 (mean ± s.d.) (Fig ??B). In the Indep. case we found that k c corresponds almost to the full dimensionality of the underlying model (L = 20), while in PWR2 k c corresponds to approximately 75% of the full underlying model dimensionality (Fig ? ?C). In the first case the divergence between k c and L (just one dimension) is explained by the presence of one low firing rate neuron (data not shown), which according to our framework, is a neglectable dimension. In the second case the divergence between k c and L is given mainly by the unobserved events, which will be shown in more detail in the next section.

L = k c + κ + N.O.,

PLOS

6/22 k c as a measure of code compression

As we previously remarked, the network size does not affect the shape of the spectrum for random PWR2 rasters (Fig S2C), but just its length, i.e. the number of eigenvalues (trivially explained by the increased dimensionality of the MEM). This shape invariant property comes from the fact that the underlying statistics lacks of linear dependencies by construction, so there should be approximately one dimension per parameter. Thus, finding k c L comes mainly from two facts: (i) unobserved events and (ii) linear dependencies between model parameters. The first case can be easily overcame by adding to k c the number of unobserved events, making rasters of different recording length comparable. On the second case we are facing compression, where given the linear dependencies between model parameters many of them are mapped onto the same dimensions, i.e. the model can be described by less dimensions than the number of parameters. So, for an experimental recording with no unobserved events, the value of k c shows how much all the variables of the MEM (i.e. our assumptions about the neural code) can be compressed on a lower dimensional space. The eigenvalue spectrum of the χ matrix, for all the cases obtained through f ac , were computed (Fig ? ?C). Increasing f ac has 3 main effects on the spectrum: (i) the first eigenvalue is decreased, (ii) the offset spectrum is decreased, and (iii) the number of eigenvalues above the minimal observed probability is reduced (1/T, the lower limit of the ordinate on the plot). The first 2 effects come from χ computation based on monomials correlations at different time lags, so the higher the probability of the monomials (density of the raster), the denser χ is and the higher the entries of it, yielding as a consequence bigger eigenvalues. The third fact comes from the unobserved events, that will be related to vanishing eigenvalues.

Scaling the Underlying Statistics

Differences on the spectra reflect differences on the raster second order statistics: (i) dense rasters are dominated by high-firing rates where almost any kind of interaction is highly probable, (ii) an intermediate situation with moderate firing rates where some interactions are highly probable and some others are not, and (iii) a very sparse activity with low firing rates where many pairwise interactions are unlikely to be observed.

Regarding k c , we see that none of the modified rasters yields a k c comparable to the original PWR2 (Fig ? ?D). In the case of f ac < 1 there are no unobserved events (Fig

??D)
, so the reduction of dimensionality should come from a different fact. The high density of the raster makes almost any kind of pattern highly probably, hiding the defined underlying statistics under all the other (not defined) observed patterns,

PLOS 7/22
resulting in a k c = N , like an independent case. Increasing f ac , but keeping it below 1, increases k c given that the activity become sparser, but still there are many pairwise interactions that are hiding under the increased raster density. In the case of f ac > 1, there are many unobserved events, consequence of the very negative underlying rates parameters, yielding smaller k c values, compared to the original PWR2. However, if we add the number of unobserved events to k c , PWR2 and the f ac > 1 versions become comparable, showing that the dimensionality reduction observed for those raster comes mainly from the unobserved events.

Then, dealing with random underlying statistics the dimensionality reduction will depend on the density of the raster: on one side, very dense rasters are dominated by the firing rates, hiding the underlying second-order statistics on the noisy activity.

Under this situation we are not able to recover the original underlying statistics dimensionality. On the other side, very sparse rasters will show many unobserved events, reducing the observed dimensionality. However, this situation can be easily overcame adding the number of unobserved events to k c , recovering the original underlying statistics dimensionality. Green square shows a slice of this raster, which is considered as 1 (or 0) if it generated (or not) activity inside a time window of width equal to the bin size. For this raster 2 types of pairwise interactions are defined: spatial interactions (blue) and temporal interactions with R = 2, i.e. one time-step between spikes (red). Spatial interaction is exemplified as w 3 (0)w 4 (0) (both neurons firing at the same time) and w 3 (1)w 4 (0) (neuron 3 firing one bin after neuron 4). B: Susceptibility matrix eigenvalue spectrum in log-log scale for Independent (black) and Pairwise Range=2 (PWR2, red). Black vertical dashed line denotes the network size (N = 20). Indep. raster shows a sharp cut-off close to the network size, while PWR2 shows a monotonic decay of the eigenvalues magnitude, without a clear cut-off. C: Illustrates a 3D system representing a volume of indistinguishable models living in x, y, z dimensions which is rotated into a new coordinate system x', y', z', that are linear combinations of the original coordinate system (original coordinate system denoted as dashed lines). Bottom shows the 2D projections, where the amplitude of the ellipsoid on each dimension is the inverse of the corresponding eigenvalue. In this case x' and y' are the stiff directions, while z' is a sloppy direction. D: Shows the value of the volume of indistinguishable models, V(k, ) as the degrees of freedom (k) increases. The volume reaches a minimum (orange dot) for certain k named k c . Increasing k upon k c makes the volume explode. Increasing/decreasing one order of magnitude the imposed accuracy ( ) yields no convex functions, i.e. no minimization.

Dimensionality reduction on retina data

We also tested the dimensionality reduction method (finding k c ) on in vitro retina data recorded with multi-electrode arrays under 3 stimuli conditions: photopic spontaneous activity (PSA), white noise (WN) and natural movie (NM), sorted in increasing level of stimuli high-order correlations. The stimuli high-order correlations increases the probabilities both of the firing rates and the pairwise interactions, i.e. increases the raster density. However, the distribution of the pairwise interactions probabilities of the empirical recordings and their shuffled version for all stimuli are not significantly different (bin 10ms, Mann-Whitney test P > 0.05, Fig S3), showing that the distributions of the pairwise interactions probabilities are given by the independent neuron firing rates. However, we are not interested in the magnitude of the pairwise A: shows the underlying parameters distribution, split by firing rates (blue) and pairwise interactions (red). The bigger the scaling factor fac, the more negative the rates parameters and the more wide the interactions parameters distribution. B: shows the corresponding monomials probabilities for the scaled rasters, split as firing rates and interactions, as A.

PLOS 8/22
We see that increasing fac has the effect of decreasing both rates and pairwise interactions probabilities, reaching the point where many pairwise interactions vanish (f ac > 1). C: are the corresponding χ eigenvalue spectra (average out of 10 rasters) for the scaled rasters. Increasing f ac decreases the first eigenvalue, the spectrum offset and the number of eigenvalues above the minimal observed probability (1/T). D: k c values for the scaled rasters. None of the scaled rasters shows k c value as the one obtained for the original PWR2 raster. Dost are averages and errorbars 1 s.d. of 10 different rasters with the same underlying parameters. E: k c values plus the number of unobserved events. For f ac > 1 adding the unobserved events yields values close to the original PWR2 raster, showing that the dimensionality reduction obtained for those rasters is given mainly by the unobserved effects. For f ac < 1 we see no unobserved events, so the dimensionality reduction obtained for those cases is given mainly by the increased density of the raster hiding the underlying statistics.

interactions, but instead, we are interested on the inner structure of the neural code, its dimensionality and compressibility, which is captured by χ and k c .

To this end, we computed the χ eigenvalue spectrum of 30 random sub-networks of N = 50 from the total number of neurons recorded in each of the four experiments, under the 3 stimuli conditions. According to our synthetic rasters experiments, for N = 50 and T ∼ 10 6 we can get good k c estimates (Fig S2B andS2B), which fits with our experimental recordings. We applied the same procedure to the shuffled version of the rasters, to compare with rasters having exactly the same firing rates distribution, but with no dependency between neurons.

Experimental versus shuffled spectrum

Similar to what we obtained for scaled synthetic rasters (Fig ? ?C), we see that the first eigenvalue, the offset and the vanishing eigenvalues (below 1/T) of χ spectrum increases with the raster density, which in this case is driven by the stimuli (Fig ? ?A). In addition, we only see a clear cut-off close to N for WN, while the other two conditions the cut-off is not clear. Counter wise, all the shuffled rasters show a cut-off close to N , suggesting that in the experimental recordings there are significant linear dependences

PLOS 9/22
between monomials that are not present in the shuffled version. Specially, NM shows a smooth decay of the spectrum, similar to the observed for PWR2, which is highly modified when the raster is shuffled, having a sharp cut-off close to N . The k c analysis on retinal data revealed that RGC activity is not random, showing almost the half of dimensionality compared to a random PWR2 with the same network size. Furthermore, k c corrected by the number of unobserved events (Fig

??D)
decreases with the bin size as a consequence of a larger window in which spikes can be correlated, increasing the interdependence of RGC population activity (both for empirical and shuffled data). Specifically, we see that k c increases with the stimuli high-order correlation for fast time scales (Mann-Whitney test, P < 10 -5 for all comparisons). This suggests that for fast time scales the retina increases the number of coding channels as the stimuli high-order correlation increases. However, for medium time scales the pictures changes, showing no significant differences between WN and NM for 10ms (Mann-Whitney test, P > 0.1) and showing higher k c values for WN than NM for 20ms (Mann-Whitney test, P < 10 -5 ). This suggest that at larger time scales the RGC activity under NM becomes more interdependent than for WN.

On the other hand, the shuffled rasters show higher k c values for all the bin sizes and conditions (Mann-Whitney test, P < 10 -4 ), except for 1ms where is significantly small (Mann-Whitney test, P < 10 -4 ) and for PSA at 5ms, where they are not significantly different (Mann-Whitney test, P > 0.1). However, when corrected by the number of unobserved events (Fig ??D), the picture is the same for all bin sizes: shuffled rasters have always higher k c values than the empirical ones (Mann-Whitney test, P < 10 -6 and P < 0.01 for WN at 1ms). This confirms that the RGC neural code has interdependences that can be mapped onto a lower dimensional space (i.e. compressed), compared to the shuffled version, that lacks of interdependences by construction. Yet, the k c obtained for the shuffled raster is still very small (almost a half) compared to the random PWR2 raster, which suggests that just the firing rates distribution introduces some non-random interdependences in the neural code, allowing compression.

In sum, RGC neural code is highly compressible compared to a random raster of same network size. Even compared to rasters with exactly the same firing rates distribution, that also reproduces the pairwise interactions probability distribution (Fig S3), the RGC neural code is more compressible. Nevertheless, a significant compression can be achieved considering just the firing rates distribution. Besides, stimuli high-order correlations increases the raster density and the number of independent coding channels increases with the stimuli spatio-temporal modulation. For fast time scales, k c increases with the stimuli high-order correlations. Finally, at medium time scales the code for NM becomes more compressible than WN, suggesting higher redundancy on the code under NM.

Fig 4. Dimensionality reduction on RGC data.

A: shows the average eigenvalue spectrum (solid line) of RGC data (N = 50) under 3 different stimuli conditions, with 10ms bin size. Stimuli high-.order statistics increases the first eigenvalue, the spectra offset and the number of eigenvalues above 1/T. Except from WN, there is no clear cut-off close to N . Shaded area is 1 s.d. out of 30 sub networks. Black line is a random PWR2 rasters of the same network size than the RGC raster. B: shows the same than A, but for the shuffled version of the empirical rasters. We see that all of them show a clear cut-off close to N and that they preserve effects induced by stimuli high-order correlations, i.e. change on the first eigenvalue, offset and number of eigenvalues above 1/T. The differences on the cut-off suggests that experimental data has an inner structures that is not merely explained by firing rates. C: show box plots for k c values for empirical and shuffled rasters. Except from 1ms, k c is higher for shuffled data in all bin sizes, under all conditions. Stimuli spatio-temporal modulation significantly increases k c . Also, neither of both type of rasters reaches the k c values obtained for a PWR2 rasters, which shows almost no dimensionality reduction. Note that k c decreases as the bin size increases, showing that for bigger bin sizes the inner structure of the code becomes more interdependent. For fast time scales (1 and 5ms), k c increases with the stimuli-high order correlations, but at 10ms WN and Nm are not significantly different and for 20 ms WN is larger than NM, suggesting that bigger bin sizes capture more redundancies on the code. D: shows box plots k c values corrected by the number of unobserved events. Now for all bin sizes the shuffled data has bigger values than the empirical one, showing that the RGC data has an inner structure with interdependences between variables that is not present on a shuffled raster. However, shuffled data has values almost a half than a PWR2 raster of the same size, suggesting that the firing rates distribution only is enough to perform compression on the neural code. The effect of stimuli spatio-temporal modulation and stimuli high-order statistics remain.

Discussion

In this paper we have proposed a method to reduce the dimensionality of MEM on artificial and biological spiking networks. It is grounded on information geometry of the matrix χ, which characterizes how a small variation of parameters impacts the statistical estimations. The χ matrix captures the interdependences between the neural code variables. After an eigendecomposition process, the eigenvalue spectrum of χ exhibit two cut-offs. The first one shows that, both in synthetic as well as in retina data, large part of statistics is "explained" by the neurons firing rates. On the opposite, the second cut-off (here called k c ) reflects a non trivial effect associated with higher order statistics. As the eigendirections on the right part of the cut-off correspond to noise, the spectrum lying between the two cut-off contains a relevant information associated to statistics of second order.

The reduction of the MEM dimensionality is directly linked with data compression, obtained from the linear dependencies between variables of the MEM, i.e., higher order interactions between neurons. For example our analysis in the case of synthetic rasters, where both the firing rates and the pairwise interactions are defined randomly, show a k c value very close to L (maximal dimensionality). This demonstrates that if there are no linear relationships between the parameters by construction, the code is not compressible and we have almost one dimension per parameter. Counter wise, retina data show a significant compression of ∼ 50%, as expected from a neural tissue where cells are driven by common inputs and cells are electrically coupled [?], increasing the level of dependency between them. This compression reduces the dimensionality of the MEM to a lower dimensional space, where each dimension is a linear combination of model parameters, characterizing the population activity by a set of independent

PLOS 11/22
dimensions representing the inner structure of the network activity.

Limits of the method

[MJ: BRUNO: from where this method comes from? or this manner to compute the susceptibility matrix is a result of this work?] The first limitation of our method comes from the numerical approximation used to compute χ matrix, which is obtained by summing monomials correlations at different time lags (Eq. 12). In general, this function decays exponentially with time, but adding more time lags on this sum adds more noise to the matrix, reaching a point where the matrix is highly dominated by noise. To truncate this approximation we need to consider a trade-off between temporal resolution and reduction of noise. To this end we use 4 time lags (i.e. 4 terms of the sum), which is equivalent to the double memory depth used in the model (R = 2). This numerical estimation of χ also imposes limits on the method, given that considering the infinite memory case (R → ∞) will generate a χ matrix that is governed by noise.

On the other hand, it is possible to compute χ from the model parameters, requiring a previous model fitting step, as done by [?]. However, for N > 20 and R > 1 this computation becomes more and more prohibitive as the network size and the memory depth of the model increases. So, despite the numerical approximation and its intrinsic errors, computing χ from the empirical monomials time correlations is the best approach we found to work with medium size networks and spatio-temporal constrains.

Furthermore, computing χ this way, instead of fitting the MEM in advance, give us geometrical information about the model and, ultimately, about the linear dependences of the neural code. So, it is possible to have insights about how a neural population modulates its activity given the stimuli, increasing or decreasing the number of independent coding channels, without a fiting procedure. Not less important, fitting the model will give information about the sign (positive or negative) and the magnitude of the interactions (weak or strong), giving detailed information about the network topology. Nevertheless, the scope of this work was not to fit different models on data and test its performance (e.g. Bayesian or Akaike information criterion that takes into account the model parameters and likelihood [?]) neither study changes in network topology under different stimuli. Instead, we focus on exploring the geometrical properties of the MEM and its meaning in terms of the neural code redundancy and compressibility.

k c is not one value, but a set of values

The main challenge this methodology introduces is the selection of k c . The selection of k c is related to the minimization of the so-called confidence volume, which not only depends on the number of dimensions given by k c , but also on the imposed accuracy . As shown on Fig S1, the confidence volume (logV(k c , )) decays monotonically as k c increases and decreases, reaching an inflection point where even if we decrease the accuracy by orders of magnitude k c remains almost unchanged, so adding more dimensions only reduces the accuracy. It could be also possible to choose the first cut-off as k c (or any other k), based on the difference between consecutive eigenvalues (another criteria to find the cut-off on eigenvalue spectra), but from our observations on synthetic data we know that the first cut-off is related mainly to the firing rates, so it misses the raster high-order statistics. We also know that even when the underlying statistics defines pairwise interactions but the firing rates are too high, the density of the raster hides the higher order statistics, yielding k c = N , considering the dimensions related to higher order moments as noise.

Even though we found two minima on the volume function, the global minima for experimental and shuffled data was always at the second minima. The only case where

PLOS 12/22
the first minima was the global minima was for the Indep. synthetic raster (Fig ? ?A, black trace), which is in agreement with the underlying statistics.

Comparison with similar analysis

To our knowledge, there is no previous work related to the analysis of the χ matrix considering spatio-temporal pairwise interactions applied to neuronal networks. The authors in [?] proposed a similar analysis considering only spatial interactions, i.e. the Fisher Information Matrix (FIM) for the Ising model on in vitro, in vivo and in silico networks. The work of [?] studies small neural networks (10 neurons), they looked for stiff neurons, which are related to stiff dimensions (the largest FIM eigenvalues),

proposing that those neurons are the ones giving stability to the network, while the neurons involved on the sloppy dimensions (dimensions where the parameters can have significant changes without affecting the model) are the ones involved on plasticity, allowing the network to remodel its connections. On our approach, which involves larger networks (N = 50), we deal with neurons and their spatio-temporal interactions, exploiting the linear dependences between them to find a minimal set of dimensions that better represents the neural code. To this end, we found analytically two set of stiff dimensions: the ones before the first cut-off, related mainly to neuron firing rates and the second set, after the first cut-off, related to spatio-temporal interactions. According to our framework, the sloppy dimension would be the ones beyond k c , but we could also interpret the first N dimensions as the stiff dimensions, the ones between N and k c as the sloppy dimensions and the ones beyond k c just noisy dimensions. This re-interpretation arises from the large magnitude difference between the first and second set of dimensions (at least one order of magnitude).

In our case, the analysis was extended to spatio-temporal interactions in larger neural networks, and more importantly, we used χ as a tool for finding linear dependences between the neural code variables, compressing the observed data and characterizing the neural activity on the basis of independent dimensions.

Comparing the inner structure of the relevant dimensions between different stimuli conditions could give us insights about which channels are invariant between stimuli and which ones are adaptive/plastic. However, this is not the scope of this article and the analysis of invariant or adaptive channels will be kept for a future work.

Recently, Battaglia et al [?], studying large scale networks between brain areas, proposed a concept called Meta Connectivity, which instead of analysing the correlation between nodes of a network (the usual functional/effective connectivity analysis), they focus on the correlations between the network interaction along time. This means focusing not on the coupling w i (0)w j (0) between neurons w i and w j , but focusing on the interactions between the couplings w i (0)w j (0) and w k (0)w l (0). This provides information about high-order correlation for at least 3 nodes of the network (e.g. case of i = k) and captures the relationships between modules of network activity. Thus, χ matrix is both a functional/effective connectivity matrix (the matrix entries related to the correlations between firing rates) and also a meta connectivity matrix (the matrix entries related to correlations between pairwise interactions), providing information about the spatio-temporal interaction between network nodes and network modules.

The extension of our analysis towards the understanding of the meta connectivity has never been applied to networks of neurons and we proposed it as a future research direction, where the focus is on the variability and dependence of the interactions.

Stimuli-induced changes on RGC population activity

Retina data has significant high-order correlations, including pairwise spatial [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF][START_REF] Shlens | The Structure of Multi-Neuron Firing Patterns in Primate Retina[END_REF], temporal [START_REF] Vasquez | Gibbs distribution analysis of temporal correlation structure on multicell spike trains from retina ganglion cells[END_REF] interactions, triplets [START_REF] Ganmor | Sparse low-order interaction network underlies a highly correlated and learnable neural population code[END_REF] and groups of neurons [START_REF] Tkacik | The simplest maximum entropy model for collective behavior in a neural network[END_REF]. This correlations have

PLOS 13/22
been widely studied under MEM, which can accurately reproduce the raster spatio-temporal patterns. Nevertheless, there has been no work devoted to reduce the MEM dimensionality considering the inner dependences between the population activity variables.

Here, as a proof of concept, we used retina data of a diurnal rodent under 3 different stimuli with different statistics, from photopic spontaneous activity (spatio-temporal uniform full field, no second order statistics), a spatio-temporal white noise (gaussian statistics) to a repeated natural movie (high-order correlations both on time and space). From our analysis, we know that this stimuli high-order correlations increases the magnitude of both the firing rates and the raster high-order correlations, even making silent cells to fire. This could be related to recruitment of specific cell types by the stimuli features (e.g. local contrast [?], optic flow, color [?], among others).

In general, most of the correlated activity observed in the retina could be attributed to the receptive field overlap between recorded RGCs and to shared common noise [?]. Additionally, electrical coupling are highly present in the O. degus retina [?] inducing fast correlations between close cells. These three effects modifies the pairwise spatio-temporal interactions observed in real data. We explored the presence of these correlated activity in the neural code. To do this we compared the results of recorded data with a suffled version of it, which preserves the firing rates distribution.

Interestingly, both data sets share the pairwise interactions probability distribution suggesting that statistics up to the second order are preserved. Nevertheless, our method exhibits significant differences in the compression capability of each data set.

RGC presents a maximal compression compared to the shuffled data, suggesting that statistics of higher order are needed to fully characterize the response of RGCs.

Compressibility of the RGC code

Then, in order to study the compressibility of the RGC code, we studied χ spectrum and k c for RGC under three stimuli conditions, finding that RGC population code adapts to stimuli conditions by changing the number of independent channels.

The first difference we found between stimuli was on χ spectrum, which in correspondence with the stimulus-induced increase of the monomials probability (Fig S3), shows an increase on the offset, i.e. the eigenvalues increases their magnitude as the stimuli high-order correlation increases. But given the way χ is computed (see Eq. 11), the shape of the spectrum comes not only from the increased monomials probabilities, but also from the dependence between firing rates and spatio-temporal interactions and the dependence between spatial and temporal interactions, all of them captured by the matrix. On the other hand, shuffled data, where there are no dependences between the neural code variables, exhibits these differences on the eigenvalue magnitudes (Fig ? ?B), but changing the first cut-off, showing that eigenvalues magnitude are closely related to the raster density (also shown for synthetic data on fig ??C), while the cut-off is related with the linear dependences between the monomials.

The second difference we found between stimuli is on k c , i.e. our approximation to the dimensionality of the neural code. As expected from the stimuli statistics, k c (and also its corrected version by the unobserved events) is always lower for PSA than for the dynamic stimuli, suggesting that the network optimizes the number of dimensions required for coding the stimuli depending on the stimuli statistics. In terms of metabolic cost, a very redundant stimulus as PSA (which has the same spatial and temporal information all over the stimuli space), should be coded with less dimensions than a stimuli with more independent components (less redundant), thus, optimizing the metabolic resources.

However, for bin sizes of 1 and 5ms we observe that k c is higher for NM than for of 10ms WN has the same number of relevant dimensions than NM, while for 20ms WN has more dimensions than NM. So, for fast time scales, we face a non-optimal situation, because NM has more redundancies than the WN. It could be possible that at this time scales we are not capturing the inter-dependences of the neural code that are relevant for the brain, given that at 20 ms bin size we see the expected effect: the system exploits the stimuli redundancies and exhibit more compression for NM than for WN.

Coincidently, many MEM on retina have been done using 20 ms as bin size [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF][START_REF] Shlens | The Structure of Multi-Neuron Firing Patterns in Primate Retina[END_REF], which in our case is the bin that allow the highest compression. In addition, synthetic data shows that if the raster is too dense the underlying statistics hides under the noisy activity, which could also be possible at large values of bin sizes. To control this situation we used shuffled rasters, which preserves the same raster density. In the shuffled rasters we see that k c increases with the raster density, discarding the effect of density on the changes of dimensionality at higher bin sizes. Thus, at large bin sizes the interdependences of the neural code are responsible for the compression effect and not the raster density hiding some events. Nevertheless, we do not know in advance what time scale(s) is(are) actually relevant for the brain and neither if our assumptions about the neural code (firing rate and spatio-temporal interactions) are right, so the choice of the bin size and code variables is still an open question and somewhat arbitrary.

Finally, our work is related to the idea that a stimuli-dependent network of in silico noisy spiking neurons adapts its code according to noise and stimuli correlations [?],

instead of using just one way of coding. On our case, the stimili-dependent network is a biological one, so we don't have access to modify the noise of each neurons nor the network noise. Instead, we can just modify the stimuli correlations, which changes the dimensionality of the code. This change in dimensionality could reflect the smooth interpolation between encoding strategies: highly redundant stimuli evokes less code dimensions than stimuli which presents high-order correlations, suggesting that the MEM dimensionality is also a measure of the code redundancies. However, the analytic relationship between dimensionality reduction and the coding strategies requires an extensive mathematical and computational research that is not developed here, but we give an indirect way of studying the interdependences of the neural code as the stimuli conditions changes.

Materials and Methods

Model Definition

A spike train, ω, is a discrete-time sequence of events that represent the activity of a neural network of size N over a time period T time-steps. The spiking activity of neuron i at discrete time t is mathematically defined by a variable:

ω i (t) = 1, if neuron i fires at time t, 0, otherwise. (2)
A spike pattern {ω(t)} i=1...N , represents the activity of the whole network at a given time t. Finally, a spike block, ω t2 t1 = {ω(t) t1≤t≤t2 }, represents the activity of the whole network between two instants t 1 and t 2 . The range of a spike block is the number of time steps t 2 -t 1 + 1; the degree is its number of spike-events.

A monomial is a function which associates to a raster a real value: , instead of the whole raster. R is called the range of the monomial. From now on we note D = R -1.

m l (ω) = m k=1 ω i k (t k ), 0 ≤ i k < N, 0 ≤ t k < R. (3)
We note π (T ) ω the empirical probability constructed from the raster ω, where the average is performed by time average over a time interval of length T . This way of averaging is relevant only if one assumes that the underlying process that generated the raster is stationary. Thus, the empirical average of a monomial m l is noted π

(T ) ω [m l ].
We want to construct a stationary Markov chain (i.e. where transition probabilities are time-translation invariant) whose invariant distribution µ "approaches at best" the empirical data. Approaching at best means here minimizing Kullback-Leibler divergence

d KL π (T ) ω
| µ between the empirical probability of the spike train ω of size, π

ω , and the invariant probability distribution of the estimated Markov chain, µ. For this, we use the Maximum Entropy Principle (MEM).

As opposite to most existing MEM, like Ising model [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF][START_REF] Shlens | The Structure of Multi-Neuron Firing Patterns in Primate Retina[END_REF][START_REF] Tkačik | Ising models for networks of real neurons[END_REF][START_REF] Tkačik | Spin glass models for a network of real neurons[END_REF][START_REF] Tkacik | Optimal population coding by noisy spiking neurons[END_REF] or extensions with triplets and quadruplets interactions [START_REF] Ganmor | Sparse low-order interaction network underlies a highly correlated and learnable neural population code[END_REF], we consider here models having a memory so as to treat time correlations. The main developments of the method and application to retina data have been published in [START_REF] Nasser | Spatio-temporal spike train analysis for large scale networks using the maximum entropy principle and Montecarlo method[END_REF][START_REF] Nasser | Parameters estimation for spatio-temporal maximum entropy distributions: application to neural spike trains[END_REF][START_REF] Vasquez | Gibbs distribution analysis of temporal correlation structure on multicell spike trains from retina ganglion cells[END_REF][START_REF] Cessac | Spike train statistics from empirical facts to theory: the case of the retina[END_REF]. The Maximum Entropy Principle corresponds to finding a stationary probability distribution µ which maximizes the statistical entropy S(µ) given the constraints:

µ [m l ] = π (T ) ω [m l ] , l ∈ S. (4)
This means the average µ [m l ] of monomial m l predicted by the model, equals the raster empirical average π

(T ) ω [m l ]
, where m l belongs to a prescribed list of monomials S, where S defines the model-type. If S contains only monomials of the form ω i (0) the model is Bernoulli. If S contains monomials of the form ω i (0) and ω i (0)ω j (0), i = j the model is Ising. If S contains monomials of the form ω i (0); ω i (0)ω j (r), i = j, r = 0 . . . R the model is pairwise with range R. Here, a model with range R = 2 is used.

It can be shown that there is a unique probability distribution µ S , depending on the set of constraints S, maximizing the statistical entropy under the constraints (4). µ S has the following properties:

φ = l∈S h l m l . (5)
is called "potential" such that µ S is a Gibbs distribution for φ. Being a Gibbs distribution means that the probability of a spike block ω n 0 obeys:

µ S ω n-1 0 ∼ e -(n-D)P e φω n-1 0 (6)
where:

P = lim n→∞ 1 n log Z n , (7)
is called free energy density and

Z n = ω n-1 0 e φω n-1 0 . ( 8
)
is called n time steps partition function

PLOS 16/22
Susceptibility Matrix, χ P is a function of h l s. One can show that the average of m l predicted by the Gibbs distribution µ S is given by the derivative of the free energy with respect to h l [START_REF] Keller | Equilibrium States in Ergodic Theory[END_REF] :

µ S [m l ] = ∂P ∂h l . (9)
Therefore, the second derivative:

χ ll = ∂ 2 P ∂h l ∂h l , (10)
tells us how much the average of the monomial m l is modified when we slightly modify the coefficient of monomial m l . We call the matrix χ with entries χ ll the susceptibility matrix.

χ is numerically computed without the need of fit the model. This is achieved by computing the correlations between monomial m l at time 0 and monomial l at time n is:

C l,l (n ) = µ[m l (0) m l (n)] -µ[m l ]µ[m l ]. (11)
In the case where m l = ω i (0), m l = ω j (0), this gives the time pairwise correlations between neurons i,j at time n. C l,l (n ) depends on the model S.

One can show that:

χ ll = C l,l (0 ) + 2 +∞ n=1 C l,l (n ) . (12)
Thus, χ ll integrates all time correlations between monomial m l and m l . χ is symmetric positive thus it has real positive eigenvalues, but given its numerical computation, as more terms are considered on the sum of (12) more noise is captured by χ, yielding several negative eigenvalues (data not shown). For that reason, n was constrained to 4, being the double of the range of the model. When the potential φ only contains range 1 monomials successive time steps are uncorrelated and (12) reduces to χ ll = C l,l (0 ) the so-called fluctuation-dissipation theorem in physics.

Volume of indistinguishable models

Functions of the form (5), parametrized by the coefficients h l , live in a L dimensional space E. Since each set of h l defines a maximum entropy (Gibbs) probability distribution, each point of E represents such a probability. This space is very huge, and a generic probability in this space have all h l s different from 0. Depending on the dynamics and specific properties of the studied system, the h l s are typically distributed onto a sub-manifold W in E, corresponding to functional relations between the h l s parameters, [START_REF] Cofre | Exact computation of the maximum-entropy potential of spiking neural-network models[END_REF]. The right object to study both the geometric structure (metric) of the manifold W and the second order statistical fluctuations associated with the central limit theorem is the Fisher matrix or Susceptibility matrix [START_REF] Nakahara | Information-Geometric Decomposition in Spike Analysis[END_REF][START_REF] Mastromatteo | On the typical properties of inverse problems in statistical mechanics[END_REF][START_REF] Sessak | Inverse problems in spin models[END_REF][START_REF] Sessak | Small-correlation expansions for the inverse Ising problem[END_REF]. Given χ properties, it has eigenvalues λ k , k = 1 . . . L ordered by decreasing value (i.e. the spectrum) without loss of generality and corresponding eigenvector v k associated with λ k . The v k s define an orthogonal eigenbasis composed of a linear combinations of model parameters. From the statistical perspective, 1

√ λ k
gives the amplitude of second-order fluctuations in the estimation of coefficients h l projected on direction v k . The smaller the eigenvalue, the larger the amplitude. When projecting back to the original basis, these errors propagate to h l s estimation. From the geometrical perspective, as χ is a metric, i.e. the inverse of a square distance, 1 √ λ k is a distance on the manifold W. This distance can be interpreted as a scale where, in direction k, W is essentially flat. This is also a characteristic distance on which models (projected on direction k) are indistinguishable. From the linear response perspective, 1 λ k tells us how much a small variation in the estimation of average of monomials (a linear combination of these variations parallel to v k ) affects the estimation of h l . Thus, following Mastromateo [START_REF] Mastromatteo | On the typical properties of inverse problems in statistical mechanics[END_REF] two distributions µ (1) , µ (2) are indistinguishable with accuracy if -log µ * h (1) = h (2) ≤ . If T is large enough, the set of indistinguishable distributions defines an ellipsoid s in E with a volume V:

V = 1 √ det χ 1 Γ( L 2 + 1) 2π T L 2 . ( 13
)
Thus, when the sample size T , the model dimension L and the accuracy are fixed the log of the volume is proportional to log V ∝ -1 2 L k=1 log λ k . Thus, model estimation is better if λ h s are larger, which implies that there exists a set of eigenvalues that can be neglected given their small magnitude, i.e. huge fluctuations impairing the model estimation. Then, instead of considering the volume of s in the space of all parameters, let us consider the volume V(k) of the projection of s in the subspace spanned by the k first eigenvectors of χ. We have:

log V(k) = 1 2 S φ (k) -log Γ( k 2 + 1) + k 2 log 2π T , (14)
with:

S φ (k) = - k i=1 log λ i , (15)
where λ i s are ordered by decreasing values. In eq. ( 14), The second term, log Γ( k 2 + 1) is purely combinatory, whereas the third term, log 2π T characterizes the effect of precision accuracy and finite sampling. Therefore, the only term which depends on the statistical model (here the potential H) is S H (k). In particular, it depends on the number of neurons N and memory depth R.

Finding the cut-off, k c

As eigenvalues k increases, there will be one or more ks at which S φ (k) is expected to become bigger than the sum of the other two terms of [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF], which are both negative. This means that, for increasing k, log V(k) will first decrease then it will increase, yielding at least one minima on the function. There is therefore a k, k c ≡ k c ( , T ), which characterizes the number of eigenvectors (i.e. dimension) for which the volume is minimal. Precisely, as we observed, belongs to some interval; outside this interval log V(k) is not convex any more. On this range over we obtain a set of k c ≡ k c (T ) minimizing the volume, which characterize the number of eigenvectors ensuring an optimal fit with the model: the model indeterminacy is minimized. Among the set of possible k c , we chose the smallest one, representing the minimum set of dimensions necessary to reliable fit the model given data. The corresponding eigenvectors span a subspace E r of E, embedded in the tangent space of W at the point corresponding to the empirical measure. This subspace corresponds to linear combination of h l s giving the most reliable model. Getting outside this space, by adding/subtracting degree of freedom/eigenvectors leads to a worse approximation of the empirical distribution.

PLOS

  (i) In the MEM each interaction (monomial m l ) has a weight h l . The set of all weights fixes the statistics predicted by the model. In particular it fixes the predicted average value of m l , µ[m l ]. A slight variation δh l of the sole weight h l induces both a variation δµ[m l ] and δµ[m l ] of µ[m l ] and µ[m l ], respectively, i.e. a change on the parameter controlling a given monomial average also affects other monomial average. One can show that δµ[m l ] = χ ll δh l . Thus, χ is a square matrix whose dimension, L, is the number of parameters fixing the model.

Finally, k c is

  conditioned to the accuracy value as shown on Fig 1D. Increasing or decreasing one order of magnitude yields no convex functions, where the minima would be at k = 1 in the former and k = L in the later. This constrains the search of the minima to a subset of values. Extending the results for different values of shows that the number of dimensions related to the minima decays monotonically with the accuracy ( ) until it reaches an inflection point; this inflection point is the k c that we use as the number of relevant dimensions (Fig S1).

  where κ is the compressibility of the code and N.O. are the Not Observed events. [MJ: What about the recording length?, Fig S2A and S2B?]

  6 k c and L are still very close to each other (Fig S2D). [MJ: missing to talk about Fig S2C] [MJ: I also suggest include Fig S2 as a main figure of the article]

Five

  different versions of the PWR2 underlying statistics were generated by scaling the parameters of the MEM by a factor f ac = {0.4, 0.6, 0.8, 1.2, 1.4}. The resulting parameters distributions (split by firing rates and pairwise interactions) and the generated monomials probabilities are shown in Fig ??A and Fig ??B, respectively.Recall that very negative parameters defines very unlikely monomials, while very positive the opposite. Parameters very close to zero define events which occurrence probability is near to 0.5.Values of f ac < 1 reduce the magnitude of the MEM parameters, but it increases both the firing rates and the pairwise interactions probabilities with negative associated parameters, compared to the original PWR2 raster, i.e. denser activity (density as the total number of spikes in a raster respect to the recording length) (Fig ??B).Multiplying the MEM parameter by values of f ac > 1 generates a sparser activity with low firing rates probabilities, as expected from very negative rates parameters, and many pairwise interactions vanish (the reduced width of the plot, see Fig ??B), given the low firing rates and the big amount of negative pairwise interactions parameters (Fig ??A).

Fig 1 .

 1 Fig 1. Dimensionality reduction framework overview.A: shows a raster representing the binary activity of N = 4 neurons (rows) over time (T ). Green square shows a slice of this raster, which is considered as 1 (or 0) if it generated (or not) activity inside a time window of width equal to the bin size. For this raster 2 types of pairwise interactions are defined: spatial interactions (blue) and temporal interactions with R = 2, i.e. one time-step between spikes (red). Spatial interaction is exemplified as w 3 (0)w 4 (0) (both neurons firing at the same time) and w 3 (1)w 4 (0) (neuron 3 firing one bin after neuron 4). B: Susceptibility matrix eigenvalue spectrum in log-log scale for Independent (black) and Pairwise Range=2 (PWR2, red). Black vertical dashed line denotes the network size (N = 20). Indep. raster shows a sharp cut-off close to the network size, while PWR2 shows a monotonic decay of the eigenvalues magnitude, without a clear cut-off. C: Illustrates a 3D system representing a volume of indistinguishable models living in x, y, z dimensions which is rotated into a new coordinate system x', y', z', that are linear combinations of the original coordinate system (original coordinate system denoted as dashed lines). Bottom shows the 2D projections, where the amplitude of the ellipsoid on each dimension is the inverse of the corresponding eigenvalue. In this case x' and y' are the stiff directions, while z' is a sloppy direction. D: Shows the value of the volume of indistinguishable models, V(k, ) as the degrees of freedom (k) increases. The volume reaches a minimum (orange dot) for certain k named k c . Increasing k upon k c makes the volume explode. Increasing/decreasing one order of magnitude the imposed accuracy ( ) yields no convex functions, i.e. no minimization.

Fig 2 .

 2 Fig 2. Dimensionality reduction by minimization of the volume of indistinguishable models.A: shows the log of volume of indistinguishable models (logV) as function of the number of dimensions (k) at fixed accuracy ( ). Thick solid lines are averages, shaded area is ±1 s.d. of the 100 temporal subsamples (of a time width equal to half of the recording length) and black vertical dashed line is the network size (N = 20). The minimum of this function (squares) is the optimal number of dimensions capturing the raster statistics at the given accuracy, i.e. k c . Inset shows a zoom-in on the first 100 dimensions, focusing on the Indep. minima. B: Summary of the number of dimensions for both rasters. We found a close relationship between the underlying model dimensionality and the more dimensions required for an optimal model. k c is 19 ± 0, and 447 ± 7.00 (mean ± std) for Indep. and PWR2, respectively. C: Is the same than (B), showing k c as a percentage of the underlying model dimensionality. Indep case show almost 100%, while PWR2 shows ± 75% which is due to the unobserved events.

Fig 3 .

 3 Fig 3. Dimensionality reduction on scaled PWR2 statistics. A: shows the underlying parameters distribution, split by firing rates (blue) and pairwise interactions (red). The bigger the scaling factor fac, the more negative the rates parameters and the more wide the interactions parameters distribution. B: shows the corresponding monomials probabilities for the scaled rasters, split as firing rates and interactions, as A. We see that increasing fac has the effect of decreasing both rates and pairwise interactions probabilities, reaching the point where many pairwise interactions vanish (f ac > 1). C: are the corresponding χ eigenvalue spectra (average out of 10 rasters) for the scaled rasters. Increasing f ac decreases the first eigenvalue, the spectrum offset and the number of eigenvalues above the minimal observed probability (1/T). D: k c values for the scaled rasters. None of the scaled rasters shows k c value as the one obtained for the original PWR2 raster. Dost are averages and errorbars 1 s.d. of 10 different rasters with the same underlying parameters. E: k c values plus the number of unobserved events. For f ac > 1 adding the unobserved events yields values close to the original PWR2 raster, showing that the dimensionality reduction obtained for those rasters is given mainly by the unobserved effects. For f ac < 1 we see no unobserved events, so the dimensionality reduction obtained for those cases is given mainly by the increased density of the raster hiding the underlying statistics.

k

  c for empirical data Following our work flow, we minimized the volume and computed k c for all the rasters (experimental and shuffled) at 2 different time scales (bin size): fast (1 and 5 ms) and medium (10 and 20 ms. Larger values of bin sizes were discarded because they significantly reduced the total bin number biasing the k c estimation). As a global picture, the stimuli spatio-temporal modulation increases k c for all bin sizes (i.e. comparing PSA with the other two stimuli) (Fig ??C).

  The eigenvectors v k of χ constitute an orthogonal basis, i.e. linear combinations of model parameters having the following interpretation. Once the set of interactions defining the model has been fixed, each set of parameters (h l )

	(iii) L l=1 defines a
	probability i.e. assigns a definite value to the predicted probability of events. The
	set of (h l ) L l=1 is a vector of L dimensions, i.e., a point in R L . Obviously, the
	closest two points are in the space, the closest are the statistics they predict. As a
	corollary, trying to fit empirical statistics, two models corresponding to "close"
	points might be indistinguishable i.e. they describe equally well the empirical
	statistics. The notion of closeness is however ambiguous here and requires to
	define a proper metric. This is precisely what χ does: it defines a metric (called
	Fisher metric).
	An example of spectra obtained for synthetic data is
	presented in Fig 1B.

  (t) which is 1 if and only if neuron i fires at time t in spike train ω. In (3), the monomial value depends in fact on the R -1 first time steps in the raster, i.e. on the block ω R-1

	0

One fixes a set of pairs (i k , t k ) (neuron, time), and m l (ω) = 1 if and only if neuron i k spikes at time t k in the raster block ω, and m l (ω) = 0 otherwise. The simplest example PLOS 15/22 is ω i

Note that the positivity holds true for the susceptibility defined exactly by[START_REF] Vasquez | Gibbs distribution analysis of temporal correlation structure on multicell spike trains from retina ganglion cells[END_REF], but its estimation from a finite raster can have small negative eigenvalues. When this happened, negative eigenvalues were removed from the spectrum.PLOS

Acknowledgments

Financial support: CONICYT-FONDECYT 1140403 and 1150638, CONICYT-Basal Project FB0008, Millennium Institute ICM-P09-022-F, ECOS-Conicyt C13E06, ONR Research Grant #N62909-14-1-N121.

Synthetic Raster Generation

Synthetic rasters (T = 2 * 10 6 time-points, N = 20 neurons) were generated using different underlying statistics: Independent, where only firing rates are defined, L = N ; Pairwise R=2 (PWR2), with firing rates and spatio-temporal correlations, L = N (3N -1)/2. Using the underlying potential of the PWR2, 5 more rasters were generated scaling the magnitude of the model parameters by a factor f ac = [0.4 0.6 0.8 1. 2 1.4]. This scaling mimics the effect of the stimuli high-order statistics, which increases both the firing rate and the spatio-temporal correlations. In addition, 6 more random PWR2 rasters were generated with N = [30 40 50 60 70 80] to study the dependence between k c estimation and the network size. Each raster was generated using the ENAS software (https://enas.inria.fr/). For the first 3 rasters, 100 random subsamples with half duration of the whole recording were taken for each raster and the χ matrix associated with a Pairwise R=2 model was computed. Then, k c was found by volume minimization, yielding 100 k c values per raster. For the scales rasters, the same procedure was applied, but using 10 temporal subsamples.

Shuffling

In order to destroy the dependencies between the empirical raster monomials, we generated random raster where the number of neurons and firing rates was exactly the same than observed on the recordings (i.e. on each retina under each condition), but the spikes times were otherwise random, avoiding violations of the refractory period (2ms).

Recordings and Stimuli

Animals and Recordings

4 Adult male and female Octodon degus (3-6 months) were maintained in the animal facility of the Universidad de Valparaiso, at 20-25°C on a 12-h light-dark cycle, with access to food and water ad-libitum. The methods of MEA recording has previouslly been described [?]. In brief, experiments were approved by the bioethics committee of the Universidad de Valparaiso, in accordance with the bioethics regulation of the Chilean Research Council (CONICYT) and international protocols. Animals were euthanized under deep isofluorano or halothane anesthesia and both eyes were extracted. Then, one of the extracted retinas was diced into quarters while the other was stored in oxygenated in oxygenated (O 2 95% CO 2 5%) AMES medium at 32°C in the dark for further experiments. The same AMES media was used for continuous perfusion during extracellular recordings. For MEA recordings (MEA USB-256, 20kHz sample, Multichannel Systems GmbH, Germany), one piece of retina was mounted onto a dialysis membrane placed into a ring device mounted in a traveling (up/down) cylinder, which was moved to contact the electrode surface of the MEA recording array. Data were processed off-line using Plexon Offline Sorter (Plexon Instruments). Further, spikes were detected using a threshold of -4.5 to -5.5 S.D. from the mean voltage value and then were manually classified using the 2D space of the first two principal components on each electrode. Only somatic spikes were kept. Refractory period violations were detected and discarded if two or more spikes of the same neuron occur in a 2ms period. We recorded on 3 stimuli conditions (see next section): i) Photopic Spontaneous Activity (PSA), ii) Spatio-temporal white Noise (WN) and iii) Natural Movie (NM), obtaining 151, 200, 246 and 270 RGC for each of the 4 experiments, respectively. For each experiment 30 random subsamples of 50 neurons were taken and χ matrix corresponding to a Pairwise R=2 model were computed considering 4 bin sizes (1, 5, 10 and 20 ms) for each subsample, yielding 30 k c values per experiment, per condition and per bin size. The shuffled version of these rasters were submitted to the same analysis.

Visual Stimuli

Visual stimuli were generated by a custom software created with PsychoToolbox (Matlab) on a MiniMac Apple computer and projected onto the retina with a LED projector (PLED-W500, Viewsonic, USA) equipped with an electronic shutter (Vincent Associates, Rochester, USA) and connected to an inverted microscope (Lens 4x, Eclipse TE2000, NIKON, Japan). The image was conformed by 380 x 380 pixels, each covering 5µm 2 . Since rodents are dichromatic (green and blue/UV cones), in our experiments only the B (blue) and G (green) beams of the projector were used, while the R (red) channel was used for signal synchronization. Dark spontaneous activity where recorded in order to monitor the stabilization of the activity. The, stimuli where applied. For PSA a space-time invariant stimuli with G and B intensities equal to the mean intensity of the NM stimulus were presented for 15 mins. WN stimulus with a block size of 50µm was used at rate of 60 fps and presented for 20 mins, with each block taking independently 0 or 255 (max value) in the pixel value scale. NM consisted on a 1800 frames movie recorded on the natural habitat of the rodent using a robotic solution to capture the natural visual environment of degus, including grass, trees, optic flow, head-like movements. This short movie were presented 40 times at a refresh rate of 60fps, yielding a total duration of 20 mins. Optical density filters in the optical path were used to control final light intensity. A CCD camera (Pixelfly, PCO, USA) attached to the microscope was used for online visualization and calibration of the light stimuli projected onto the recording array.

Supporting Information