Functional two-way analysis of variance and bootstrap methods for neural synchrony analysis
- Aldana M González Montoro^{1}Email author,
- Ricardo Cao^{1},
- Nelson Espinosa^{2},
- Javier Cudeiro^{2} and
- Jorge Mariño^{2}
DOI: 10.1186/1471-2202-15-96
© González Montoro et al.; licensee BioMed Central Ltd. 2014
Received: 18 December 2013
Accepted: 31 July 2014
Published: 12 August 2014
Abstract
Background
Pairwise association between neurons is a key feature in understanding neural coding. Statistical neuroscience provides tools to estimate and assess these associations. In the mammalian brain, activating ascending pathways arise from neuronal nuclei located at the brainstem and at the basal forebrain that regulate the transition between sleep and awake neuronal firing modes in extensive regions of the cerebral cortex, including the primary visual cortex, where neurons are known to be selective for the orientation of a given stimulus. In this paper, the estimation of neural synchrony as a function of time is studied in data obtained from anesthetized cats. A functional data analysis of variance model is proposed. Bootstrap statistical tests are introduced in this context; they are useful tools for the study of differences in synchrony strength regarding 1) transition between different states (anesthesia and awake), and 2) affinity given by orientation selectivity.
Results
An analysis of variance model for functional data is proposed for neural synchrony curves, estimated with a cross-correlation based method. Dependence arising from the experimental setting needs to be accounted for. Bootstrap tests allow the identification of differences between experimental conditions (modes of activity) and between pairs of neurons formed by cells with different affinities given by their preferred orientations. In our test case, interactions between experimental conditions and preferred orientations are not statistically significant.
Conclusions
The results reflect the effect of different experimental conditions, as well as the affinity regarding orientation selectivity in neural synchrony and, therefore, in neural coding. A cross-correlation based method is proposed that works well under low firing activity. Functional data statistical tools produce results that are useful in this context. Dependence is shown to be necessary to account for, and bootstrap tests are an appropriate method with which to do so.
Keywords
Cross-correlation analysis Bootstrap Spike-trains Dependence Low firing-rate Functional dataBackground
The nervous system consists of a very large number of neurons—and glial cells—that are connected in complex networks. Neurons convey information by means of electrical pulses, called action potentials or spikes, consisting of a very fast and transient depolarization of their membrane potential. Sequences of spikes, called spike trains, are believed to serve as an information code in the brain. Pairwise synchrony between spike trains is widely studied in neuroscience as a way to understand these codes. Although neuronal interactions are probably not described only by pairwise associations [1], they have been shown to provide important information about neural coding [2–5]. On the other hand, the definition of neural synchrony is a matter of debate and different conceptions of it exist, such as “exact spiking coincidence” or “firing rate association”. There are several existing models for synchrony between simultaneous spike trains, most of them based on cross-correlation analysis. Common methods to measure synchrony are, for example, the cross-correlogram or the joint peristimulus time histogram [6]. Other methods not based on cross-correlation analysis include unitary event analysis [7, 8], conditional synchrony measure [9] and a method based on the distances between the closest spikes [10], among others. Most commonly, association measures are used to test for the presence/absence of synchrony. Several methodologies exist for this aim [11, 12]; however, in this paper we address a different problem. We use a cross-correlation based method, called the Integrated Cross-correlation Synchrony Index (ICCSI), which allows us to obtain a synchrony curve as a function of time and search for differences in the strength of pairwise associations regarding different factors. A similar problem was considered by Faes et al. [9], where these authors compare the synchrony between neurons under experimental conditions with synchrony at baseline.
During deep sleep, neurons in the cerebral cortex present a highly oscillatory global activity that can be observed by means of an electroencephalogram (EEG) or an electrocorticogram (ECoG). Furthermore, in the awake state, these global oscillations disappear giving rise to a less synchronized global activity. The transition from the sleep state to the awake state is regulated by the activating ascending pathways, which are afferent neuronal pathways originating in nuclei located at the basal forebrain and the brainstem [13, 14]. Under general anesthesia, the EEG is characterized by the presence of a number of patterns (spindles, K-complexes, delta waves) that show a progressive increase in low-frequency, high-amplitude activity [15–17]. In this scenario, the transition to the awake-like state can also be reproduced by means of microelectrical stimulation of some activating ascending areas [18]. On the other hand, an important property of neurons in the primary visual cortex (V1) is orientation selectivity; i.e., the neurons respond better (with a higher firing rate) to a specific orientation of the visual stimulus, the so-called preferred orientation [19]. This property is important in neurophysiological studies because it can provide meaningful clues regarding the physiology and microanatomy of the striate cortex.
The aim of the present study is to investigate the differences in synchrony strength regarding different experimental conditions given by anesthesia and awake-like activity, as well as the orientation selectivity of neurons in V1: we study whether similarities among neurons, such as the similarity in preferred orientation, affect the strength of correlated activity. We perform the analysis in a group of neurons recorded from V1 of an anesthetized adult cat. A method is proposed to investigate whether the mode of activity (anesthesia or awake-like), and the affinity in orientation selectivity of a given pair of neurons, have a determinant influence in how neuronal synchronization evolves.
The data are synchrony measurements computed at any time point. Therefore, the data are curves varying continuously over time. This is the reason why, in this paper, the word functional will be used to denote the nature of the data and not to make reference to the neurophysiology of the neurons under study. This term comes from statistics, where “functional data analysis” is a widely developing research field.
A functional two-way analysis of variance is proposed. Functional data analysis tools, based on Cuesta-Albertos and Febrero-Bande [20], are used. The method is adapted to consider the dependence that exists among the data because of the experimental setting. A parametric bootstrap is proposed for hypothesis testing.
Methods
In this section, the data are presented. Also, the synchrony measure used to obtain the functional data is described, as well as the statistical methodology used to cope with the functional analysis of variance (ANOVA) model.
Dataset
Data were recorded from an anesthetized and paralyzed adult cat. A microelectrode array with eight independent movable electrodes was introduced into the primary visual cortex of the animal for neuronal recording. Another two microelectrodes were introduced into the brainstem and basal forebrain for electrical stimulation. These stimulations, which we denote as bs (when the brainstem is stimulated) and bf (when the basal forebrain is stimulated), provoked a change in cortical activity from anesthesia to an awake-like pattern. All experiments followed the guidelines of the International Council for Laboratory Animal Science and the European Union (statute nr 86/809) and the protocols were approved by the University of A Coruña Committee on Animal Care.
At the beginning of each recording, neurons were characterized regarding their preferred orientation. Drifting gratings were used to visually stimulate the cat while the firing activities of a group of neurons were recorded. Each grating corresponded to an angle, which we call orientation, with a specific direction of movement. Orientation (and direction) are continuous variables; however, owing to the nature of the experiments, they will here be considered as discrete. Sixteen possible orientation-direction gratings were used: eight orientations with two possible directions each. For example: a drifting grating at 90° (the lines composing the grating are, therefore, vertical) that moves from right to left is a possible orientation-direction stimulus; another moving from left to right is a different one. Although the use of the two possible directions is also of interest in the study of other properties of V1 neurons (for example, the selectivity to direction), in this work we focus our analysis on the orientation selectivity. Hence, there were eight possible values for orientation: 0°,22°,45°,67.5°,90°,112.5°,135° and 157.5°. So, each recorded neuron was associated with one orientation (the preferred one), corresponding to its maximum firing rate. However, we still need to go one step further, as the objective of the study is to evaluate the effect of orientation selectivity on neural synchrony. To achieve this aim, each pair of neurons is identified with a value of a variable, G, which is defined as the difference between the preferred orientations of the neurons that form the pair. That is, if O_{1} and O_{2} are the preferred orientations of two neurons, we define G=min{|O_{1}−O_{2}|,180°−|O_{1}−O_{2}|}. Given the previous considerations, G can take one of these five possible outcomes: 0°,22.5°, 45°, 67.5° and 90°.
Throughout the paper, we will denote the number of neurons in a simultaneously recorded group by n and the number of possible pairs by $r=\frac{n(n-1)}{2}$. The number of experimental conditions will be denoted by K, the number of trials in each of these conditions will be L and N will denote the total sample size: N=K r L. For our data, n=8, r=28, K=2, L=4 and N=224. The firing rates of the test neurons varied from less than 1 Hz to around 7 Hz; however the most common values ranged from 1 to 3 Hz, denoting typically low firing activity.
Synchrony measure
To measure synchrony, we used a cross-correlation based measure defined as the area under the cross-correlation function in a neighborhood around zero. Data were recorded under spontaneous activity, which is characterized by its very low firing activity. Exact spiking coincidences hardly ever occur, although the global activity of the brain is highly synchronized. By considering the area under the normalized cross-correlation function in an interval around zero we allow for a more relaxed definition of synchronous spike trains than that of just exactly coincident events.
Integrated cross-correlation synchrony index (ICCSI)
represents the probability density function of the time from an event in train _{1} to an event randomly chosen in train _{2}[21].
where δ<1 is an arbitrary small number chosen by the researcher. In this way, we allow synchrony to be based on delayed firing and not only on simultaneous firing. Integrating $\widehat{\mathcal{T}}({\mathcal{X}}_{1},{\mathcal{X}}_{2};\tau )$ in a neighborhood of zero we account for spikes that occur closely in time, though not exactly at the same time.
To study the evolution of synchrony in time, sliding windows can be used. That is, if the total observational time of the spike trains is [0,T], synchrony at time t∈[w,T−w] can be estimated by computing Y(_{1},_{2}) in the time window (t−w,t+w]. Using a time grid, 0<t_{1}<…<t_{ M }<T, and computing the ICCSI at each time point of the grid, synchrony becomes a function of time: Y(_{1},_{2};t). Details on functional data analysis are not given here, though we outline the basic notions that are necessary to understand our analysis. For details and theory on this subject, we refer the reader to, for example, the books [22–24].
We search for population differences in the dynamics of the awake-like period induced by each stimulus, taking into account the possible effect of the other factor: difference between orientation selectivity. This problem can be dealt with as a functional two-way ANOVA with a two-level factor: stimulus and a five-level factor: G.
Functional ANOVA
As already mentioned, the aim of this work is to search for differences in synchrony dynamics relative to two factors: stimulus and difference in orientation selectivity. As the difference in orientation selectivity can take only values in a given finite set of values, the problem is one of a two-way analysis of variance with fixed effects in which the response variable is functional. In the following subsection, we outline the methods.
The random projection method for the ANOVA model
The random projection approach in the functional data context is based on the ideas of Cuesta-Albertos et al. (2007) [25]. These authors give an extension, on Hilbert spaces, to the Cramer-Wold theorem, which characterizes a probability distribution in terms of one-dimensional projections. Their Theorem 4.1 states that if two distributions are different, and we choose a marginal of them at random, those marginals will almost surely be different. Based on this fact, Cuesta-Albertos and Febrero-Bande [20] propose a method for hypothesis testing in infinite dimensional spaces. We will state their result more formally, particularizing it to our problem. Let us assume the data belong to a Hilbert space,, with a scalar product 〈, 〉, and let μ_{ G } be a Gaussian distribution in. Suppose the hypothesis to be tested is whether certain parameters, say γ_{1} and γ_{2}, are equal (H_{0}:γ_{1}=γ_{2}). If γ_{1}≠γ_{2}, then the set of random directions, ν, from μ_{ G } in, for which 〈ν,γ_{1}〉=〈ν,γ_{2}〉, has probability zero. That is, if H_{0} fails, then it also fails in its projected version for μ_{ G }—almost every ν∈. Therefore, a test at level a in the one-dimensional space is a test at the same level to test H_{0}.
- 1.If ∃k∈{1,2}, such that α _{ k }≠0, then${\mu}_{G}(\{\nu \in \mathcal{\mathscr{H}}:\u3008\nu ,{\alpha}_{k}\u3009=0\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\forall k\in \{1,2\}\})=0$
- 2.If ∃g∈{1,2,3,4,5}, such that β _{ g }≠0, then${\mu}_{G}(\{\nu \in \mathcal{\mathscr{H}}:\u3008\nu ,{\beta}_{g}\u3009=0\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\forall g\in \{1,2,\dots ,5\}\})=0$
- 3.If ∃k∈{1,2} and g∈{1,2,3,4,5} such that γ _{ k g }≠0, then$\phantom{\rule{-4.0pt}{0ex}}{\mu}_{G}(\{\nu \in \mathcal{\mathscr{H}}:\u3008\nu ,{\gamma}_{\mathit{\text{gk}}}\u3009=0\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\forall k\in \phantom{\rule{0.3em}{0ex}}\{1,2\},\forall g\in \phantom{\rule{0.3em}{0ex}}\{1,2,\dots ,5\}\})=0$
The tests defined in the one-dimensional response case are clearly conditional on the random projection used, ν. To reduce the error introduced by the choice of the random projection, we will use the correction that implies controlling the false discovery rate (FDR) introduced by Benjamini and Hochberg (1995) [26]. This method is also recommended by Cuesta-Albertos and Febrero-Bande [20]. In particular, we use the FDR procedure that arises from the work of Benjamini and Yakutieli (2001) [27]. Given the ordered p-values p_{(1)}<…<p_{(s)} obtained using s random projections, we will choose the corrected p-value as the quantity $\text{min}\left\{\frac{s}{i}{p}_{(i)},i=1,\cdots \phantom{\rule{0.3em}{0ex}},s\right\}$, where min stands for the minimum of a set.
So far, this test can help in the search for global differences between the two groups of curves, although we would rather study how these differences change in time. To do this, we propose the use of moving windows along time. For each time point, t, we consider an interval of time, centered at t, and project the pieces of curves that correspond to that interval, to perform the ANOVA test.
The hypotheses in (6) can be tested with any regular ANOVA approach. Nevertheless, some caution is needed, as the errors in our model cannot be assumed to be independent, as we discuss in the next subsection.
ANOVA model with dependent errors
In this section we introduce the problem of the dependence that is present in the data. The dependency comes from the fact that the data are observed at the neuron-pair level. So, it is only fair to think that the curves obtained from two pairs of neurons with one cell in common could be correlated.
With this definition, Σ is positive definite for small values of ρ. The eigenvalues of the matrix in (10) result in only three different values: (1−2ρ)σ^{2}, (1+4ρ)σ^{2} and (1+12ρ)σ^{2}, which, for ρ∈[0,1], are all greater than zero at the same time, whenever ρ<0.5. This means that this correlation structure for the data is plausible only for ρ∈ [ 0,0.5).
The Fstatistic
and, in both cases, h is a vector of zeros. Let us denote by q the total number of parameters in the complete model and let q_{0} be the number of parameters, different from zero, under the null hypothesis.
where d_{1} denotes the degrees of freedom of the numerator, that is, q−q_{0}, and d_{2} the degrees of freedom of the denominator: N−(q+1).
The test statistic, F, under assumptions of independence, normality and homoscedasticity, follows an F distribution with q−q_{0} and N−(q+1) degrees of freedom. In our context, as the errors of the model cannot be assumed to be independent, the test statistic does not have an ${F}_{{d}_{1},{d}_{2}}$ distribution and, therefore, we need to calibrate the null distribution of the test statistic. We propose a parametric bootstrap procedure to calibrate the distribution of the ANOVA test statistic under the null hypothesis.
Estimation of the correlation coefficient
Direct bootstrap
As the data are normally distributed, we propose the use of a parametric bootstrap to calibrate the distribution of the ANOVA test statistic under the null hypothesis. We now describe the procedure for a general null hypothesis ${H}_{0}^{i}:{\theta}_{i}=0$.
Once the linear model has been fitted to obtain $\widehat{\theta}$ and the classical ANOVA statistic has been computed (denoted by F^{ o b s }), we estimate σ^{2} and ρ from the residuals to build the estimated covariance matrix, $\widehat{\mathbf{\Sigma}}=\mathbf{\Sigma}({\widehat{\sigma}}^{2},\widehat{\rho})$. We proceed with the following bootstrap algorithm:
- 1.
Replace the i-th parameter in the estimated $\widehat{\mathbf{\theta}}$ by zero (null hypothesis). This set of parameters will be denoted by ${\widehat{\mathbf{\theta}}}_{0}$. Build a bootstrap sample: ${\mathbf{y}}^{\ast}=\mathbf{X}{\widehat{\theta}}_{0}+\mathbf{z}$ with $\mathbf{z}\sim N(0,\widehat{\mathbf{\Sigma}})$.
- 2.
Fit the linear model to the resample, obtain the bootstrap version of the estimated parameters ${\widehat{\theta}}^{\ast}$ and compute F^{∗}, the bootstrap version of F.
- 3.
Repeat Steps 1–2 B times to obtain F^{∗1},…,F^{∗B}.
- 4.
Compute the desired (1−a)-quantile of the bootstrap versions, ${\mathbf{F}}_{1-a}^{\ast}$. We reject the null hypothesis if ${\mathbf{F}}^{\mathit{\text{obs}}}>{\mathbf{F}}_{1-a}^{\ast}$.
Results and discussion
On the other hand, we would like to have directions without tendency and such that their variability throughout the trajectory does not change too much. For this aim, we define the random trajectories as the sum of two Brownian motions, as just defined, where one of them has been “flipped" so as to be equal to zero in the last time point t_{ M }. That is, let ν_{1} and ν_{2} be defined as in (14) and let ν_{3}(t_{ k })=ν_{2}(t_{M−k+1}). The final directions we use are defined as ν(t)=ν_{1}(t)+ν_{3}(t).
When large correlation is present (ρ>0.5) an alternative, nonparametric, bootstrap can be carried out. In the following procedure, each set of bootstrap residuals is defined as the set of original residuals of a trial chosen at random with equal probability from all possible trials:
1. For each k∈{1,2} and l∈{1,2,3,4} draw the bootstrap pair (k^{∗},l^{∗}) with equal probability from {1,2}×{1,2,3,4}, that is, $P({k}^{\ast}={k}^{\prime},{l}^{\ast}={l}^{\prime})=\frac{1}{8}$ for all k^{′}∈{1,2,} and all k^{′}∈{1,2,3,4}.
2. Define ${\widehat{\epsilon}}_{k,l}^{(i,j)\ast}={\widehat{\epsilon}}_{{k}^{\ast},{l}^{\ast}}^{(i,j)}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\forall (i,j)$.
This bootstrap procedure has the drawback that, in our case, the vector of bootstrap residuals can only take eight possible values. A possible improvement is to use a smoothed version. To achieve this, a smoothing parameter h, typically small with respect to the standard deviation of the residuals, is chosen and Step 2 is replaced by:
2. Define ${\widehat{\epsilon}}_{k,l}^{(i,j)\ast}={\widehat{\epsilon}}_{{k}^{\ast},{l}^{\ast}}^{(i,j)}+{\mathit{\text{hZ}}}_{k,l}^{(i,j)\ast}$ with ${Z}_{k,l}^{(i,j)\ast}\sim N(0,1)$ iid ∀(i,j).
The green curve in Figure 2 represents the p-values obtained with an alternative bootstrap algorithm based on an approximation to the real distribution of the F statistic. This alternative method, which we call the χ^{2}-based bootstrap, reduces computational time considerably. It involves some theoretical insight that we describe in the Appendix. The main result is based on the fact that the test statistic F can be expressed as a ratio of quadratic forms on the errors of the model. There are plenty of results concerning the exact distribution of quadratic forms on normal vectors. For example, Provost et al. [28] derive the exact density function of an indefinite quadratic form in noncentral vectors, which allows us to derive the distribution of ratios of quadratic forms as those involved in ANOVA tests. Nevertheless, the closed formulas for the density functions of quadratic forms are complex and not practical. On the other hand, the distribution of the numerator and the denominator that give shape to the F statistic are quite easy to approximate by Montecarlo, as is also described in the Appendix. We use this approach to carry out, in the next two sections, an evaluation of the test. First, we compare the distribution of the test statistic when dependence is either taken into account or not. Finally, we perform a simulation study to evaluate the performance of the bootstrap test.
Distribution of the test statistic F
Performance of the test
where θ=(m,α_{1},β_{1})^{ t }. The covariance matrix, Σ, was defined as in (9) for four values of ρ. The values used for the simulations are m=0, α_{1}=0,0.25,0.5, β_{1}=0,0.25,0.5 and ρ=0,0.05,0.15,0.35.
Proportion of rejections for ${H}_{0}^{\phantom{\rule{0.4pt}{0ex}}\alpha}$ at level a =0 . 05
ρ=0 | ρ=0 .05 | ρ=0 .15 | ρ=0 .35 | ||||||
---|---|---|---|---|---|---|---|---|---|
Boot. | F | Boot. | F | Boot. | F | Boot. | F | ||
β_{1}=0 | 0.0456 | 0.0512 | 0.0804 | 0.1156 | 0.0828 | 0.2166 | 0.0816 | 0.3692 | |
α_{1}=0 | β_{1}=0.25 | 0.0378 | 0.0464 | 0.0764 | 0.1102 | 0.0874 | 0.2258 | 0.0806 | 0.3680 |
β_{1}=0.5 | 0.0464 | 0.0536 | 0.0756 | 0.1092 | 0.0862 | 0.2178 | 0.0848 | 0.3780 | |
β_{1}=0 | 0.7570 | 0.7906 | 0.6482 | 0.7476 | 0.4864 | 0.7044 | 0.3164 | 0.6762 | |
α_{1}=0.25 | β_{1}=0.25 | 0.7564 | 0.7920 | 0.6590 | 0.7598 | 0.5022 | 0.7124 | 0.3204 | 0.6644 |
β_{1}=0.5 | 0.7548 | 0.7886 | 0.6342 | 0.7390 | 0.4860 | 0.7050 | 0.3164 | 0.6810 | |
β_{1}=0 | 0.9996 | 0.9998 | 0.9940 | 0.9982 | 0.9430 | 0.9882 | 0.7794 | 0.9602 | |
α_{1}=0.5 | β_{1}=0.25 | 0.9996 | 1.0000 | 0.9944 | 0.9988 | 0.9478 | 0.9970 | 0.7858 | 0.9596 |
β_{1}=0.5 | 0.9998 | 0.9998 | 0.9950 | 0.9986 | 0.9460 | 0.9898 | 0.7666 | 0.9568 |
Proportion of rejections for ${H}_{0}^{\beta}$ at level a =0 . 05
ρ=0 | ρ=0 .05 | ρ=0 .15 | ρ=0 .35 | ||||||
---|---|---|---|---|---|---|---|---|---|
Boot. | F | Boot. | F | Boot. | F | Boot. | F | ||
α_{1}=0 | 0.0588 | 0.0532 | 0.0558 | 0.0462 | 0.0444 | 0.0238 | 0.0376 | 0.0032 | |
β_{1}=0 | α_{1}=0.25 | 0.0566 | 0.0538 | 0.0512 | 0.0448 | 0.0462 | 0.0244 | 0.0414 | 0.0032 |
α_{1}=0.5 | 0.0516 | 0.0460 | 0.0484 | 0.0422 | 0.0448 | 0.0248 | 0.0374 | 0.0040 | |
α_{1}=0 | 0.7936 | 0.7844 | 0.8094 | 0.7926 | 0.8838 | 0.8358 | 0.9758 | 0.8858 | |
β_{1}=0.25 | α_{1}=0.25 | 0.7938 | 0.7864 | 0.8280 | 0.8072 | 0.8716 | 0.8152 | 0.9796 | 0.9020 |
α_{1}=0.5 | 0.7830 | 0.7740 | 0.8118 | 0.7940 | 0.8824 | 0.8298 | 0.9802 | 0.9038 | |
α_{1}=0 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 0.9998 | 1.0000 | 1.0000 | |
β_{1}=0.5 | α_{1}=0.25 | 0.9998 | .9998 | 0.9998 | 0.9996 | 1.0000 | 1.0000 | 1.0000 | 1.0000 |
α_{1}=0.5 | 1.0000 | 1.0000 | 0.9998 | 0.9996 | 1.0000 | 1.0000 | 1.0000 | 1.0000 |
Proportion of rejections for ${H}_{0}^{\gamma}$ at level a =0 . 05
ρ=0 | ρ=0 .05 | ρ=0 .15 | ρ=0 .35 | ||||||
---|---|---|---|---|---|---|---|---|---|
Boot. | F | Boot. | F | Boot. | F | Boot. | F | ||
β_{1}=0 | 0.0558 | 0.532 | 0.0554 | 0.452 | 0.0454 | 0.0254 | 0.0382 | 0.0036 | |
α_{1}=0 | β_{1}=0.25 | 0.0568 | 0.532 | 0.0506 | 0.426 | 0.0458 | 0.0248 | 0.0386 | 0.0052 |
β_{1}=0.5 | 0.0500 | 0.474 | 0.0476 | 0.394 | 0.0478 | 0.0264 | 0.0410 | 0.0038 | |
β_{1}=0 | 0.0610 | 0.0552 | 0.0510 | 0.0414 | 0.0482 | 0.0274 | 0.0396 | 0.0032 | |
α_{1}=0.25 | β_{1}=0.25 | 0.0530 | 0.0500 | 0.0494 | 0.0402 | 0.0442 | 0.0260 | 0.0334 | 0.0036 |
β_{1}=0.5 | 0.0502 | 0.0500 | 0.0496 | 0.0414 | 0.0472 | 0.0280 | 0.0374 | 0.0028 | |
β_{1}=0 | 0.0516 | 0.0464 | 0.0476 | 0.0408 | 0.0424 | 0.0236 | 0.0330 | 0.0024 | |
α_{1}=0.5 | β_{1}=0.25 | 0.0582 | 0.0540 | 0.0516 | 0.0430 | 0.0440 | 0.0242 | 0.0414 | 0.0056 |
β_{1}=0.5 | 0.0524 | 0.0480 | 0.0552 | 0.0438 | 0.0500 | 0.0298 | 0.0370 | 0.0052 |
Finally, we show similar simulation results when testing for the interaction between the two factors. Table 3 shows the proportion of rejections under the null hypothesis for different values of ρ and 5000 Montecarlo replications. We can observe that the nominal level is successfully met by the bootstrap test for small correlation values. Although for the case of large correlation values the test seems to be conservative, it outperforms the F test remarkably. It is important to note that, regarding the power of the test, in these simulations the test rejected the null hypothesis 100% of the time (results not shown), using values of γ=0.5,1 to simulate data in the alternative.
Conclusions
In this work we proposed a functional two-way ANOVA for neural synchrony curves, using random projection techniques. These methods are very easy to implement and interpret, which makes them appealing for applying to many problems. The method was shown to be useful as it allowed the significance effects of the factors under study to be addressed.
The model under study involves synchrony curves obtained by a cross-correlation based method. The curves were separated into groups given by the stimuli and the difference in preferred orientation between the two neurons involved in each curve. Differences between the levels of this second factor were also of interest. Although there were groups with very few elements, several conclusions can be established.
The importance of including the dependence between curves in the analysis was shown. The distribution of the test statistic was approximated using a parametric bootstrap on the residuals of the model, allowing for dependence. The classical F test statistic can lead to false positives, as the distribution of the test statistic has a heavier tail than the F distribution. Two algorithms were presented to carry out the bootstrap: a direct resampling plan, which resamples from the model, and a second one, based on the fact that the F statistic has the same distribution as that of a ratio of particular linear combinations of ${\chi}_{1}^{2}$ variables. The second algorithm has a substantially lower computational cost than the first one.
Regarding the real data, the interactions between the factors were not statistically significant. An effect of the factor stimulus can be found at the beginning of the awake-like period (a few seconds after the stimulus onset). To make a statement regarding the differential effect that bs and bf have on pairwise synchrony, the analysis of more groups of neurons would be necessary. However, this work shows that the question is worth asking, as we were able to find evidence of these differences after the switch from anesthesia to the awake-like mode. In fact, taking into account the ability of bs and bf to promote wakefulness with a similar effect on ECoG power spectral density [29], the robustness of the proposed statistical method (to be applied in low-activity cortical single units recorded simultaneously) allows us to find some differences in neural pattern synchrony, consistent with physiological data as follows. Activation of the parabrachial nucleus in the bs enables thalamic relay neurons to disrupt cortical synchronization via glutamate release [30]. In contrast, bf stimulation induces not only ACh release, but also GABA, glutamate and NO in the cortex [18, 31]. Thus, although the effect that bs and bf activation has on ECoG dynamics has been characterized as regulatory on cortical activation, the different operational mechanism of each system could be reflected in the temporal coordination between neurons. Finally, a relevant issue in neurophysiology relates to the type of information processing carried on by primary visual cortex neurons [32–34]. Briefly, some basic properties of the visual processing—like orientation selectivity—could be based either on highly refined and specific anatomical connections or, in contrast, could be carried out by distributed computational processes at the cortical level. The results presented here show that the effect of the difference in orientation selectivity was significant throughout all the generated awake-like activity, suggesting that the strength of the connectivity was dependent on the orientation selectivity of primary visual cortex neurons, thus favoring the second hypothesis.
Overall, this work is a contribution to the development of statistical tools for neuroscience. Although the methods we propose here have been focused on a particular problem, they are applicable in many other contexts. It is worth to mention that differences were found even though the firing activity of the test group was considerably low, showing the method can be useful when low firing rates are present. Also, bootstrap techniques are very powerful and easy to implement (although, admittedly, they can be computationally expensive), and might be a good alternative to parametric inference using minimal mathematical assumptions.
Appendix - χ^{2}-based bootstrap
where ε∼N(0,Σ).
where ${V}_{i}^{2}\sim {\chi}_{1}^{2}$ for i=1,…,N.
has the same distribution as F. The distribution of (20) can be approximated by Montecarlo.
In practice, Σ is not known, so the λ_{ i } and μ_{ i } cannot be computed. However, Σ can be replaced by $\hat{\Sigma}=\mathbf{\Sigma}({\widehat{\sigma}}^{2},\widehat{\rho})$ as in (10). Thus, ${\widehat{\lambda}}_{i}$ and ${\widehat{\mu}}_{i}$ (the eigenvalues of $\hat{\Sigma}{A}_{1}$ and $\hat{\Sigma}{A}_{2}$) can be used in (20). Observe that this is equivalent to the original bootstrap resampling plan as presented above. This gives an alternative method for performing the bootstrap that is much less time consuming. For instance, in the example shown in Figure 2, the results obtained with the direct bootstrap are reproduced, except for minor differences due to Montecarlo error, and the computational time is at least 30 times smaller with this alternative method.
Declarations
Acknowledgments
This work was supported by grant MTM2011-22392 from the Spanish Ministry of Economy and Competitiveness. AMGM was supported by a “Formación de Personal Investigador” (FPI) grant (BES-2009-017772) from the Spanish Ministry of Economy and Competitiveness. AMGM, NE and JM were supported by Xunta de Galicia (INCITE09 137 272 PR). The authors would like to thank Professor Juan Cuesta-Albertos for his useful comments during the early stage of this research.
Authors’ Affiliations
References
- Kuhn A, Aertsen A, Rotter S: Higher-order statistics of input ensembles and the response of simple model neurons. Neural Comput. 2003, 15: 67-101.View ArticlePubMed
- Steinmetz PN, Roy A, Fitzgerald PJ, Hsiao SS, Johnson KO, Niebur E: Attention modulates synchronized neuronal firing in primate somatosensory cortex. Nature. 2000, 404: 187-190.View ArticlePubMed
- Stopfer M, Bhagavan S, Smith BH, Laurent G: Impaired odour discrimination on desynchronization of odour-encoding neural assemblies. Nature. 1997, 390: 70-74.View ArticlePubMed
- Vaadia E, Haalman I, Abeles M, Bergman H, Prut Y, Slovin H, Aertsen A: Dynamics of neuronal interactions in monkey cortex in relation to behavioral events. Nature. 1995, 373: 515-518.View ArticlePubMed
- Zohary E, Shadlen M, Newsome W: Correlated neuronal discharge rate and its implications for psychophysical performance. Nature. 1994, 370: 140-143.View ArticlePubMed
- Gerstein G, Perkel D: Simultaneously recorded trains of action potentials - Analysis and functional interpretation. Science. 1969, 164: 828-830.View ArticlePubMed
- Grün S, Diesmann M, Aertsen A: Unitary events in multiple single-neuron spiking activity: I. Detection and significance. Neural Comput. 2002, 14: 43-80.View ArticlePubMed
- Grün S, Diesmann M, Aertsen A: Unitary events in multiple single-neuron spiking activity: II. Nonstationary Data. Neural Comput. 2002, 14: 81-119.View ArticlePubMed
- Faes C, Geys H, Molenberghs G, Aerts M, Cadarso-Suárez C, Acuña C, Cano M: A flexible method to measure synchrony in neuronal firing. J Am Stat Assoc. 2008, 103: 149-261.View Article
- González-Montoro AM, Cao R, Faes C, Molenberghs G, Espinosa N, Cudeiro J, Mariño J: Cross nearest-spike interval based method to measure synchrony dynamics. Math Biosci Eng. 2014, doi:10.3934/mbe.2014.11.27.
- Grün S: Data-driven significance estimation for precise spike correlation. J Neurophysiol. 2009, 101: 1126-1140.PubMed CentralView ArticlePubMed
- Ventura V, Cai C, Kass RE: Statistical assessment of time-varying dependence between two neurons. J Neurophysiol. 2005, 94: 2940-2947.View ArticlePubMed
- Steriade M: Sleep oscillations and their blockage by activating systems. J Psychiat Neurosci. 1994, 19: 354-358.
- Steriade M, Gloor P, Llinás RR, Lopes da Silva FH, Mesulam MM: Basic mechanisms of cerebral rhythmic activities. Electroencephalogr Clin Neurophysiol. 1990, 76: 481-508.View ArticlePubMed
- Brown EN, Lydic R, Schiff D: General anesthesia, sleep, and coma. N Engl J Med. 2010, 363: 2638-2650.PubMed CentralView ArticlePubMed
- Franks NP: General anaesthesia: from molecular targets to neuronal pathways of sleep and arousal. Net Rev Neurosci. 2008, 9: 370-386.View Article
- Sleigh JW, Voss L, Steyn-Ross ML, Steyn-Ross DA, Wilson MT: Modelling sleep and general anaesthesia. Springer Ser Comp Neurosci. 2011, 15: 21-41.
- Mariño J, Cudeiro J: Nitric oxide-mediated cortical activation: A diffuse wake-up system. J Neurosci. 2003, 23: 4299-4307.PubMed
- Hubel DH, Wiesel TN: Receptive fields, binocular interaction and functional architecture in the cat’ s visual cortex. J Physiol. 1962, 10: 106-154.View Article
- Cuesta-Albertos J, Febrero-Bande M: A simple multiway ANOVA for functional data. Test. 2010, 19: 537-557.View Article
- Perkel DH, Gerstein GL, Moore GP: Neuronal spike trains and stochastic processes. II. Simultaneous spike trains. Biophys J. 1967, 7: 419-440.PubMed CentralView ArticlePubMed
- Ramsay JO, Silverman BW: Applied Functional Data Analysis: Methods and Case Studies. 2002, New York: SpringerView Article
- Ramsay JO, Silverman BW: Functional Data Analysis. 2005, New York: Springer Series in StatisticsView Article
- Ferraty F, Vieu P: Nonparametric Functional Data Analysis: Theory and Practice. 2006, New York: Springer
- Cuesta-Albertos JA, Fraiman R, Ransford T: A sharp form of the Cramer-Wold theorem. J Theor Prob. 2007, 20: 201-209.View Article
- Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995, 57: 289-300.
- Benjamini Y, Yekutieli D: The control of the false discovery rate in multiple testing under dependency. Ann Stat. 2001, 29: 1165-1188.View Article
- Provost SB, Ridiuk EM: The exact distribution of indefinite quadratic forms iin noncentral normal vectors. Ann Inst Statist Math. 1996, 48: 381-394.View Article
- Fuller PM, Sherman D, Pedersen NP, Saper CB, Lu J: Reassessment of the structural basis of the ascending arousal system. J Comp Neurol. 2011, 519: 933-956.PubMed CentralView ArticlePubMed
- Steriade M: Arousal: revisiting the reticular activating system. Science. 1996, 272: 225-226.View ArticlePubMed
- Henny P, Jones BE: Projections from basal forebrain to prefrontal cortex comprise cholinergic, GABAergic and glutamatergic inputs to pyramidal cells or interneurons. Eur J Neurosci. 2008, 27: 654-670.PubMed CentralView ArticlePubMed
- Schummers J, Mariño J, Sur M: Synaptic integration by V1 neurons depends on location within the orientation map. Neuron. 2002, 36: 969-978.View ArticlePubMed
- Mariño J, Schummers J, Lyon DC, Schwabe L, Beck O, Wiesing P, Obermayer K, Sur M: Invariant computations in local cortical networks with balanced excitation and inhibition. Nat Neurocsi. 2005, 8: 194-201.View Article
- Stimberg M, Wimmer K, Martin R, Schwabe L, Mariño J, Schummers J, Lyon DC, Sur M, Obermayer K: The operating regime of local computations in primary visual cortex. Cereb Cortex. 2009, 19: 2166-2180.PubMed CentralView ArticlePubMed
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.