Coupled variability in primary sensory areas and the hippocampus during spontaneous activity

The cerebral cortex is an anatomically divided and functionally specialized structure. It includes distinct areas, which work on different states over time. The structural features of spiking activity in sensory cortices have been characterized during spontaneous and evoked activity. However, the coordination among cortical and sub-cortical neurons during spontaneous activity across different states remains poorly characterized. We addressed this issue by studying the temporal coupling of spiking variability recorded from primary sensory cortices and hippocampus of anesthetized or freely behaving rats. During spontaneous activity, spiking variability was highly correlated across primary cortical sensory areas at both small and large spatial scales, whereas the cortico-hippocampal correlation was modest. This general pattern of spiking variability was observed under urethane anesthesia, as well as during waking, slow-wave sleep and rapid-eye-movement sleep, and was unchanged by novel stimulation. These results support the notion that primary sensory areas are strongly coupled during spontaneous activity.

Project, which enabled online view of the electrophysiological recordings. The off-line raw extracellular recordings were stored (locally and at http://osf.io) for later processing, including automatic and manual spike sorting.

Spike sorting
Off-line spike sorting from raw extracellular recordings was performed with the software provided by Klusta-Team 6,7 . Since the shanks were >200um, the activity between channels in different shanks was considered independent for spike sorting purposes, therefore spike sorting was performed in parallel, one per shank. The independence among spike data from different shanks was confirmed by the coherence matrices based on segments of raw data after bandpass filtering 300-500 Hz, based on customized Matlab code from Kenneth Harris. Spike sorting began with spike detection (filtered from 500 Hz-11.5KHz, with a 3rd order Butterworth), followed by the automatic spike sorting using KlustaKwik 7 , which is an "masked" expectation-maximization (EM) algorithm for automatic clustering action potential waveforms. The outcome of the automatic part of the spiking sorting is a set of clusters of waveforms and the correspondent values of quality for each cluster. Using the Klustaviewa 7 the set of clusters was refined by merging and splitting operations, and selection adopted the following criteria: (a) low contamination (< 5%) within the refractory period (2ms); (b) spiking stability across the period of interest; (c) minimum firing rate (0.5Hz). The final outcome of the spike sorting procedure was a set of spike trains, in which the i-th spike train was defined by the timestamps of action potentials assigned to i-th neuron within a given neuronal population. The automatic part of the spike sorting was executed in local computers and also in a computer cluster provided by the Portuguese initiative for scientific grid computing (INCD/LIP).

Experiments in freely behaving animals
Behaviours were continuously recorded under visible or infrared light with a CCD camera connected to a video tracking system (Cineplex) synchronized to the neural recordings. Illumination in the visible spectrum was measured inside the recording chambers, and the absence of visual stimuli in the visible spectrum was verified using a Minipa Electronics luximeter, model MLM -1011. Recordings of large neuronal populations were performed using a 96-channels Plexon recording system (MAP, Plexon Inc, Dallas, TX). The animal's behavioural states were classified based on a previously described method, based on hippocampal and cortical LFPs 8 . Briefly, a two-dimensional state space was defined by two spectral amplitude ratios calculated by dividing energy at selected frequency bands from LFPs simultaneously recorded in the areas of interest. The states were classified based on amplitudes ratios of the LFPs within selected bands.

Data analyses
We used the same approach used to asses the coefficient of variation in a previous study 9 , which we will briefly describe herein. For our purpose, the spike train of a given single-or multiunit is represented by the ordered sequence of time-stamps u i = [t i ] in which were detected spikes from that single-or multi-unit. Given a set of N population spike trains S = u i , the correspondent population activity spike train is represented for u i ∈S u i . Some of a spike calculations has been done only over part of spike trains, therefore it would be convenient define a slice of a spike train u k in a given time interval I as the ordered sequence [u k ] I = [t k ;t k ∈ I] or just u I k . It also convenient to define the time-series based on single spike trains, and the most known is the spike count z k = h(u k , ∆t), which is given by the histogram h over a spike train uk using a bin size t; and similarly would be done over slice of spike trains z I k = h(u I k , ∆t); and even over the population activity in S, z I S = h(u I S , ∆t).

Coefficient of variation
The coefficient of variation, CV, is an example of calculation over slices of spike trains: the CV of the neuronal activity found in neuronal population, S, during a time interval I, in a given time scale, is given as ratio between the standard deviation and mean found in that spike count, such as is defined in Equation (1).
where z I S is the corresponding spike count for u I S using a time scale ∆t. Herein, the default parameters to calculate CVs were ∆t = 50ms, over consecutive non-overlapping 10s-long time periods. Whereas in Fig S2 were used different time-scales, such indicated in its horizontal axis.

Spiking correlations
Usually the spiking data can be modulated by exogenous signals, which would distort the estimation of spiking correlations, usually biased by slow fluctuations in exogenous signals 10 . To avoid this issue and to be able to specify, a priori, the time scale for the spiking correlations, we took advantage of an approach applied in previous studies 9 , which convolved a Mexican-hat kernel with a spike count version of raw data, such as is briefly described in this section.
Given a set of spike trains, S = {u i }, and a corresponding set of spike count time-series Z = {z i } at millisecond resolution (∆t = 1ms). Before to quantify the spiking correlation is demanded use a time-filter in order to specify the time-scale in which the spiking correlation should be estimated by using a Mexican hat filter, which in its turns is derived from the difference between two Gaussians, with zero mean and different standard deviations, s 1 and s 2 , whose define the time-scale of hat filter, h s 1 ,s 2 (t). Therefore, given a time-scale and its corresponding Mexican hat filter, h s 1 ,s 2 (t), the first step in order to quantify the spiking correlation is convolve it with each time-series in Z , and then produce an equivalent set of rate (discrete) functions {n i (t)}, where the i-th rate function is given by the convolution of the filter h(t) with i-th spike count time-series, z i (t), in Z, such as is defined in the Equation (2).
where s 1 and s 2 are the different standard deviations used to build the Mexican-hat filter, with s 1 > s 2 , which define the time-scale to be observed into the neuronal activity data; in the current data analysis we have used s 2 = √ 17s 1 , with the following default values: s 1 = 100ms and s 2 = 400ms.
Given the set of rate functions F = {n i (t)} the ultimate step was calculate the covariance matrix. The covariance between a pair of rate functions, n i and n j is defined by Equation (3).
where L is the (integer) number of milliseconds in each rate functions. Once the kernel used in Equation (2) produces zero mean rate functions, the Equation 3 conveniently becomes a dot product such as in Equation (4).
Given that each rate function is a finite discrete time function, with millisecond resolution, each one of those rate function can be represented as finite sequence, and based on those sequences we defined the rate matrix, B = [n i, j ] N×N , where each n i, j is the instantaneous firing rate of the i-th neuron on the j-th millisecond, N is the number of recorded neurons. Based on rate matrix B, the covariance matrix is given by the Equation (5).
The correlation matrix, C = [r i, j ] N×N can be calculated from the covariance, where r i, j = σ i, j σ i ·σ j , and s k is the standard deviation of k-th rate function, n k , and σ i, j is given the Equation (3). The eigenvectors of the covariance matrix, C, defined in Equation (5), are called the principal components (PC) of this rate data, and the coordinates of a given PC are called its the coefficients or loadings, which define the linear combination of the rate functions in the direction to that PC, for instance the first principal component (PC1) over time is given the Equation (6).
where λ i, j is the j-th loading of the PC1, and n j (t) is the j-th rate function such as has been defined in the Equation (2). In order to avoid sensitivity to heterogeneity in the individual variance usually found in neuronal populations, we used PC based on correlation matrices instead of based on covariance matrices 11 .

Variability around the mean activity
Given that mean population activity in neuronal population S with N units is n(t) = 1 N ∑ N i=1 n i (t), and taking σ 2 n as its variance, from basic statistics:

3/11
where σ 2 i is variance of the i-th rate function, and σ 2 i, j is the covariance between a pair of rate functions n i , n j . Given that the average covariance is σ 2 (7) can be rewritten: Therefore, when N 1, σ 2 n ≈ σ 2 × . Thereby, in a large neuronal network, the size of fluctuation over the mean activity is in the same order of the average covariance in that population. Therefore, just by looking at the level of variability around global fluctuation in local neuronal population is possible estimate the kind of spiking correlation structure can be found in that neuronal substrate, and consequently its dynamics: a synchronized neuronal activity is related to a covariance/correlation structure with positive mean, and thus large fluctuations around the mean activity; whereas desynchronized neuronal activity is related to a near zero covariance/correlation structure, and thus small fluctuations around the mean activity 9, 12-15 .