Bivariate and Multivariate NeuroXidence: A Robust and Reliable Method to Detect Modulations of Spike–Spike Synchronization Across Experimental Conditions

Synchronous neuronal firing has been proposed as a potential neuronal code. To determine whether synchronous firing is really involved in different forms of information processing, one needs to directly compare the amount of synchronous firing due to various factors, such as different experimental or behavioral conditions. In order to address this issue, we present an extended version of the previously published method, NeuroXidence. The improved method incorporates bi- and multivariate testing to determine whether different factors result in synchronous firing occurring above the chance level. We demonstrate through the use of simulated data sets that bi- and multivariate NeuroXidence reliably and robustly detects joint-spike-events across different factors.


IntroductIon
Synchronous neuronal spiking has been proposed as a code that carries, groups, or binds information (von der Malsburg, 1981;Victor and Purpura, 1996;Fetz, 1997;Singer et al., 1997;Singer, 1999;Pillow et al., 2008;Uhlhaas et al., 2009); supports information processing and selection (Fries et al., 2001(Fries et al., , 2002Womelsdorf and Fries, 2007); and maintains and recollects information for shortterm memory (Compte et al., 2000;Miller et al., 2003). However, even though this synchrony-coding hypothesis has been investigated intensively for the last two decades, its actual role and importance are still controversially debated. One of the reasons for this ongoing debate has been the lack of tools to identify synchronous activity and to decipher its relative importance compared to complementary features of neuronal spiking, such as the firing rate (Barlow, 1972;Bialek et al., 1991;Bialek and Rieke, 1992;Newsome, 1994, 1998;Brody, 1999;Oram et al., 1999;Baker and Lemon, 2000;Nakahara and Amari, 2002). A second reason for the continuing debate is that one must be able to demonstrate that synchronous activity is indeed involved in information processing.
Recently, however, very important steps have been taken to address these two points. In order to identify relevant synchronous spiking activity, recently developed tools enable the tracking of the temporal modulation of neuronal synchronization (Brown et al., 1998(Brown et al., , 2004Grün et al., 1999Grün et al., , 2002aGrün et al., ,b, 2003Martignon et al., 2000;Aertsen et al., 2001;Kuhn et al., 2003;Kass et al., 2005;Truccolo et al., 2005;Schneider et al., 2006) and the comparison of this temporal modulation to the modulations of other properties of neuronal firing, such as the firing rate (Riehle et al., 1997;Pipa et al., 2007). The use of these tools allows one to demonstrate that the synchronization of spiking activity is often tightly linked to relevant temporal periods of an experiment, which indicates that neuronal synchronization is correlated with processes that underlie neuronal information processing. In 2008, the method NeuroXidence was stochasticity of spike times, a JSE is defined by a model pattern of precise spike times and a millisecond-scale temporal window, which accounts for the maximal uncertainty around each of these exact spike times (Figure 1). Surrogate data, derived by jittering each individual spike train in the original data, is used to estimate the frequencies of chance JSEs under the H 0 assumption that neurons are not coupled on a fine temporal time scale. For each trial and each experimental/ behavioral factor, the differences in the frequencies of JSEs between the original and the surrogate data sets are computed. Finally, bi-and multi-variate NeuroXidence tests whether the mean or median of the frequencies of JSEs is significantly different across factors.

MaterIals and Methods
In order to determine the relationship between the frequency of a certain JSE and the state of the neuronal system, the bi-and multi-variate NeuroXidence contains two main steps: detection and comparison.
The detection step identifies the frequency of JSEs in the original and surrogate data sets. The comparison step tests whether the frequencies of a particular JSE, which are estimated in the detection step, vary across different experimental/behavioral factors. Here, the null hypothesis H 0 states that the precise timing of neuronal firing is irrelevant to the state of the neuronal system. The alternative hypothesis (H 1 ) assumes that there are more (or less) coordinated firing events according to different conditions of the system. Therefore, different statistical tests (t-test, Mann-Whitney U test, etc.) are applied to the frequencies of JSEs from different factors to test the null hypothesis.

detectIon of Jses
In order to detect JSEs in simultaneously recorded spike trains, we define two time-scale parameters, t c and t r . The time scale of synchronous firing determines the parameter t c , which ranges between 1 and 10 ms and defines the expected systematic precision of a JSE. Time scale t c is, therefore, equivalent to the maximal deviation of the individual spikes in the assumed model pattern (Figure 1). The parameter t r defines the lower bound of the rate modulation and is distinguishably slower than t c by a factor of h, which ranges from 2 to 5. Next, surrogate data sets are generated by jittering the original individual spike trains by a random interval t t, n (t indexes the trials, n indexes the neurons), which is determined by either a uniform or Gaussian distribution. The extent of the jittering is on the slow time scale t r , which destroys any fine temporal crossstructure between the spike trains that is less than t r , but maintains the auto-structure and other features of each individual spike train.
We define the identity of a JSE by the set of neurons that exhibits the coordinated firing, and we refer to this set of neurons as a jointspike pattern (JS-pattern). Here we distinguish the complexity of a JSE (Shadlen and Newsome, 1994), which is defined by the number of neurons participating in the spike pattern, from the order of correlation (Singer et al., 1997;Shadlen and Newsome, 1998;Singer, 1999), which is defined by the number of neurons that are directly coupled and are the primary cause of the coordinated firing. Due to an excess or lack of one or more chance events, the order is not necessarily the same as the complexity.
Between the original and the surrogate data sets, we compute the difference in the raw counts of how many times a specific JS-pattern k appears in trial t of experimental/behavioral factor j. This difference is expressed as j ∈ 0 and ∆f Q t k ,j ∈ for trial t = 1,2, …,T; JS-pattern k = 1,2, …,K; surrogate s = 1,2, …,S; and experimental/behavioral factor j = 1,2, …., M.

coMparIson of Jses across dIfferent factors
There exists an underlying distribution of frequencies of JSEs, Φ j k , for each specific JS-pattern k and experimental/behavioral factor j. Under the bivariate condition, we wish to compare the distributions for factors j = 1, 2. For each factor j, we approximate the distribution Φ j k with the set of mean frequencies defined by equation (1) taken across all trials t, where ∆f t k ,j is a sample from trial t taken from the set ∆F k j for a given experimental/behavioral factor. Because each sample in a given set is independent and identically distributed, the values of ∆Θ t bi k , computed for all the trials can form a set ∆Θ Θ ∆Θ ∆ Θ bi k bi k bi …, }. Modified from a previously published figure (Pipa et al., 2008).
In order to study the effect of rate covariation on the estimation of JSEs, two scenarios are used to generate the rate profile for each Poisson model. Scenario I: Two simulated Poisson models exhibit rate covariation, such that the rate profiles of the two simulated data sets are identical. Scenario II: Two simulated Poisson models exhibit no rate covariation, such that the rate profiles of the two simulated data sets are different. For example, the mean firing rate of one simulated data set is chosen as 15 spikes/s, while the other changes from 7 to 90 spikes/s.
For both scenarios, bivariate NeuroXidence is applied to detect an excess or lack of JSEs between the two sets of model data (Figures 2  and 3). Based on the parameter analysis and Appendix 6 recommendations from the original NeuroXidence publication by Pipa et al. (2008), the standard parameters for both scenarios include a sliding window with a duration of l = 200 ms, T = 50 trials, a mean spike rate of 15 spikes/s, and h equals 5. Scenario I utilizes 20 surrogates samples, and 15 different combinations of parameters are derived for the model data by taking all combinations of the number of trials (T = 20, 50, 100) and the covarying mean spiking rates (r 1 = r 2 = 7, 10, 30, 60, 90 spikes/s). Scenario II uses 1 surrogate sample, and 15 different parameter combinations are derived by varying the number of trials (T = 20, 50, 100) and the second mean spiking rate (r 1 = 15, r 2 = 7, 10, 30, 60, 90 spikes/s). In both scenarios, the simulated data sets are generated in order to determine the false-positive rate for JS-patterns of complexities 2-6.
In both scenario I and II, in order to compute the false-positive rate, eight simulated spike trains are generated for each condition, which are taken to be signals from eight channels. The number of JS-patterns for each condition is defined by combinations of the different spike trains. In this case, there are 247 JS-patterns, including all second-order to eighth-order JS-patterns. By comparing the occurrence frequencies of these 247 JS-patterns across two conditions with a selected test level, some of the JS-patterns shows a significant difference between the two conditions. After gathering this information for JS-patterns of the same complexity, one can determine the false-positive rate for each complexity of JS-pattern. This procedure is repeated 200 times resulting in the final percentile plots for scenarios I and II (Figures 2 and 3).
In both scenarios, none of the results from the model data sets for any of the parameter combinations exceed the chance level, for either a test level of 5% (Figures 2A,B and 3A The false-positive rates are much smaller for the high-complexity JS-patterns compared to the rates for low-complexity JS-patterns. When there is no rate-covariation between two conditions, as in scenario II, increasing the number of surrogates above 1 leads to a biased estimate of the excess or lack of JSEs due to the skewness of the distribution of JSEs. Setting S = 20 for non-covarying rates results in a false-positive rate that is above the chance level (data not shown). Therefore, if there is no rate-covariation between two conditions, the surrogate number should be set to 1 to ensure a reasonable false-positive rate. By properly adjusting the number is not zero, there is an excess in the frequency of occurrence of JS-pattern k for trial t in one experimental/behavioral factor or the other. The statistical significance of the differences in the means or medians of the sampled distributions ∆Θ bi k is evaluated by t-test or Mann-Whitney U test, respectively. To account for distributions that might have different variances, the Behrens-Fisher t-test is applied instead of the normal Student's t-test. All future references to t-tests in this manuscript refer to t-tests with unequal variances.
In analogy with the bivariate condition, under the multivariate condition the underlying distribution of frequencies of JSEs Φ j k is approximated by forming the sets ∆F k j . The statistical significance of the differences in the means and medians of the sampled distributions Φ j k across factors is evaluated by analysis of variance (ANOVA) or a Kruskal-Wallis one-way test, respectively. To account for distributions that might have different variances, the Kruskal-Wallis one-way test is highly recommended over an ANOVA. In some simulations, we apply an ANOVA test during the method calibration in order to verify whether or not the test will fail.
Finally, after gathering all of the significance levels for each JS-pattern k, one can evaluate which factor plays the most important role in modulating the strength of synchronous firing.

results
We demonstrate the robustness and sensitivity of this bi-and multivariate extension of NeuroXidence in two stages. We first present results for the bivariate case that tests the modulation of spike-spike synchronization across exactly two conditions or factors. Secondly, we examine the multivariate case, which tests for modulations of spike-spike synchronization across more than two conditions.

part I: BIvarIate neuroXIdence
We establish the robustness of the bivariate method, by showing false-positive rates estimated for different spike-train models. Next, we reveal the sensitivity and test power of the presented method and discuss the detectability of sub-patterns and supra-patterns of synchronized neuronal ensembles. In the last step, we demonstrate the robustness of the method when applied to simulated spike trains that contain various properties, which often occur in real neuronal firing activities, such as rate covariations, and regular or bursty firing.

False positives for two stationary processes
Two simulated data sets are generated for the bivariate case. The generation of simulated data is described in detail in Pipa et al. (2008), but we provide a brief description here. Each of the two simulated data sets contains eight neurons, and each spike train is generated by an independent and homogeneous Poisson process. Such spike trains are generated for each of the two experimental conditions being examined. For each condition, we apply NeuroXidence to detect the occurrence frequencies of JSEs for each JS-pattern and for S surrogates. Next, the test level is set for comparing the corrected frequencies of JSEs across the two experimental conditions. Therefore, NeuroXidence yields an estimate of whether or not JS-patterns show a significant difference across the two conditions.  The standard parameter set is chosen as 50 trials (T), a background spike rate of 15 spikes/s (r), 20 surrogates (S), and h equals 5. To study how the test power of bivariate NeuroXidence is affected by the number of trials (T = 20, 50, 100), the background spike rates (r 1 = r 2 = 10, 15, 30, 60, 90 spikes/s), and the number of surrogates (S = 1, 20, 50), each parameter is varied in turn. The test power is derived for JS-patterns of complexity 2-6, with a test level set to 5%. It is clear in Figure 4 that an increase in the number of trials leads to an increase of the test power, as exhibited by the leftward shift and increase in steepness of the curves. With an increase in the background rate, more JSEs are required to achieve a high test power, as demonstrated by the rightward shift and decrease in steepness of the curves. However, for higher complexity patterns, this effect is reduced. For lower complexity patterns, a larger number of surrogates (for example, S = 20) is required to reliably detect the excess of surrogates in this way, bivariate NeuroXidence provides a conservative method for detecting the excess or lack of JSEs between two groups of neurons, no matter whether rate-covariations exist between the groups or not.

Test power for two stationary processes
To assess the test power of bivariate NeuroXidence, one simulated data set is generated by a single-interaction process based on a Poisson process (Fries et al., 2001), while the other is generated by an independent and homogenous Poisson process. Each data set contains eight simultaneous spike trains. For the single-interaction process, correlated spike trains were characterized by a background rate, which corresponds to the independent spiking of neurons, and by a JSE rate, which defines the expected frequency of the JS-pattern of interest.

False positives for a non-stationary process
In order to study how non-stationarity might affect the ability of bivariate NeuroXidence to detect differences in JSE frequencies between two conditions, two simulated data sets are generated by the same non-stationary processes. Different processes are used to generate the data in 11 separate periods, each 2 s in length, in order to model some common features that are observed in real data sets, such as random spiking modeled by an inhomogeneous and independent Poisson process (period 1), regular spiking modeled by a g-process (period 2), low rates (period 3, 4), rate modulation (periods 5-10), latency co-variation of rate responses across neurons (periods 8, 10), and trial-by-trial variability of firing rates (period 11; Figure 6). Each of the simulated data sets contains 50 trials of 18 simultaneous spike trains and are used to calculate percentages of false positives of JSEs.
Bivariate NeuroXidence is applied to the simulated data sets with a sliding window of length 800 ms. The null hypothesis (H 0 ) is defined as there being no difference in the occurrence frequencies of JSEs between two conditions for complexities 2 to 6. Based on the results of hypothesis tests for each occurring JS-pattern, we derive the percentage of JS-patterns that show significantly different occurrence frequencies between the two conditions. In order to make comparisons across complexities, the number of significant JS-patterns for each complexity value is normalized by the total number of JS-patterns for the corresponding complexity.
of JSEs between the two conditions, while for higher complexity patterns, the number of surrogates has no effect on the test power. Thus, S = 1 is a reasonable choice for higher complexity JS-patterns.

Sub-and supra-patterns of induced JS-patterns
To study the different occurrence frequencies of sub-and suprapatterns between two conditions, a single-interaction process and an independent and homogenous Poisson process were used to generate two simulated data sets. Each data set contained 50 trials of eight simultaneous spike trains. The JS-patterns induced by the single-interaction process are called mother-patterns. Bivariate NeuroXidence is then applied to evaluate the occurrence frequencies between the two conditions of any JS-patterns caused by mother-patterns.
For mother-patterns with a complexity greater than 2, subpatterns are expected to be detected. However, if the mother-pattern is the only reason for the existence of sub-patterns, then the test power of the sub-patterns should decrease as the complexity of the sub-patterns decreases. Supra-patterns are also induced by a mother-pattern, but these are generated by the occurrence of spikes from additional neurons, which by chance coincide with the mother-pattern. Thus, the occurrence frequency of supra-patterns is at most that of the mother-pattern. Moreover, the test power of supra-patterns is expected to decrease as the complexity of the supra-patterns increases (Figure 5).  5 | Test power of bivariate NeuroXidence for an induced motherpattern and its sub-and supra-patterns. Two simulated data sets were generated as two different conditions. Each sub-figure shows the gray-coded test power of a certain mother-pattern, all sub-patterns of lower complexities, and all supra-patterns of higher complexities. The excess rate of JSEs in one condition, which corresponds to the mother-pattern, is given on the x-axis. The standard parameters were chosen as T = 50 trials, background spike rate r = 15 spikes/s, S = 20 surrogates, and h = 5. (A-D) Shows the variations of mother-patterns with complexity 2-7.

part II: MultIvarIate neuroXIdence
Next we present results demonstrating the robustness and sensitivity of the multivariate version of the new method for testing modulations of spike-spike synchronization across more than two experimental conditions. We start with results demonstrating the robustness of this new multivariate NeuroXidence, by showing false-positive rates estimated for different spike-train models. Lastly, we demonstrate the sensitivity and test power of the presented method.

False positives for stationary processes
Six simulated data sets are generated to represent six different conditions, where each data set contains eight neurons, and each spike train is generated by an independent and homogenous Poisson The simulated spike trains for the two conditions are generated by the same non-stationary processes for all of the time periods, which means that H 0 should not be rejected, or, put another way, any excess or lack of JSEs between the two conditions should stay under the chance level. The actual-false-positive rate, which is the percentage of false rejections of H 0 (Bialek and Rieke, 1992), is clearly below test level (5%) for all complexities and throughout all sliding windows, during all of the periods ( Figure 6D). Therefore, bivariate NeuroXidence is a conservative way to detect the excess or lack of JSEs between two conditions for non-stationary processes. The false-positive rates are not affected by low rates, rate modulation, latency variability, or trial-to-trial rate changes. Moreover, different processes (g-process, Poisson process), which are used to generate the spike trains, do not cause extra false positives to be detected by bivariate NeuroXidence. Two simulated data sets consist of 11 periods, each 2 s in length. Each period is modeled by different features, which are used to generate the spike trains. An inhomogeneous and independent Poisson process was used as a standard model, and three additional features were added as modifications. Feature a (period 11) describes the changing rates across trials and neurons, such that neurons 1-9 were modeled by a homogenous Poisson process with a background rate of 15 spikes/s, while the rates of neurons 10-18 changed from trial-to-trial from 15 to 30 spikes/s. Feature b (periods 8, 10) indicates latency covariations, where the latency for each trial, for all of the neurons, varies by a random amount between 0 and 100 ms. Feature c (periods 2, 4, 7) represents inhomogeneous gamma processes, with shape factor g = 7, instead of Poisson processes. (B) The peristimulus time histogram (PSTH) displays the rate profile of the non-stationary processes. During period 5, the spike rate is modulated from 5 to 50 spikes/s with a Gaussian shape with s = 250 ms, and during periods 6, 7, and 8, s t = 50 ms. The spike rates in periods 9 and 10 were modulated from 5 to 30 spikes/s by a step function. (C) The number of individual JS-patterns of complexities 2-6 that were detected in each sliding window (t c = 5 ms, "SW" = sliding window = 800 ms). (D) The percentage of JS-patterns that show significantly different occurrence frequencies between the two conditions (test level 5%).
100) and the second through sixth mean spiking rates (r 1 = 15, r 2 = r 3 = r 4 = r 5 = r 6 = 7, 10, 30, 60, 90 spikes/s). In both scenarios, the simulated data sets are generated in order to determine the false-positive rate for JS-patterns of complexities 2-6. In both scenarios, none of the results from the model data sets for any of the parameter combinations exceed the chance level, for either a test level of 5% (Figures 7A,B and 8A,B) or 1% (Figures 7C,D and  8C,D) according to both ANOVAs (Figures 7A,C and 8A,C) and Kruskal-Wallis tests (Figures 7B,D and 8B,D). The false-positive rate decreases with the increase of the number of surrogates. As is the case for bivariate NeuroXidence, when the surrogate number is set to 1, multivariate NeuroXidence provides a conservative method for detecting the excess or lack of JSEs among six groups of neurons, no matter whether rate-covariations exist among the groups or not.

Test power for the multivariate case
To assess the test power of multivariate NeuroXidence, one simulated data set is generated by a single-interaction process based on a Poisson process (Fries et al., 2001), while five other data sets are generated by independent and homogenous Poisson processes. Each data set contains eight simultaneous spike trains. For the single-interaction process, correlated spike trains were characterized by a background rate, which corresponds to the independent spiking of neurons, and by a JSE rate, which defines the expected frequency of the JS-pattern of interest.
processes. Multivariate NeuroXidence is then applied to evaluate the false positives due to an excess or lack of JSEs among these six conditions. Two scenarios are applied to generate the rate profiles of the Poisson spike trains for all of the conditions. Scenario I: Six simulated Poisson models exhibit rate covariation, such that the rate profiles of the six simulated data sets are identical. Scenario II: Six simulated Poisson models exhibit partial rate covariation, such that the rate profile of one simulated data set is chosen differently from the other five, which covary.
For both scenarios, multivariate NeuroXidence is applied to six simulated data sets to detect an excess or lack of JSEs among the six sets of model data (Figures 7 and 8). As with the bivariate analysis, the standard parameters for both scenarios include a sliding window with a duration of l = 200 ms, T = 50 trials, a mean spike rate r = 15 spikes/s, and h = 5. Scenario I utilizes S = 20 surrogates samples, and eight different combinations of parameters are derived for the model data by varying, in turn, the number of trials (T = 20, 50, 100) and the mean spiking rates (r 1 = r 2 = r 3 = r 4 = r 5 = r 6 = 7, 10, 30, 60, 90 spikes/s). Scenario II uses S = 1 surrogate sample, and eight different combinations of parameters are derived by varying, in turn, the number of trials (T = 20, 50, lack of JSEs across different experimental or behavioral factors. This new variation of NeuroXidence considers the full auto-structure of each spike train and the trial-to-trial variability; moreover, these extensions correct for the effects of rate-covariation among groups of neurons, low firing rates, fast rate modulation, regular spiking defined by a g-process, and latency co-variation of rate responses across neurons. After calibrating bi-and multi-variate NeuroXidence on simulated data sets, the new extension reliably detects an excess or lack of JSEs due to different experimental or behavioral factors, for both stationary and non-stationary processes. Moreover, this extension to NeuroXidence maintains a high sensitivity for detecting different modulations of the JSEs. The procedure for applying bi-and multi-variate NeuroXidence consists of two main parts: detection of the JSEs, whose occurrence frequencies are above the chance level; and comparison of the frequency of JSEs across different experimental/behavioral factors. The latter finds the relationship between the strength of synchronous firing and the neuronal system state. In order to achieve this goal, different statistical tests are applied to the simulated data sets at each step of the procedure. The statistical test applied during the detection step is described in the original NeuroXidence article (Pipa et al., 2008), which focuses on the distinction between the chance level of coordinate firing and the underlying synchronous firing events. Because different features of spike trains, such as the auto-structure The standard parameter set is chosen as 50 trials (T), a background spike rate of 15 spikes/s (r), and h equals 5. Six parameter sets are used to study how the test power of multivariate NeuroXidence is affected by the number of trials (T = 20, 50, 100) and the number of surrogates (S = 1, 20). The test power is derived for JS-patterns of complexity 2-6, with a test level set to 5%. It is clear in Figure 9 that an increase in the number of trials leads to an increase of the test power, as exhibited by the leftward shift and increase in steepness of the curves. When the number of surrogates changes from 20 to 1, the detection of extra JSEs for either higher or lower complexity values shows only subtle changes, which means multivariate NeuroXidence is sensitive enough to detect small differences in the JSE frequencies for each JS-pattern, among several conditions. Therefore, the number of surrogates can be set to 1 without harming the test power of multivariate NeuroXidence.

dIscussIon
NeuroXidence reliably and robustly detects synchronous firing patterns, beyond those expected by chance, for the univariate case (Pipa et al., 2008). However, the originally proposed procedure is not suitable for comparing the modulation of synchronous firing strength across experimental or behavioral factors, nor is it robust against features that might induce false positives to such an estimate. Here we have extended NeuroXidence to analyze bi-and multi-variate cases, which enables a comparison of an excess or surrogates need to be generated to estimate the frequencies of JSEs that are expected by chance. This parameter is the most important one, when the mean firing rate is not the same for different conditions. Based on equation (1), surrogate data sets are generated by jittering the original spike trains to estimate the chance level of frequencies of JSEs [expressed as 1 ]. The amount of jitter applied to the original spike train is defined by t r . However, the distribution of the frequencies of JSEs in the S surrogate data sets is not symmetrical. The skewness of the distribution, which will decrease as the mean firing rate increase, leads to a biased estimate of the underlying synchronous firing events. Therefore, by setting the number of surrogates to 1, both bi-and multi-variate NeuroXidence are able to maintain a high test power and low false-positive rate, thereby ensuring a non-biased estimate of the JSEs for each factor.
Based on the test power of simulated data sets, an average of 2 JSEs/s difference in the occurrence frequencies of JSEs, among several conditions, can be detected with bi-and multi-variate of each spike train, the rate-covariation among groups of neurons, the low firing rate, etc., might effect an estimate of the rate of coordinated firing events, the chance level for coordinated firing is estimated directly from the surrogate data sets. The statistical test in the comparison step focuses on the different frequencies of the JSEs, which result from different experimental/behavioral factors. In order to evaluate the differences, we compare the occurrence frequency of the same JS-pattern k, for each trial t, across the different factors j. Different statistical tests (t-tests, Mann-Whitney U tests, ANOVAs, and Kruskal-Wallis one-way tests) are applied to evaluate whether coordinated firing events are the same across conditions or not.
There are three important parameters for bi-and multi-variate NeuroXidence: the time scale t c , the time scale t r , and the surrogate number S. The first time scale, t c , defines the precision of synchronous events to be within the 1-to 10-ms range. The second time scale, t r , establishes the lower bound of rate covariations, and it is defined to be h times slower than t c . The parameter S defines how many across different experimental or behavioral factors. After being calibrated on simulated data sets, the extension exhibits its reliable and robust detection of synchronous firing patterns across different conditions, for both stationary and non-stationary processes. Thus, bi-and multi-variate NeuroXidence provides a new approach for evaluating different strengths of synchronous firing across experimental data sets.
NeuroXidence. Given more trials and more conditions, the test power surpasses 90% when the JSE rate is as low as 0.6 JSEs/s. In general, the conservative and sensitive properties of bi-and multivariate NeuroXidence ensure that any differences in the modulation of synchronous firing strength, among several experimental or behavioral factors, can be reliably detected.

conclusIon
Expanding upon the previously published univariate NeuroXidence method (Pipa et al., 2008), here we present bi-and multi-variate NeuroXidence, which is designed to detect the lack or excess of JSEs