Vincent Gouldieff

Jacques Palicot
email: jacques.palicot@supelec.fr

Steredenn Daumont
email: steredenn.daumont@zodiacaerospace.com

BLIND DIGITAL MODULATION CLASSIFICATION BASED ON M-TH POWER NONLINEAR TRANSFORMATION

Keywords: Automatic Modulation Classification, M th -Power nonlinear Transform, Blind Detection

Automatic Modulation Classification (AMC) has received a major attention last decades, as a required step between signal detection and demodulation. In the fully-blind scenario, this task turns out to be quite challenging, especially when the computational complexity and the robustness to uncertainty matter. AMC commonly relies on a preprocessor whose function is to estimate unknown parameters, filter the received signal and sample it in a suitable way. Any preprocessing error inherently leads to a performance loss. To improve the robustness of the blind AMC, we propose to proceed almost directly on the received signal -with neither matched-filtering step nor synchronization step. In this paper, Analytical M th -Power nonlinear Transformation (AM PT) is considered for its robustness towards timing, phase and frequency uncertainty. The generated feature-vector then feeds a Minimum Distance classifier. Numerical simulations show the effectiveness of the proposed method for a 7-class problem of low-order modulations.

INTRODUCTION

The demand for reconfigurable radio technologies -such as Cognitive Radio (CR) -has led to an increasing need for blind signal processing, such as spectrum sensing [START_REF] Axell | Spectrum Sensing for Cognitive Radio: State-of-the-Art and Recent Advances[END_REF], standard recognition [START_REF] Palicot | A New Concept for Wireless Reconfigurable Receivers[END_REF] and feature detection [START_REF] Khalaf | Blind Cyclostationary Feature Detector based on Sparsity Hypotheses for Cognitive Radio Equipment[END_REF]. Since the early 90', Automatic Modulation Classification (AMC) has been extensively studied as it represents a key challenge in Cognitive Radio. In recent years, several techniques for AMC have been reported in the literature [START_REF] Dobre | Survey of Automatic Modulation Classification Techniques: Classical Approaches and New Trends[END_REF]. Most of these methods rely on a two-steps architecture: a signal preprocessing task followed by a recognition algorithm (see e.g., Fig. 1).

The first step generally consists of a demodulator which feeds the recognition algorithm with a filtered, normalized and symbolspaced sampled signal. Therefore, numerous inherent parameters of the received signal -such as the symbol rate, pulse shaping, carrier frequency, signal and noise powers -are generally assumed to be perfectly known at the receiver. This strong assumption is in contradiction with the blind property needed in a smart receiver.

The second step (the modulation recognition algorithm itself) may rely on a Likelihood-Based (LB) or a Feature-Based (FB) method [START_REF] Dobre | Survey of Automatic Modulation Classification Techniques: Classical Approaches and New Trends[END_REF]. While the former is optimal in the Correct Classification Rate (CCR) sense, the latter is often favoured for its near-optimal performance at a lower complexity. As a consequence, we propose in the scope of this article a new Feature-Based classifier. Obviously, the quality of the whole blind algorithm depends on both the performance of the preprocessor and of the recognition algorithm. FB classifiers rely on two sub-steps, namely a feature computation step and a decision making step. Signal statistics-based features, for instance Higher Order Statistics (HOS), have been widely studied for their ability to differentiate between digital modulations. Cyclostationarity is also a good candidate in the blind scenario for its robustness to frequency, phase and time offsets [START_REF] Hazza | An Overview of Feature-based Methods for Digital Modulation Classification[END_REF]. Numerous decision algorithms were proposed, as Artificial Neural Network (ANN), Support Vector Machine (SVM) or Decision Tree (DT) [5][6].

Swami et al. proposed a classifier relying on Higher Order Cumulants (HOCs) and DT to recognize PAM, QAM and PSK modulations [START_REF] Swami | Hierarchical Digital Modulation Classification using Cumulants[END_REF]: this method is indeed efficient but unsuited in the blind scenario. As a step towards blindness, some improved HOS-based methods were proposed to face unknown signal power [START_REF] Dai | Joint Power Estimation and Modulation Classification Using Second-and Higher Statistics[END_REF] or carrier frequency offset [START_REF] Chaithanya | Blind Modulation Classification in the Presence of Carrier Frequency Offset[END_REF] [START_REF] Flohberger | Modulation Classifier for Signals used in Satellite Communications[END_REF]. As expected, the CCR is improved with these methods when the preprocessing task is assumed non-ideal. Cyclostationarity was also proposed to detect spectrum holes [START_REF] Khalaf | Blind Cyclostationary Feature Detector based on Sparsity Hypotheses for Cognitive Radio Equipment[END_REF] and classify analog and digital modulations [START_REF] Dobre | Cyclostationarity-Based Blind Classification of Analog and Digital Modulations[END_REF]. Its main drawback is the complexity entailed by the computation of the cyclic function.

Recently, several fully blind AMC techniques based on the M th -Power nonlinear Transformation (M PT) were suggested [START_REF] Xie | An Efficient and Simple Method of MPSK Modulation Classification[END_REF][START_REF] Lim | Automatic Modulation Recognition for Spectrum Sensing using Nonuniform Compressive Samples[END_REF][START_REF] Xing | Automatic Modulation Recognition of PSK Signals using Nonuniform Compressive Samples based on HOS[END_REF]. M PT was first studied as a tool for carrier and timing synchronization [START_REF] Je Mazo | Jitter Comparison of Tones Generated by Squaring and by Fourth-Power Circuits[END_REF] and then used for AMC in the 90' [START_REF] Reichert | Automatic Classification of Communication Signals using Higher Order Statistics[END_REF]. M PT is of major interest since the produced features are insensitive to time, frequency and phase offsets: thus, no elaborate preprocessing step is needed. In most of the previous work, the recognition only relies on the number of peaks generated by the M PT for different M . Then, the classifier is unsuited for an efficient classification if the test-set contains some constellations of the same order.

We propose a novel fully-blind method for AMC, based on the theoretical power of the peaks generated by the M PTs. This technique -we called Analytical M PT (AM PT) -is suitable for PSK, PAM, QAM and hybrid modulations. Moreover, the preprocessing task is the simplest possible to improve the robustness of the whole classifier to estimation uncertainty. A simple MD classifier is considered for a better fusion of the different powers in the feature-space.

The paper is organized as follows: In Section 2 the problem under study is formulated, while in Section 3 the proposed blind AMC method is described. In Section 4 the theoretical development of the references is carried out. The performance is studied via extensive simulations in Section 5.

PROBLEM FORMULATION

Let us assume that the signal-of-interest consists of a pseudobaseband pulse-shaped digital signal in AWGN channel. This scenario is relevant, for instance, for fixed satellite communication systems. The signal is then modelled as

x(n) = a • e i•(2πfr •n+φ) • k s(k) • h(nTe -kT -τ ) + ω(n) (1)
where a is the amplitude of the signal, fr is the residual carrier frequency offset, φ is the initial phase, s(k) is the unknown sequence of i.i.d. symbols, h is the pulse shaping function, Te is the sampling period, T is the symbol period, τ is the initial delay and ω is the sampled and unfiltered white Gaussian noise.

We consider that h is a root raised cosine (RRC) filter with unknown roll-off β and that the oversampling factor ρ = T Te respects the Shannon's theorem i.e., ρ ≥ 1 + β. Remark that we do not assume any time or frequency synchronization here. We denote by N the number of considered symbols. Thus, we dispose of Ne = ρN samples for the classification. We assume that all the parameters in (1) are constant on [0; Ne -1].

We consider the 7-class problem of C = {BPSK, QPSK, 8PSK, 8AMPM, R8QAM, C8QAM, 16QAM} in order to show the effectiveness of the proposed method on quite similar zero-mean and unitvariance constellations (see Fig. 2).

Fig. 2. Representation of the considered zero-mean and unitvariance constellations in the complex plan

DESCRIPTION OF THE PROPOSED AMC METHOD

Basic definitions

In the scope of the paper, we define the M PT of x(n) as the M throot of the Power Spectral Density (PSD) of x M (n)

M PTx(f ) = Ns -2 M • Ne-1 n=0 x M (n) • e -2iπnf 2 M (2)
This transform may show spectral lines around the null frequency for some values of M : their existence and power mainly depend on the emitted constellation.

While other M PT-based techniques in the literature only exploit the existence property [START_REF] Xie | An Efficient and Simple Method of MPSK Modulation Classification[END_REF][START_REF] Lim | Automatic Modulation Recognition for Spectrum Sensing using Nonuniform Compressive Samples[END_REF][START_REF] Xing | Automatic Modulation Recognition of PSK Signals using Nonuniform Compressive Samples based on HOS[END_REF][START_REF] Reichert | Automatic Classification of Communication Signals using Higher Order Statistics[END_REF] -and then fail to differentiate between quite similar constellations such as QPSK and 16QAM, or BPSK and R8QAM -we instead propose to exploit the power of the spectral lines as strong features. We then dispose of more useful information, which also further enables to derive a soft decision.

Proposed preprocessing unit

In the fully-blind scenario, some preprocessing tasks have to be carried out. First, a spectrum sensor (see e.g., [START_REF] Axell | Spectrum Sensing for Cognitive Radio: State-of-the-Art and Recent Advances[END_REF]) followed by a standard recognizer (see e.g., [START_REF] Palicot | A New Concept for Wireless Reconfigurable Receivers[END_REF]) respectively locate the bands where the AMC have to be performed and provide information for further detection stages. Thereby, these tasks are assumed correctly performed without affecting the blind property of the proposed method. Then, a simple normalizer and parameter estimator proceed as follows:

Normalization

Assume that no interferer is present in the -3dB bandwidth of the received signal, it is clear that

P 3dB = Pω • ρ -1 + Pu • P H,3dB
, where Pω is the power of the noise, Pu is the power of the noiseless signal-of-interest u(n) = x(n)ω(n), and P H,3dB is the power of the pulse-shaping filter in its -3dB bandwidth. With a root-raised cosine filter, we have P H,3dB = 1 -β 2 + β π and we eventually get

Pu = P 3dB -Pω • ρ -1 1 -β 2 + β π ( 3
)
If not known, remark that β, ρ and Pω can be easily estimated using the formerly computed PSD of x(n).

Residual Carrier Frequency Estimation

The residual carrier frequency fr also has to be estimated in order to locate the spectral lines -which occur at frequencies fM = M fr.

Such an estimation is generally performed by a maximum-tomean detector applied on the computed M PTs:

M fr = argmax f,M M PTx(f ) -M PTx(f ) (4)
where M PTx(f ) is the mean value of M PTx(f ) around the peak.

Output of the preprocessor

The preprocessor feeds the feature-computer with the normalized version of x(n) and with the set of estimated parameters (Fig. 1):

y(n) = x(n) • Pu -1 2 and η = β, ρ, Pω (5)
In realistic conditions, band-pass pre-filtering is considered as well for a better rejection of out-of-band noise and adjacent signals.

Decision Making via Minimum Distance classification

We use a Minimum Distance (MD) classifier to fusion the features of set MiPTy(Mi fr), Mi ∈ M and compare them with the theoretical references. Denote by Π the feature-vector

Π = [M1PTy(M1 fr) • • • MmPTy(Mm fr)] (6)
and by Π C th the theoretical feature-vector if C is emitted (please refer to the next Section for the derivation of the references)

Π C th = [M1PT C th • • • MmPT C th ], (7)
the distance is simply computed with norm-2 as

∆ C = ||Π -Π C th ||2 (8)
We finally get the estimated constellation C as

C = argmin C∈C ∆ C (9)

Description of the proposed algorithm

The proposed fully-blind whole AMC algorithm is summed up in Fig. 3, where the general structure of Fig. 1 is conserved.

THEORETICAL REFERENCES

In this section, we derive the theoretical power of the spectral lines occurring in the M PT function. Since the PSD of a cyclostationary signal is insensitive to phase and time offsets [START_REF] Reichert | Automatic Classification of Communication Signals using Higher Order Statistics[END_REF], we assume that φ = 0 and τ = 0 in (1). Moreover, without loss of generality and for ease of notation, we assume that the residual frequency fr is null.

We first work on a noise-free, continuous-time and infinite length signal. The derivation of the references in realistic conditions is carried out in future work. In this section, we then consider:

u(t) = k s(k) • h(t -kT ) (10)
where the notations were formerly defined in Section II.

For the sake of a better understanding, we derive the theoretical references when M = 2 and M = 4. References for other values of M follow the exact same principle.

Computation of 2P T C th

As developed in [START_REF] Je Mazo | Jitter Comparison of Tones Generated by Squaring and by Fourth-Power Circuits[END_REF], u 2 easily splits into two sums

u 2 (t) = u20(t) + u11(t) (11)
where

u20(t) = k s 2 (k) • h 2 (t -kT ) (12)
u11(t) = k =k ′ s(k) • s(k ′ ) • h (t -kT ) • h t -k ′ T (13)
The presence of a frequency peak in the PSD of u 2 (t) at null frequency is made obvious by applying the Poisson Summation Formula on fictive signal u20(t). In fact, we have

E [u20(t)] = E[s 2 ] T m e i•m•t • H2 (m) (14)
where m corresponds to the multiples of the symbol frequency and H2(f ) corresponds to the Fourier Transform of h 2 (t). Note that u11(t) doesn't produce any spectral line for i.i.d. symbol sets since E[s(k)s(k + n)] = 0 whatever n = 0. In the scope of the paper, h 2 (t) corresponds to a raised cosine (RC) filter. Then H2(0) = T and the theoretical reference when M = 2 is finally derived as:

2PT C th = E[s 2 ] ( 15
)
Theoretical values of 2PT C th are tabulated as follows:

C → BPSK QPSK 8PSK 8AMPM R8QAM C8QAM 16QAM 2PT C th 1 0 0 0.2 2/3 0 0 Table 1. Tabulated values of 2PT C th for C ∈ C.

Computation of 4P T C th

The same development is carried out for the fourth power of u(t), which is plainly developed as

u 4 (t) = u4000(t) + u3100(t) + u2200(t) + u2110(t) + u1111(t)
where the notation in the right-hand side terms is straightforward when compared to ( 11) -( 13). Now, both terms u4000(t) and u2200(t) have an impact on the power of the spectral lines while the remaining terms don't (they actually show null expectations). Then, following the method of IV. A and after some simplifications in the indices, we may get:

4PT C th = 1 T E[s 4 ]H4(0) + 6 E[s 2 ] 2 k>0 H (k) 22 (0) 1 2 (16)
where H4(f ) stands for the Fourier Transform of h 4 (t) and

H (k) 22 (f ) stands for the Fourier Transform of h 2 (t) • h 2 (t -kT ).
Note that no simple analytical formula can be found neither for H4(0) nor for H (k) 22 (0) and that they depend on β. However, they can be easily computed numerically. A few values of 4PT C th as a function of β are gathered in Table 2.

β ↓ C BPSK QPSK 8PSK 8AMPM R8QAM C8QAM 16QAM 0.

Representation in the feature-space

Fig. 4 shows the theoretical and simulated feature-vectors in the feature-space. Note that a proper noise-corrective term is added to the references so that they better match the simulated features. Obviously, some sets of constellations will be easier to distinguish than other sets, as Fig. 4 and Table 3 show. Note that Table 3 gives the distance between constellations C and

C ′ as ∆ C,C ′ = ||Π C th -Π C ′ th ||2.
In future work, it would be interesting to study the probabilistic behaviour of the features. Then, we could derive the theoretical Correct Classification Rate (CCR) of the proposed AMC method. Herein, we present numerical results in the next section.

NUMERICAL RESULTS

Extensive simulations were carried out to show the effectiveness of the proposed method. The performance of our algorithm is also compared with the CCRs obtained with Swami's C40-based classifier ( [START_REF] Swami | Hierarchical Digital Modulation Classification using Cumulants[END_REF], Ex. 1) in both the blind and the non-blind scenario. The experiments were performed on 10000 trials per constellation. Recall that N is the number of symbols used for the classification.

Performance of the proposed classifier

Table 4 shows the Confusion Matrix (CM) for the proposed 7-class problem in the blind scenario. For a better reading, " • " means 0%.

As expected, the worst CCR is obtained for QPSK/C8QAM since these constellations are the closest in the 2PT/4PT feature-space.

Considering other features -such as higher-order or conjugated M PTs -would improve the distance between these constellations. However, increasing M may also involve a degradation of the CCR at low SNRs since the features would show a more stochastic behaviour as M grows and as the SNR weakens. Fig. 5 shows the CCR in the blind scenario for the proposed 7class problem as a function of the SNR. In this context, we have set N = 1000 and ρ = 2. As expected, set {BPSK, R8QAM, 8PSK} is the easiest to classify, while the CCR for {QPSK, C8QAM} remains slightly lower for medium-to-high SNR (see especially Table 4).

Performance comparison with the C40-based classifier

We compare the CCR obtained by Swami's C40-based method in [START_REF] Swami | Hierarchical Digital Modulation Classification using Cumulants[END_REF] -commonly used as a benchmark in the literature -with the proposed AM PT-based classifier, in the same conditions as Fig. 5.

Mean performances of the methods in both the blind and the non-blind scenario are depicted in Fig. 6. The A2PT/A4PT-based classifier shows a slight increase in the performance of the classification in the non-blind scenario. Besides, in the fully-blind scenario, the A2PT/A4PT-based classifier clearly outperforms Swami's method. This is mainly due to a better robustness of the proposed features to inherent estimation errors at the preprocessor.

Complexity comparison with the C40-based classifier

Last, we compare the 2PT/4PT-based method with the C40-based classifier in terms of complexity. In both cases, the greediest operation is the computation of the features.

The computation of M PT(M fr) (2) asymptotically requires M Ne complex multiplications and Ne complex additions if the exponential was precomputed. Then, the 2PT/4PT-based classifier requires 6Ne complex multiplications and 2Ne complex additions. In comparison, the computation of C40 is slightly lower since it requires 4Ne complex multiplications and 2Ne complex additions.

CONCLUSION

Throughout this article, a new efficient way to blindly classify digital modulations is described. The method is based on an Analytical study of the M th -Power nonlinear Transformation. The featurevector used for the classification is insensitive to time, phase and frequency offsets, which makes the whole classifier more robust to estimation uncertainty in the blind scenario.

The key features of the proposed method are the absence of complex preprocessing and its low complexity. These properties make it a great tool in smart receivers for wide-band monitoring and for the estimation of the features of the signal before the demodulation step.

Further work shall be carried out, for instance regarding the statistical distribution of the feature-vector in realistic conditions. The effects of multi-path channels on the analytical development of the references shall also be investigated.

Fig. 1 .

 1 Fig. 1. A common Feature-Based AMC architecture

Fig. 3 .

 3 Fig. 3. The proposed fully-blind AMC algorithm

Fig. 4 .

 4 Fig. 4. Theoretical and simulated feature-vectors, for 100 realizations per constellation, SNR = 15 dB, β = 0.3, ρ = 2 and N = 1e 4

Fig. 5 .Fig. 6 .

 56 Fig. 5. Correct Classification Rate as a function of SNR

Table 2 .

 2 Tabulated values of 4PT C th for C ∈ C and for different β.

	15 1.233 0.856 0.000	0.734	0.770	0.929	0.713
	0.25 1.200 0.881 0.000	0.754	0.743	0.961	0.737
	0.35 1.173 0.907 0.000	0.774	0.719	0.991	0.760

Table 3 .

 3 Distance Matrix in the proposed feature-space

						Constellation
			BPSK QPSK	8PSK 8AMPM R8QAM C8QAM 16QAM
		BPSK	0	1.004	1.419	0.905	0.564	0.988	1.055
		QPSK		0	0.761	0.219	0.656	0.068	0.159
	Constellation	8PSK 8AMPM R8QAM C8QAM			0	0.639 0	0.861 0.468 0	0.829 0.268 0.677 0	0.602 0.161 0.629 0.227
		16QAM							0

Table 4 .

 4 Confusion Matrix with SN R = 10 dB and N = 1000 for the proposed 7-class problem and in the blind scenario


Wahiba U Taouali
email: wahiba.taouali@univ-amu.fr

Giacomo U Benvenuti

Pascal U Wallisch

Frédéric U Chavane

Laurent U Perrinet

Testing the Odds of Inherent versus Observed Over-dispersion in Neural Spike Counts Odds of Inherent versus Observed Over-dispersion

Keywords: Spike counts, over-dispersion, Negative-Binomial distribution, Decoding, Tuning function

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

Introduction

One of the most common ways to study the central nervous system is to use a stimulus-response paradigm. A system, which could be a neuron or a population of neurons, is provided with an input (a stimulus) and the resulting state change (response) of the system is quantified in order to better characterize functional associations. A main problem with this approach is that the system's response to repeated presentations of an identical visual stimulus often exhibits high variability, particularly in the spike count statistics (see Figure 1). As a consequence, statistical methods are essential to characterize the information carried by neural populations.

The Fano factor [START_REF] Fano | Ionization yield of radiations. II. The fluctuations of the number of ions[END_REF] is commonly used to quantify the variability of spike counts [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF][START_REF] Alitto | A comparison of visual responses in the lateral geniculate nucleus of alert and anaesthetized macaque monkeys[END_REF][START_REF] Zhuang | Layer 4 in primary visual cortex of the awake rabbit: contrasting properties of simple cells and putative feedforward inhibitory interneurons[END_REF][START_REF] Tolhurst | The statistical reliability of signals in single neurons in cat and monkey visual cortex[END_REF][START_REF] Softky | The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs[END_REF][START_REF] Shadlen | The Variable Discharge of Cortical Neurons: Implications for Connectivity, Computation, and Information Coding[END_REF][START_REF] Koch | Predicting the visual world: silence is golden[END_REF]. It is defined by F F (∆t) = σ 2 (∆t) µ(∆t) where µ and σ 2 are the average and the variance over multiple trials of the number of spikes in some time window of length ∆t, respectively. This variability is most often characterized by the Poisson spiking Model (PM). By the definition of this model, the spike count follows a homogeneous Poisson process with a rate parameter λ (the expected number of spikes that occur per unit of time). If k is the observed number of spikes fired by a neuron in response to a stimulus in a time window ∆t, the probability of such an observation is given by a Poisson distribution of parameter f = λ∆t, which gives the expected number of spikes per sample in the considered time window:

P (k) = f k e -f k! (1)
It follows that µ = f and σ 2 = f and thus that the Fano factor under the Poisson model has a value of 1, regardless of λ or ∆t. Such a property is in accordance with a wide range of observations of variability in cortical responses to visual stimuli [START_REF] Softky | The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs[END_REF][START_REF] Shadlen | The Variable Discharge of Cortical Neurons: Implications for Connectivity, Computation, and Information Coding[END_REF][START_REF] Britten | Responses of neurons in macaque MT to stochastic motion signals[END_REF][START_REF] Buracas | Efficient Discrimination of Temporal Patterns by Motion-Sensitive Neurons in Primate Visual Cortex[END_REF][START_REF] Mcadams | Effects of Attention on the Reliability of Individual Neurons in Monkey Visual Cortex[END_REF] showing variance-to-mean ratios of spike counts that equal or just exceed unity. However, for a large majority of reported evoked responses, the variability is higher than expected from a Poisson model, especially in mammalian visual cortex [START_REF] Tolhurst | The statistical reliability of signals in single neurons in cat and monkey visual cortex[END_REF][START_REF] Heggelund | Response variability and orientation discrimination of single cells in striate cortex of cat[END_REF][START_REF] Dean | The variability of discharge of simple cells in the cat striate cortex[END_REF][START_REF] Schiller | Short-term response variability of monkey striate neurons[END_REF][START_REF] Vogels | The response variability of striate cortical neurons in the behaving monkey[END_REF], whereas Sub-Poisson variability has been observed only under specific conditions [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF].

In fact, the variability of neural responses depends on different parameters. Some studies reported that the Fano factor values can vary depending on the brain region [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF][START_REF] Gur | High Response Reliability of Neurons in Primary Visual Cortex (V1) of Alert, Trained Monkeys[END_REF][START_REF] Kayser | Visual enhancement of the information representation in auditory cortex[END_REF]. Indeed, neural responses in higher cortical areas [START_REF] Maimon | Beyond Poisson: Increased Spike-Time Regularity across Primate Parietal Cortex[END_REF] or in other non-cortical sensory areas such as the retina [START_REF] Berry | The structure and precision of retinal spiketrains[END_REF] and the thalamus [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF] show substantially lower variability. It has been also shown that some neural features such as the refractory period [START_REF] Kara | Low response variability in simultaneously recorded retinal, thalamic, and cortical neurons[END_REF][START_REF] Berry | The structure and precision of retinal spiketrains[END_REF] lead to a decrease of variability. Moreover, the nature of the stimulus itself could contribute to the variability of neuronal responses. It has been suggested that increasing the complexity of the stimuli by using natural-like statistics increases the reliability [START_REF] Baudot | Animation of natural scene by virtual eye-movements evokes high precision and low noise in V1 neurons[END_REF][START_REF] Haider | Synaptic and network mechanisms of sparse and reliable visual cortical activity during nonclassical receptive field stimulation[END_REF][START_REF] Borst | Information theory and neural coding[END_REF]. In addition, it has been reported that the variability of evoked responses is not stationary even during a trial time-course and decreases after the presentation onset of most sensory stimuli [START_REF] Churchland | Stimulus onset quenches neural variability: a widespread cortical phenomenon[END_REF]. Finally, a recent study [START_REF] Nienborg | Decision-related activity in sensory neurons: correlations among neurons and with behavior[END_REF] on behaving animals explored the role of decision and other cognitive processes in shaping variability. In this paper, we mainly explore over-dispersion, that is, cases where there is greater variance in the spike count than might be expected by the Poisson model. Our hypothesis is that the Poisson process only accounts for the intrinsic noise induced by the spiking mechanism itself. Indeed, under the Poisson process hypothesis, the trial-to-trial variability only results from the stochastic properties of the neuron [START_REF] Carandini | Amplification of Trial-to-Trial Response Variability by Neurons in Visual Cortex[END_REF][START_REF] Mainen | Reliability of spike timing in neocortical neurons[END_REF][START_REF] Schneidman | Ion Channel Stochasticity May Be Critical in Determining the Reliability and Precision of Spike Timing[END_REF]. Thus, the PM can not account for other significant sources of variability including hidden contextual variables such as attention, varying cortical states [START_REF] Buracas | Efficient Discrimination of Temporal Patterns by Motion-Sensitive Neurons in Primate Visual Cortex[END_REF][START_REF] Arieli | Dynamics of Ongoing Activity: Explanation of the Large Variability in Evoked Cortical Responses[END_REF][START_REF] Tsodyks | Linking Spontaneous Activity of Single Cortical Neurons and the Underlying Functional Architecture[END_REF][START_REF] Kenet | Spontaneously emerging cortical representations of visual attributes[END_REF], perceptual effects [START_REF] Ress | Neuronal correlates of perception in early visual cortex[END_REF] and overt behaviors (such as the precise eye position) [START_REF] Gur | Response variability of neurons in primary visual cortex (V1) of alert monkeys[END_REF]. Thus, we consider an alternative compound distribution, called Negative-Binomial, that generalizes the Poisson distribution with a dispersion parameter that directly controls the ratio between the mean and variance. Such a doubly-stochastic model was recently proposed to model the variability of neural dynamics showing more complex variability than what is expected from a Poisson model [START_REF] Churchland | Variance as a signature of neural computations during decision making[END_REF][START_REF] Ponce-Alvarez | Stimulus-dependent variability and noise correlations in cortical MT neurons[END_REF][START_REF] Ponce-Alvarez | Comparison of local measures of spike time irregularity and relating variability to firing rate in motor cortical neurons[END_REF][START_REF] Goris | Partitioning neuronal variability[END_REF]. Our aim is to further explore the potential suitability of this doubly-stochastic model to account for response variability in neurophysiological data.

Θ =0 • Θ =22.5 • Θ =45 • 0 500 Time(ms) Trials Θ =67.5 • Θ =90 • Θ =112.5 • Θ =135 • 0 500 Θ =157.5 • Θ =180 • Θ =202.5 • Θ =225 • 0 500 Θ =247.5 • Θ =270 • Θ =292.5 • Θ =315 • 0 500 Θ =337.
The paper is organized as follows. First, we describe the Negative-Binomial encoding model (NBM) and two possible implementations of over-dispersion sources. Second, we propose a statistical test for evaluating whether the observed over-dispersion is significantly better described by the Negative-Binomial model than the Poisson model. To assess the generality of this model, we then apply this test to three data sets of spiking responses in different animals, different states and different visual stimuli: in the lateral geniculate nucleus (LGN) neurons of anaesthetized mice presented with drifting gratings, in the primary visual cortex (V1) of awake macaque monkeys viewing oriented moving bars and in the middle temporal cortex (MT) of anaesthetized macaque monkeys presented with drifting gratings. Finally, we compared the efficiency of PM and NBM models on a population decoder using the most dispersed data set (MT).

1

The Negative-Binomial model: encoding of inter-trial overdispersion Several models have been developed to account for over-dispersion such as generalized Poisson [START_REF] Consul | A generalization of the Poisson distribution[END_REF], zeroinflated Poisson [START_REF] Lambert | Zero-inflated Poisson regression, with an application to defects in manufacturing[END_REF][START_REF] Miaou | The relationship between truck accidents and geometric design of road sections: Poisson versus negative binomial regressions[END_REF][START_REF] Shankar | Modeling accident frequencies as zero-altered probability processes: an empirical inquiry[END_REF] or quasi-Poisson [START_REF] Mccullagh | Generalized linear models[END_REF][START_REF] Wedderburn | Quasi likelihood functions, generalized linear models, and the GaussNewton method[END_REF] models. However, only a few such as the Negative-Binomial model (NBM) are well known and widely available. The Negative-Binomial model has been studied and tested in different fields such as epidemiology [START_REF] Byers | Application of negative binomial modeling for discrete outcomes: a case study in aging research[END_REF], accident statistics [START_REF] Poch | Negative binomial analysis of intersection accident frequencies[END_REF], biology [START_REF] Anscombe | Sampling theory of the negative binomial and logarithmic series distributions[END_REF][START_REF] Bliss | Fitting the negative binomial distribution to biological data[END_REF], etc. However, it has only recently been applied to evoked neural responses [START_REF] Goris | Partitioning neuronal variability[END_REF][START_REF] Scott | Fully Bayesian inference for neural models with negative-binomial spiking[END_REF] where the variability was traditionally supposed to be well modeled by a Poisson model (PM).

A Negative-Binomial model depends on a distributional form equivalent to a compound stochastic process. It has recently been represented as a Polya-Gamma mixture of normals [START_REF] Polson | On the half-Cauchy prior for a global scale parameter[END_REF], but its standard representation is the Gamma-Poisson mixture. It corresponds to a doubly stochastic process where the Poisson process has a probability distribution P with a parameter λ which is itself a random variable generated from a Gamma distribution G(φ, β) with a shape parameter φ and a scale parameter β. When we consider a random variable x generated by this compound process (i.e. x ∼ P(λ) with λ ∼ G(φ, β)), the resulting distribution of this variable corresponds to a Negative-Binomial distribution (NB) which is usually defined by its shape parameter φ and its scale parameter p =1 β+1 :

x ∼ N B(φ, p = 1 β + 1 ) (2)
Using the Negative-Binomial model, the probability of getting k spikes from a given cell with parameters φ and p in response to a specific stimulus is given by:

P (k; φ, p) = Γ(k + φ) • p φ • (1 -p) k Γ(φ)k! (3)
We deduce the mean (µ), the variance (σ) and the Fano factor (F F ) as:

µ = φ • β, σ 2 = µ + µ 2 φ and FF = σ 2 µ = 1 + µ φ (4)
This distribution approaches the Poisson probability of parameter µ for large values of φ, but allows to reach larger-than-Poisson variances for small φ values. We end up with a multivariate distribution that is -in practice -driven by two parameters: the mean spike count (µ) and the parameter (φ) that we call the "inverse-dispersion parameter" 1 . Under the NBM hypothesis, the distribution of inter-spike intervals (ISIs) is given by a Loamax distribution (equivalent to an exponential distribution with a random parameter generated from a Gamma distribution). The mean and the variance of such a distribution are given by the following expressions: mean(ISIs) = β (φ-1) defined for φ > 1, variance(ISIs) =

β 2 φ (φ-1) 2 (φ-2) defined for φ > 2.
The resulting squared coefficient of variation is given by :

CV 2 = φ (φ -2) = 1 + 1 φ 2 -1 (5)
defined for φ > 2. CV 2 is higher than 1, which implies than a Negative-Binomial process is more irregular than a Poisson process and approaches Poisson regularity for large values of φ. However, unlike under the Gamma distribution hypothesis made in [START_REF] Nawrot | Measurement of variability dynamics in cortical spike trains[END_REF], the Negative-binomial model do not allow us to systematically relate over-dispersion to the positive deviation from the null-hypothesis F F = CV 2 . Thus, to better qualify this over-dispersion parameter, we need to focus on how to generate the Gamma input drive that feeds the Poisson generator.

Over-dispersion may result from extrinsic sensory noise: NBM as a mixture of Poisson and LogNormal distributions. Many studies assume a Log-Normal or a Gamma distribution when analyzing spike count distributions as it most often produces similar results [START_REF] Mccullagh | Generalized linear models[END_REF][START_REF] Atkinson | Regression diagnostics, transformations and constructed variables[END_REF]. Indeed, a Log-Normal process might be a good approximative model for a Gamma process. The underlying spiking process could result from a simple encoding model: Let us assume that a noisy sensory input is transformed to a response by a linear correlation measure with a cell's receptive field (for instance, in the early visual system, selective to a given orientation) and then through a non-linearity into a stochastic spike count (Figure 2). If we consider additive extrinsic Gaussian noise as the sensory noise and an exponential function as the non-linearity, the resulting input drive follows a Log-normal distribution that is closely fitted by a Gamma distribution. Using such a model, Figure 2 shows that over-dispersion could simply result from an additive normal sensory noise and the non-linearity of the spiking response. However, even if log-Normal and Gamma distributions are often used to model the same phenomena, looking at their log functions -in general -shows clear differences in their respective skewness. Furthermore, there is no exact morphing between the two distributions that allows to link their respective parameters analytically. As such, this approximation does not allow a precise quantification of the relationship between the noise source and the over-dispersion parameter.

Over-dispersion may result from redundancy within a neural population: NBM as a mixture of Poisson and Exponential distributions. The Gamma distribution could also result from a mixture of distributions. Indeed, it is well known that the sum of n independent and identically distributed exponential random variables E( 1λ ) (the parameter of an exponential distribution is by definition the inverse of its mean, here called λ) follows a Gamma distribution G(n, λ), where n and λ are its shape and scale parameters respectively. As such, exponential distributions have previously been used as a simple model of variability in formal neural networks [START_REF] Treves | What determines the capacity of autoassociative memories in the brain ?[END_REF][START_REF] Levy | Energy efficient neural codes[END_REF][START_REF] Baddeley | An efficient code in V1?[END_REF]. They correspond to a model of the distribution of time between two Poisson events. In addition, such a distribution maximizes the entropy of spike counts (knowing the average count) and thus maximizes the carried information [START_REF] Shannon | A Mathematical Theory of Communication[END_REF].

From a biological point of view, to drive a Gamma input -with an inverse-dispersion parameter φ equal input RFs

X φ 1 E( 1 λ )
Poisson(G(φ, λ)) This measures the correlation between the input Gabor and the Gabor representing the orientation θ characteristic of their receptive field (RF). Let us further assume that the filter cells exhibit an exponential distribution of spike counts E( 1 λ ) with mean λ, thus maximizing the entropy of the input drive. Then, the pooling of the resulting spike counts results in an input drive G(φ, λ) equal to the sum of exponential variables having the same mean equal to λ. It can be shown that the Poisson spiking process driven by this input drive is equivalently described by a Negative-Binomial probability distribution function with parameters λ (mean) and φ (inverse-dispersion) (see equation 2). Such a distribution describes observed trial-to-trial variability well and in this model, the parameter of dispersion is directly correlated with the redundancy, i.e., the convergence of similar input cells.

to n -to the output neuron considered as a Poisson generator, this neuron should receive as input the sum of n exponential variables having the same mean λ (see Figure 3). Thus, φ reflects the amount of convergence to a neuron which is -for instance -consistent with a model of the convergent connectivity from V1 to area V5/MT [START_REF] Simoncelli | A model of neuronal responses in visual area MT[END_REF][START_REF] Rust | How MT cells analyze the motion of visual patterns[END_REF][START_REF] Wang | Responses to second order texture modulations undergo surround suppression[END_REF]. In addition, it was observed that the simulation of networks based on biologically realistic parameters [START_REF] Voges | Phase space analysis of networks based on biologically realistic parameters[END_REF] could lead to complex dynamics, showing in particular an excess of variability [START_REF] Voges | Complex dynamics in recurrent cortical networks based on spatially realistic connectivities[END_REF]. Thus, we could expect an over-dispersion in V5/MT resulting from dimensionality reduction (i.e. projecting inputs from a high dimensional space to a lower dimensional one [START_REF] Haykin | Neural Networks: A Comprehensive Foundation[END_REF]).

To summarize, we have focused on two main possible sources of over-dispersion in spike counts. The first, is an extrinsic source that accounts for sensory noise, conveyed along the thalamo-cortical pathway. The second is a source of noise that emerges from convergent inputs to a population, including feedforward, lateral and feedback. We have shown that the convergence of inputs between two data spaces may both result in an over-dispersion. In the next section, we will focus on quantitatively testing evidence for this over-dispersion.

2 The Fano Gamma test: How to quantify over-dispersion evidence?

The Fano factor test is misleading Most studies concerned with the characterization of spike count variability using different methods and formal tests start with the Fano factor (FF). However, most of the methodological studies measuring trial-to-trial variability in firing rates [START_REF] Koch | Predicting the visual world: silence is golden[END_REF][START_REF] Teich | Fractal character of the neural spike train in the visual system of the cat[END_REF][START_REF] Dayan | Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems 1st[END_REF] do not take the estimation uncertainty introduced by the limited number of trials into account. Each trial may be considered as a realization of the random process and when accumulating evidence from N trials, the average and the variance converge with an expected error of 1/ √ N . Instead, the computed mean and variance are often treated as real features of the distribution whereas they actually correspond to estimations given the (necessarily) limited set of measurements. Thus, considering only the FF value may induce errors in conclusions regarding the assumption of probability distribution. In particular, when the observed FF deviates from 1, is this the result of the limited number of trials or due to an over-dispersion? Or when a value is measured close to 1, can it be safely attributed to a Poisson model?

The chi-square test is problematic To evaluate how likely the observed data was generated by a Poisson model, the most common approach is to use a chi-square goodness of fit test (χ 2 ). This statistical tool can be used to estimate how closely an expected distribution matches an observed distribution of a discrete quantitative variable (having only finite possible values). However, it has been shown that small expected frequencies may result in a loss of power, that is, a tendency to not reject a false null hypothesis (or to reject a true null hypothesis) [START_REF] Cochran | Some methods for strengthening the common χ 2 tests[END_REF][START_REF] Yates | Practice Of Statistics Advanced Placement 1st[END_REF][START_REF] Campbell | Chi squared and Fisher Irwin tests of two by two tables with small sample recommendations[END_REF]. This is a potential problem when we deal with small sample sizes and its outcome also depends on how data is pooled into categories. Indeed, it often makes sense to pool the small frequencies together. Moreover, as this test is one-tailed , it does not distinguish evidence for over-dispersion from evidence for under-dispersion.

The Fano Gamma test To better quantify the over-dispersion while accounting for sample size and without dealing with pooling problems, we propose a statistical test based on an empirical sampling distribution of the Fano factor under the Negative-Binomial assumption. Here, we opt for an analytic procedure that we call Fano Gamma test (F G). It allows for a simple computation (or tabulation) of probability bounds for specific distributions hypotheses. That is not the case for other existing procedures [START_REF] Nawrot | Measurement of variability dynamics in cortical spike trains[END_REF][START_REF] Amarasingham | Spike count reliability and the Poisson hypothesis[END_REF][START_REF] Jackson | Network dynamics of hippocampal cell-assemblies resemble multiple spatial maps within single tasks[END_REF]. This procedure characterizes the probability distribution of the estimated FF. Thus, it allows one to compute FF bounds with respect to p-value precision limits. Starting from previous analyses [START_REF] Hoel | On Indices of Dispersion[END_REF][START_REF] Kathirgamatamby | Note on the Poisson Index of Dispersion[END_REF][START_REF] Eden | Drawing inferences from Fano factor calculations[END_REF] about the distribution of the Fisher dispersion index, we propose here that, for a Negative-Binomial distribution, the FF asymptotically follows a Gamma distribution with a shape parameter n-1 2 and a scale parameter

2( µ φ +1) (n-1) : F F N egative-Binomial ∼ Γ( n -1 2 , 2( µ φ + 1) (n -1) ) ( 6
)
where n is the number of spike count observations. We tested the accuracy of this approximation through simulation by comparing it to the empirical estimates of the true distribution of the Fano factor for any (µ, φ, n) triplet. We found that the p-value upper bounds resulting from the approximative formula are slightly over-conservative. Thus, it is worthwhile mentioning that using it to reject the NBM hypothesis is based on a strong assumption (i.e. a rejected realization could be a good one because it's p-value is between the approximation and the empirical upper bound) but we do not detail this analysis here. For large values of φ, the Negative-Binomial distribution approaches a Poisson distribution and the previous formula gives:

F F P oisson ∼ Γ( n -1 2 , 2 n -1 ) (7)
Eden & Kramer [START_REF] Eden | Drawing inferences from Fano factor calculations[END_REF] studied the FF distribution under the Poisson distribution hypothesis, and have shown that the convergence to this asymptotic distribution is fast. Thus, this analysis provides a formal test for both hypotheses (PM and NBM) given the observed variability in the spiking. We first compared the power of this test to the chi-square test for the general NBM hypothesis and for the specific PM hypothesis on surrogate data of spike counts. Figure 4 shows that the Fano Gamma test did better than a chi-square goodness of fit test at a significance level of 0.025. This latter test has a larger probability to make a type II error for both distributions (dashed red lines in Figure 4.C and dashed blue lines in Figure 4.D).

Over-dispersion evidence from neurophysiological data We applied these tests to evaluate evidence against2 the Poisson model on three data sets.

The first data set results from recording signals from mouse LGN neurons (72 cells) in response to moving gratings. The stimulus was presented several times in 12 different directions. The second data set comes from extra-cellular recordings in area V1 (67 cells) of two awake macaque monkeys in response to an oriented moving bar. The task was repeated several times for 12 different directions (orthogonal to the bar's orientation). The third data set results from recording in area MT (40 cells) of three anaesthetized macaque monkeys in response to moving gratings. The stimulus was presented several times in 16 different directions. See the methods section for more details on the data collection.

We show the results in Figure 5. For the three data sets, we found that the Fano Gamma test gives lower evidence against the PM and the NBM compared to the FF and χ 2 tests (as expected given our discussion above). The Fano Gamma test shows, as do the other tests, that the evidence of over-dispersion is only significant in MT (with a rate of 64%) among cell/condition pairs, compared to LGN and V1 (with rates of 19% and 26%, respectively). Similarly, it shows that evidence against the NBM is significantly lower in the MT dataset (15%) compared to LGN and V1 data sets (52% and 41%, respectively). For these particular datasets where the number of trials is limited (mostly less than 20 trials per condition), we have shown that in the LGN and V1 data sets, the high values of Fano factor observed could be simply related to the sampling of the stochastic process itself. Indeed, looking only at the FF values induces a false belief that these datasets exhibit a significant degree of over-dispersion, while the FG test reliably rejects this hypothesis. The χ 2 test shows slightly more evidence against the Poisson model than the FG test which could be explained by categorization issues. It shows also very strong evidence against the NBM model known to be suitable for over dispersed data. Finally, we have shown that the situation is different for the MT data set where we could indeed show evidence of a higher variability than expected from a Poisson model and less evidence against a Negative Binomial model. It is worth mentioning that our statistical test (as validated on surrogate data) allows us to make this conclusion, while other tests fail to do so. Since we mostly observe over-dispersion in the MT data set, we propose to compare the performance of an NBM-based to a PM-based direction decoder applied to this data set.

NBM-based direction decoder: what if the over-dispersion is tuned ?

To evaluate the gain in decoding when accounting for the over-dispersion using NBM compared to PM, we extended the classical probabilistic decoding approach proposed by Jazayeri & Movshon [START_REF] Jazayeri | Optimal representation of sensory information by neural populations[END_REF] and tested the population decoder on the MT data set (presenting the highest over-dispersion probability). Based on a Poisson model of trial-to-trial variability applied to single neurons and a model of tuning for a population Each subplot shows the measures of variability for the different cells and conditions as points on a (mean, variance) plot. The straight line corresponds to FF=1 (Poisson), while dots correspond to observations. Blue dots correspond to observations with no evidence against Poisson, red dots correspond to no evidence against the Negative-Binomial model (NBM) and black dots correspond to strong evidence against the Poisson model (PM) on the top and against NBM on the bottom based on the empiric distribution of the Fano factor (FG test). The results show that there is more evidence of over-dispersion in MT and that the proportion of cells showing over-dispersion seems to increase along the visual hierarchy. (B) The percentage of observations with evidence against PM and NBM are given with respect to three criteria: the Fano factor test discards the observations with a Fano factor value higher or lower than its formal expected value (equal to one for the Poisson model), the χ 2 test and the F G test with a significance level of 0.025. The bar plots correspond to the percentage of observations with evidence against the NBM (left) and the PM (right). The tests were applied to 3 data sets (from the left to the right, respectively): LGN cells (mouse), V1 cells (awake macaque monkey) and MT cells (anaesthetized macaque monkeys). The results show that the Fano factor values are misleading when the number of trials is not taken into account. Looking only at the FF gives strong evidence against PM and NBM that is not confirmed by the other tests. The χ 2 test shows slightly more evidence against the Poisson model than the FG test which could be explained by categorization issues. It shows also very strong evidence against the NBM model known to be suitable for over dispersed data. In summary, the FG test gives more coherent results: weak evidence of over-dispersion in LGN and V1 data sets and low evidence against the NBM in MT data set. coded variable (direction θ). This decoding approach allows us to infer the stimulus (the hidden variable θ) given a particular neuronal response (an observed spike raster Y).

The decoding algorithm consists of maximizing the posterior probability P (θ|Y ) function of the estimated direction θ given a distribution hypothesis. Bayes' rule gives: P (θ|Y ) = P (Y |θ)P (θ)

P (Y )
, where the evidence term P (Y ) is a normalization term independent of θ. So, if we assume that there is no prior knowledge on θ (P (θ) constant), maximizing the posterior under the PM hypothesis is equivalent to maximizing the following loglikelihood function LP (θ):

LP (θ) = N i=1 k i log[f i (θ)] - N i=1 f i (θ) - N i=1 log[k i !], (8)
where f i is a function of the stimulus parameter θ and corresponds to the mean response of each cell i, k i is the number of trials and N is the number of cells. Similarly, under the NBM hypothesis, we maximize the following log-likelihood function LN B(θ) in order to estimate θ, where p i (θ) = φi(θ) φi(θ)+fi(θ) and φ i corresponds to the inverse-dispersion parameter of each cell i as a function of the stimulus parameter θ:

LN B(θ) = N i=1 log[ Γ(k i + φ i (θ)) Γ(φ i (θ))k i ! ] + N i=1 φ i (θ) log[p i (θ)] + N i=1 k i log[1 -p i (θ)], (9)
Given that objective, let us first define the tuning functions of the cells. These correspond to the mean spike count (the mean parameter of the Negative-Binomial model) with respect to direction. Thus, we defined for the MT data set a generative tuning function as von Mises functions peaking on the preferred direction. The parameters corresponding to the scaling, maximum and minimum values (κ o , κ d , R 0 , R 0 +R m ) are free to vary. The main tuning function is given by:

f M T (θ) = R 0 + R m • e κ d •cos(θ-θ0) -e -κ d e κ d -e -κ d (10)
Figure 6.B shows the estimated tunings of f for 16 different cells. Then, we should consider the estimation of inverse-dispersion parameter φ. We tested two hypotheses. First, having no prior on the tuning of φ as a function of the stimuli direction, we supposed that this parameter is stimulus independent (constant chosen as the average value of its estimates over the different conditions) and varying between cells. Second, we tested the hypothesis of a stimulus dependent inversedispersion parameter. We considered a von Mises tuning of φ centered on the preferred direction for each cell. Figure 6.D shows the estimated tunings of φ for 16 different cells. Please note that our aim is not to have an optimal fit but to show that a simple tuning hypothesis of φ can improve the quality of decoding, even if the estimation of φ values is difficult (see discussion).

Finally, to perform population decoding, we used a leave-one-out (LOO) cross validation method. In this method, we randomly pick a single trial from the data set for each neuron and each condition. Then, we apply the decoding approach previously described: we infer an estimate of the encoded direction (stimulus direction) that maximizes the log-likelihood functions based on the Poisson model (eq. 8) or the Negative-Binomial model (eq. 9) using the spike count of each cell given by the picked trial (as k i ) and the estimated values of mean and inverse-dispersion using the remaining trials(as φ i and f i ). Then, we compare the estimated direction angle to the stimulus direction for the PM model and the NBM model under tuned and un-tuned dispersion hypotheses.

The resulting decoding error values for both models applied to MT data show that the precision of decoding of the PM and the NBM under the non-tuned dispersion hypothesis are very similar (Figure 7.A). However, we found that there is a significant gain in decoding when using the NBM under the tuned dispersion hypothesis (Figure 7.C). To test if this is related to some data or estimation bias, we performed the same test on surrogate data. We generated Negative-Binomial data using the estimated parameter and applied the two decoding schemes. Similarly, the results of decoding on surrogate data show that the precision of decoding is not improved when using non-tuned dispersion (Figure 7.B). However, when considering additive knowledge about over-dispersion tuning, there is a clear gain in accuracy and precision (mean and variance Figure 7.D) which seems to follow decoding performance in real data.

Moreover, comparing the goodness of fit of these models to MT data favors strongly the tuned NBM. Indeed, we used the Bayesian information criterion (BIC) that penalizes models with additional parameters [START_REF] Box | Time Analysis, Forecasting and Control[END_REF] to determine which model best explained our data. With the BIC measure, 37 (respectively 28) over In summary, we have shown two main results in this section. First, the Poisson model performs as well as the Negative-Binomial model in the decoding task most of the time but not for good reasons as the data set is not well described by a Poisson variability model. Second, the Negative-Binomial model with tuned over-dispersion performs better and could be a generic model of neural computation.

Discussion

In this paper, we explored over-dispersion in neural spike counts using a doubly stochastic model, namely the Negative-Binomial model (NBM), that could disentangle different sources of over-dispersion. We have first explained that some commonly used procedures to quantify over-dispersion (the Fano factor and chisquare goodness of fit tests) could be problematic. Then, we have shown that the alternative analytical procedure that we propose (the Fano Gamma test) better characterizes over-dispersion in three different data sets (different animals, tasks, conditions and areas). Using the Fano Gamma test we show that there is weak evidence of over-dispersion in LGN and V1 data sets and low evidence against the NBM in MT data set, especially for high mean spiking values, a trend expected from [START_REF] Goris | Partitioning neuronal variability[END_REF].

In the last section, we compared the performances of the traditional Poisson model to the more recent Negative-Binomial model on a population decoding algorithm (direction decoding from spike recording in MT cells of an anesthetized macaque monkey). Under the Poisson hypothesis, a linear decoding model is optimal, as proposed by Jazayeri & Movshon [START_REF] Jazayeri | Optimal representation of sensory information by neural populations[END_REF]. We show here that this decoder is still efficient (mean decoding error < 10 degrees) even though there is strong evidence against the Poisson hypothesis. This is partly due to the tuning profile of the likelihood functions: the direction corresponding to the maximum likelihood will converge toward the stimulus' direction even if the noise model is not optimal.

Intuitively, one could expect that using the compound model that accounts for the observed overdispersion, one should obtain better decoding results. We first considered the hypothesis of non-tuned cells' dispersion (i.e. each cell dispersion parameter is constant). The Poisson model and The Negative-Binomial model with non-tuned dispersion show similar direction decoding results, not only on real data, but also on surrogate data. However, note that, similarly to what was reported before [START_REF] Goris | Partitioning neuronal variability[END_REF], we found that the likelihood values of the estimated directions are always higher for the Negative-Binomial based decoding compared to the Poisson based decoding (not shown), but with similar bandwidths. This leads to similar decoding accuracies, despite the fact that the NBM fits better to the data. Indeed, the spike counts corresponding to the response of the cell to its preferred direction are more likely to be generated by an NBM with higher inverse-dispersion values, i.e. Poisson-like. This could explain why PM performance could be comparable to NBM performance in decoding.

In addition, the estimation of the inverse-dispersion parameter is known to be difficult. Indeed, the estimation of this parameter has been widely studied in different fields [START_REF] Mccullagh | Generalized linear models[END_REF][START_REF] Gourieroux | Pseudo maximum likelihood methods: Applications to Poisson models[END_REF][START_REF] Lawless | Negative binomial and mixed Poisson regression[END_REF][START_REF] Dai | Maximum likelihood estimate for the dispersion parameter of the negative binomial distribution[END_REF] and several estimation methods were provided such as the maximum likelihood method [START_REF] Piegorsch | Maximum likelihood estimation for the negative binomial dispersion parameter[END_REF], or the method of moment and maximum extended quasi-likelihood method [START_REF] Clark | Estimation of the negative binomial parameter κ by maximum quasilikelihood[END_REF]. A wide range of studies have shown that the estimation of this parameter is very challenging and can be significantly biased, especially when dealing with small data sets. It is worthwhile to point out this problem here: Since neurophysiological experiments are costly and time consuming, we often deal with small data sets. Thus, having more knowledge and control about the dispersion tuning is crucial.

Therefore, we tested a simple possible form of dispersion tuning, that is, the tuning of the inversedispersion having similar preferred direction as the tuning of the spike count (using a von Mises function). We found that under the hypothesis of a stimulus-dependent inverse-dispersion, the NBM performs qualitatively better.

Moreover, such a tuning function could explain complex Fano factor tuning profiles that were recently reported [START_REF] Ponce-Alvarez | Stimulus-dependent variability and noise correlations in cortical MT neurons[END_REF]. Our results lead us to propose that the reported trial-by-trial variability among MT neurons -showing a directional tuning that is not trivially explained by firing rate variations alone -may simply be explained by a bell-shaped tuning of the over-dispersion (see Figure 8). However, contrary to what they expect, we did not observe a decrease of over-dispersion with respect to the spontaneous level. Indeed, applying the Fano-Gamma test to MT spontaneous spiking data (results not shown) yields lower evidence against a PM (47.5%) compared to that for evoked spiking (64.1%, Figure 5.D) but similar evidence against NBM (15%) compared to evoked spiking (16.3%, Figure 5.D). This is more in line with an increase of The resulting FF profile could not be explained only by the mean response tuning as in a Poisson model. under-dispersion, but which still in accordance with the known fact that FF is decreased when a stimulus is applied.

Where could such a tuning of the over-dispersion come from? One first hypothesis could be that over-dispersion results from the sensory noise. Then, as suggested by Goris, Movshon & Simoncelli [START_REF] Goris | Partitioning neuronal variability[END_REF], over-dispersion could also reflect an excitability state of the cell. However, both of these hypotheses cannot account for stimulus dependent variability. Alternatively, we suppose that this property is better apprehended at the population scale. Going back to the encoding paradigm described in the first section (Figure 3), an alternative hypothesis would be that the over-dispersion arise from a stimulus-dependent degree of network convergence (within MT or from V1). Such a scheme could result from a center-surround connectivity profile in the direction domain: more excitatory input for a stimulus close to the preferred direction of the recipient cell.

This study calls for an extension over a wider range of experimental observations. First, we would like to extend this work to more datasets with more trials to explore if our predictions are still valid, both for teasing out sources of variability but also to evaluate population decoding. It would be mostly interesting to compare this extra variability along the sensory pathways to associative areas to confirm the tendency that we observed with a gradual increase of over-dispersion from the lower to associative areas. A second important perspective is to consider the role of experimental conditions. Indeed, it is likely that results would greatly vary depending on the arousal state of the animal, the nature of the stimulation or the behavioral task. In particular, behaviorally more relevant stimulations (such as the model-based synthesis of textures [START_REF] Sanz-Leon | Motion Clouds: model-based stimulus synthesis of natural-like random textures for the study of motion perception[END_REF] may lead to a more precise response [START_REF] Baudot | Animation of natural scene by virtual eye-movements evokes high precision and low noise in V1 neurons[END_REF].

Materials and Methods

In our analyses, we used three different data sets from neurophysiological recordings.

LGN data set This data set [START_REF] Scholl | Spiking responses of LGN neurons in the mouse to drifting gratings CRCNS[END_REF] contains the spiking responses of [START_REF] Box | Time Analysis, Forecasting and Control[END_REF] LGN neurons in the anesthetized mouse to drifting grating (square wave (n=24) and sinusoidal (n=48)). Many of these neurons are already contained in the data set published in [START_REF] Scholl | Emergence of orientation selectivity in the Mammalian visual pathway[END_REF]. Contrarily to the two other data sets, it is an openly available data set.

V1 data set Experiments were conducted at the INT on two adult male rhesus macaque monkeys (Macaca mulatta). Experimental protocols have been approved by the Marseille Ethical Committee in Neuroscience (approval #A10/01/13, official national registration #71-French Ministry of Research). All procedures complied with the French and European regulations for animal research, as well as the guidelines from the

  the Odds of Inherent versus Observed Over-dispersion in Neural Spike Counts}, author = {Taouali, Wahiba and Benvenuti, Giacomo and Wallisch, Pascal and Chavane, Frederic and Perrinet, Laurent U.}, journal = {Journal of Neurophysiology}, year = {2015}, doi = {10.1152/jn.00194.2015}, issn = {1522-1598}, volume = {115}, number = {1}, pages = {434-444}, url = {http://dx.doi.org/10.1152/jn.00194.2015}, url = {http://invibe.net/LaurentPerrinet/Publications/Taouali15}, pmid = {26445864}, publisher = {American Physiological Society}, }

Figure 1 :

 1 Figure 1: Trial-to-trial variability of neural responses. (Left) We show the spiking activity of a representative cell of the area MT of a macaque monkey (see Methods section) in response to oriented drifting gratings. Each subplot corresponds to a raster plot of the responses to 16 equally distributed motion directions. There is a high trial-to-trial variability in the spiking responses. (Right) To quantify this variability, we plot the directional tuning curves in response to the same 16 directions (mean spike count, black lines) along with the range of ± the standard deviation (red lines) for a population of 16 cells from the same area. The subplot with a star corresponds to the representative cell on the left. The area corresponding to the mean ± the standard deviation expected under Poisson distribution hypothesis (σ = √ µ) is plotted in blue. In this experiment, most of the cells show a higher variability than expected from the Poisson model.

Figure 2 :

 2 Figure2: Over-dispersion may result from extrinsic sensory noise: NBM as a mixture of Poisson and LogNormal distributions. We show here a simple linear / non-linear spiking model with three processing steps: (A) The correlation between an noisy Gabor image (input) and a population of Gabor filters at different orientations models the linear processing in the receptive field (RF) of a visual neuron. The green, black and red lines correspond respectively to one trial (simulation), the average and the average ± the standard deviation (within 1000 trials). (B) This input is transformed using a static, exponential non-linearity multiplied by a scaling factor corresponding to the maximum spike count. The input drive resulting from a Gaussian input noise and the non-linearity follows a Log-normal distribution that is closely fitted by a Gamma distribution. (C) This drive is finally transformed into spike counts using a Poisson point process with the mean as the input drive. The resulting tuning curve (mean spike count) is plotted in black along with the range of ± the standard deviation (red line). The area corresponding to the mean ± the standard deviation expected given the Poisson model (σ = √ µ) is plotted in blue. The resulting spike counts variance is higher than what is expected by a Poisson model, as can be observed in neurophysiological data (see Figure1-Right).

Figure 3 :

 3 Figure 3: Over-dispersion may result from redundancy within a neural population.NBM as mixture of Poisson and Exponential distributions. Let us consider, within a neural population, a subset of φ cells with a similar selectivity and thus driven by a similar correlation measure λ= C(θ).This measures the correlation between the input Gabor and the Gabor representing the orientation θ characteristic of their receptive field (RF). Let us further assume that the filter cells exhibit an exponential distribution of spike counts E( 1 λ ) with mean λ, thus maximizing the entropy of the input drive. Then, the pooling of the resulting spike counts results in an input drive G(φ, λ) equal to the sum of exponential variables having the same mean equal to λ. It can be shown that the Poisson spiking process driven by this input drive is equivalently described by a Negative-Binomial probability distribution function with parameters λ (mean) and φ (inverse-dispersion) (see equation 2). Such a distribution describes observed trial-to-trial variability well and in this model, the parameter of dispersion is directly correlated with the redundancy, i.e., the convergence of similar input cells.

Figure 4 :Figure 5 :

 45 Figure 4: Evaluating the power of chi-square test (χ 2 ) test versus Fano Gamma test (F G) when testing evidence against Poisson (PM) model (in blue) and Negative-Binomial (NBM) model (in red). (Top) Two examples of spike count histograms resulting from artificially generated data (surrogate data) using the NBM (A) and the PM (B). These examples were chosen such that, when using the χ 2 test, the model used to synthesize the data is rejected. For each histogram, we plotted in blue the expected distribution from a PM and in red the expected distribution from a NBM. For these examples, the sample size is 50. (Bottom) χ 2 test (dashed lines) and F G test (solid lines) were applied to Negative-Binomial (C) and Poisson (D) surrogate data for different sample sizes with a significance level of 0.025. The mean is fixed to 10 spikes, the dispersion parameter is equal to 1 spike and the test was repeated 1000 times for each sample size. (E) and (F) show the results in (C) and (D), respectively, corresponding to small sample sizes. The rejection rate is higher for χ 2 test when the null hypothesis is the correct distribution (the model generating the data) and lower for the incorrect distribution than for the proposed Fano Gamma test. The results thus show this test is more powerful than the classically used χ 2 test.

44 *Figure 6 :

 446 Figure 6: Mean and dispersion tunings (Left) The box plots represent the distribution of cells' mean spike count (A) and the distribution of their inverse-dispersion values (C) function of the different directions centered around the preferred direction of each cell. The "whiskers" for each direction indicate the 25thand the 75th percentiles. The colored dashes (blue for mean and red for inverse dispersion) represent the median values. The inverse-dispersion parameter (C) shows a tuning less regular than the mean (A) but that could, as a first approximation, be fitted by a bell-shaped function. (Right) Tuning curves of the mean spike count (B, in blue) and the inverse-dispersion (D, in red) for 16 cells resulting from, respectively, the fit of a von Mises function (equivalent to a circular Gaussian) to each cell estimated parameters from their repeated responses to the 16 directions used for training (gray points). Each subplot corresponds to a cell. The tuning functions of the inverse-dispersion are supposed to be centered on the cell's preferred direction deduced from the mean fit.

Figure 7 :

 7 Figure 7: Population decoding : Poisson model versus Negative-Binomial model. Variation of decoding error values using the Negative-Binomial model (NBM) function compared to that using the Poisson model (PM) both applied to real MT data (A and C), then to surrogate data generated using the parameters inferred from the real data (B and D). We tested two alternatives of the NBM : one with an overdispersion parameter stimulus-independent (non-tuned dispersion) and the other fitting this parameter by a von Mises function of the direction centered on the preferred direction (tuned dispersion). Points correspond to the median error values and bars to 25%-75% confidence interval within 100 trials (horizontal for PM and vertical for NBM). The Poisson model and The Negative-Binomial model with non-tuned dispersion show similar results. However, the Negative-Binomial model with stimulus-dependent dispersion parameters shows better results compared to the two others. This suggests that tuned over-dispersion could play a role in the neural code.

Figure 8 :

 8 Figure 8: Complex tuning of Fano factor (FF) values for different motion directions. (A) An example of a mean response tuning curve. (B) An example of an inverse-dispersion tuning curve. (A) and (B) were analytically generated by Von Mises (circular Gaussian) functions centered on a given preferred direction (PD). (C) Let us suppose that the spiking count of a cell is driven from a Negative-Binomial model characterized by mean values given in (A) and inverse dispersion values given in (B). The expected resulting Fano factor values are given by the relationship equation (F F = 1 + f φ ). For example, the mean value (7.14) and the inverse dispersion value (21.47) corresponding to the PD result in an FF value of 1.33. The resulting FF profile could not be explained only by the mean response tuning as in a Poisson model.

The inverse of the parameter φ, usually denoted by α, has been often called in the literature as the "dispersion parameter". It has been occasionally named the "over-dispersion parameter" because it does not account for FF values smaller than one.

Indeed, an hypothesis test is a statistical analogy to proof by contradiction. Thus, we used the "evidence against" term that is the result of the statistical test and that means the rejection plausibility. Importantly, the evidence against a claim H0 does not imply the evidence (plausibility) of the opposite claim H1, i.e., this is not the same as a likelihood based test that allow explicitly to quantify the weight of evidence for one hypothesis.

Acknowledgments

WT, FC and LP were supported by EC FP7-269921, "BrainScaleS", GB by FACETS ITN project (EU funding, grant number 237955), a 'Marie-Curie Initial Training Network' and PW was supported by F32-EY019833 from the NIH. We wish to thank Philippe Foundation for a travel grant support to F.C. The authors are indebted to T. Movshon and his laboratory for giving the opportunity to F.C. to perform MT experiments at CNS/NYU.

Society for Neuroscience. Surgical preparation : Monkeys were chronically implanted with a head-holder and a recording chamber located above the V1 and V2 cortical areas. In a second surgery, a search coil was inserted below the ocular sclera to record eye movements [START_REF] Robinson | A method of measuring eye movemnent using a scieral search coil in a magnetic field[END_REF]. Behavioral task and training : Monkeys were trained to fixate -within a window of 1 to 2 degrees of diameter -on a red target presented at the center of the screen during the entire duration of the trial (1.5 s). The eye position was monitored by the scleral search coil technique [START_REF] Robinson | A method of measuring eye movemnent using a scieral search coil in a magnetic field[END_REF][START_REF] Collewijn | Precise recording of human eye movements[END_REF] and the animal's behavior controlled using the REX package [START_REF] Hays | A UNIX-based multiple process system for real-time data acquisition and control[END_REF]. Unit Recording : The two monkeys were chronically implanted with a head-holder and a recording chamber located above cortical area V1. A computer controlled microdrive (MT-EPS, Alpha Omega, Israel) was used to trans-durally insert a micro-electrode (FHC, 0.5-1.2 MΩ at 1 KHz) in the right hemisphere. Spikes were sorted online using a template matching algorithm (MSD, Alpha Omega). Visual Stimulation : Visual stimulation protocols have been produced using in-house software (developed by Gérard Sadoc, Acquis1-Elphy, Biologic CNRS-UNIC/ANVAR). Stimuli were back-projected on a translucent screen covering 80 • ×80 • of the monkey's visual field, at a distance of 1 m, using a retro-projector (resolution: 1280×1024 pixels at 60 Hz). The mean luminance of the motion stimulus was 22.2 cd/m 2 and the background was kept constant to about 2.24 cd/m 2 . The display was gamma calibrated by means of a lookup table. Once a cell is isolated, a sparse noise (SN) stimulus was first used to quantitatively map out the RF (a grid of 10x10 squares with a side length of 0.6 degrees, 10.9 cd/m 2 darker or brighter than the background of 11.1 cd/m 2 , ∼15 trials). Once the RF was properly located and estimated, a direction tuning paradigm was launched. In this paradigm, a moving bar (0.5×4 • ) translated across the receptive field (the displacement covers 3 • on either side of the RFc) in 12 possible directions (spaced by 30 degrees of polar angle, at a speed of 6.6 • /s).

MT data set Experiments were conducted at NYU-CNS on 3 anesthetized, paralyzed, adult macaque monkeys (Macaca nemestrina) of either sex. All procedures were conducted in compliance with the National Institute of Health Guide for the Care and Use of Laboratory Animals, and with the approval of the New York University Animal Welfare Committee. Surgical preparation : The standard procedures for the surgical preparation of animals and single-unit recordings have been reported in detail previously [START_REF] Cavanaugh | Nature and interaction of signals from the receptive field center and surround in macaque V1 neurons[END_REF]. Briefly, experiments typically lasted 5 to 6 days, during which we maintained anesthesia with infusion of sufentanil citrate (6-30 µg kg -1 h -1 ) and paralysis with infusion of vecuronium bromide (Norcuron; 0.1 mg kg -1 h -1 ) in isotonic dextrose-normosol solution. We monitored vital signs (heart rate, lung pressure, end-tidal pCO2, EEG, body temperature, urine flow, and osmolarity) and maintained them within appropriate physiological ranges. Pupils were dilated with topical atropine. The eyes were protected with gas-permeable contact lenses, and refracted with supplementary lenses chosen through direct ophthalmoscopy. At the conclusion of data collection, the animal was sacrificed with an overdose of sodium pentobarbital. Unit Recording : Extracellular recordings were made with quartz-platinum-tungsten microelectrodes (Thomas Recording), advanced mechanically into the brain through a craniotomy and small durotomy. To record from MT we passed microelectrodes through a small durotomy centered roughly 16 mm lateral to the midline and 3 mm posterior to the lip of the lunate sulcus at an angle of 20 degrees from horizontal in a ventro-anterior direction. Area MT was identified by the brisk direction-selective responses of isolated neurons. We made recordings from every single unit with a spike waveform that rose sufficiently above noise to be isolated. Stimuli were presented in random order. Visual Stimulation : We presented visual stimuli on a gamma-corrected CRT monitor (Eizo T966; mean luminance: 33 cd/m 2 ) at a resolution of 1.280×960 with a refresh rate of 120 Hz. Stimuli were presented using Expo software (http://corevision.cns.nyu.edu) on an Apple Macintosh computer. For each isolated unit, we first determined its ocular dominance and occluded the non-preferred eye. We presented circularly windowed sinusoidal grating stimuli (16 directions around the clock, 22.5 degrees apart) to map each cell's receptive field, determined its preferred size and speed, and then measured selectivity for orientation or spatial frequency.

The time window considered to calculate the spike counts was the duration of stimulation for the grating stimuli (LGN and MT data sets) and 400 ms covering approximatively the time when the bar stimulus crossed the cell's receptive field for V1 data set. All data were analyzed using Python and Matlab (The MathWorks Inc., Natick, MA, USA).