Establishing a Statistical Link between Network Oscillations and Neural Synchrony

Pairs of active neurons frequently fire action potentials or “spikes” nearly synchronously (i.e., within 5 ms of each other). This spike synchrony may occur by chance, based solely on the neurons’ fluctuating firing patterns, or it may occur too frequently to be explicable by chance alone. When spike synchrony above chances levels is present, it may subserve computation for a specific cognitive process, or it could be an irrelevant byproduct of such computation. Either way, spike synchrony is a feature of neural data that should be explained. A point process regression framework has been developed previously for this purpose, using generalized linear models (GLMs). In this framework, the observed number of synchronous spikes is compared to the number predicted by chance under varying assumptions about the factors that affect each of the individual neuron’s firing-rate functions. An important possible source of spike synchrony is network-wide oscillations, which may provide an essential mechanism of network information flow. To establish the statistical link between spike synchrony and network-wide oscillations, we have integrated oscillatory field potentials into our point process regression framework. We first extended a previously-published model of spike-field association and showed that we could recover phase relationships between oscillatory field potentials and firing rates. We then used this new framework to demonstrate the statistical relationship between oscillatory field potentials and spike synchrony in: 1) simulated neurons, 2) in vitro recordings of hippocampal CA1 pyramidal cells, and 3) in vivo recordings of neocortical V4 neurons. Our results provide a rigorous method for establishing a statistical link between network oscillations and neural synchrony.


Introduction
A leading theory of current neuroscience is that synchronous firing of neurons driven by network-wide oscillations may encode and transmit information within and across brain regions [1][2][3][4][5][6][7][8][9]. Supporting this theory, a number of studies have suggested that synchronous firing of action potentials or "spikes" may indeed occur in conjunction with oscillations in local field potential (LFP) [10][11][12][13][14]. However, a missing link in this theory has been the ability to dissociate enhanced spike synchrony due to network-wide oscillations from enhanced spike synchrony that may be due to other measured or unmeasured sources. Recently, we developed a statistical framework in which the association between spike synchrony and measured covariates may be assessed [15,16]. Here we show how this approach may be applied to describe the relationship between spike synchrony and oscillatory activity.
Using point process regression models, which take the form of generalized linear models (GLMs), our statistical framework compares the observed number of synchronous spikes within a small time window (here, 5 ms) to the number predicted by chance, under varying assumptions about the factors that affect the firing of each individual neuron [15,16]. The number of synchronous spikes predicted "by chance" refers here to the number predicted under conditional independence after conditioning on the various measured factors that have been hypothesized to affect individual-neuron spiking. For example, two neurons having fluctuating stimulus-driven firing rates will produce some number of synchronous spikes even if they are acting independently. The point process regression method fits fluctuating firing rate functions for each neuron separately, then predicts the number of synchronous spikes under conditional independence (i.e., after conditioning on these fluctuating firing rates), and compares the prediction to the observed number of synchronous spikes. In this way, a single factor may be either included or excluded from the regression model in order to quantify that factor's ability to explain the observed spike synchrony.
In this article, we consider the contribution of network-wide oscillations by comparing observed and predicted spike synchrony after conditioning on the phase of an LFP representing a network-wide oscillation. Thus, we predict spike synchrony with and without inclusion of LFP phase as an explanatory variable for each neuron separately. To demonstrate that increased spike synchrony is associated with a network-wide oscillation, we would begin by establishing that, without considering LFP phase, the observed number of synchronous spikes is greater than the predicted number by a statistically significant magnitude, after conditioning on both stimulus-driven firing rates and recent post-spike history effects. This would indicate a failure of the phase-free model to accurately account for spike synchrony. We would then include the LFP phase in the model, and if it succeeded in predicting spike synchrony, then we would conclude that LFP phase can explain the remaining spike synchrony. Furthermore, we could estimate the proportion of excess synchronous spikes accounted for by the LFP phase.
The same procedure could be used, instead to demonstrate the role of network-wide oscillations in suppressing spike synchrony.
In order to carry out this general procedure, we first need to model an individual neuron's spiking probability in terms of LFP phase. We follow [17], which recently described and assessed point process regression models that include a sinusoidal phase term. We enhance their approach by weakening the sinusoidal assumption, allowing the phase relationship to be nonparametric as in [18], and we add to the favorable results of [17] by showing that, in estimating phase relationships, the point process regression model can reduce bias and meansquared error in comparison with the more familiar spike phase histogram approach. Using this point process regression model, we are then able to quantify the dependence of synchronous spiking on an oscillatory modulation. We illustrate the method using simulated neurons, in vitro recordings of hippocampal CA1 pyramidal cells, and in vivo recordings of neocortical V4 neurons from a behaving monkey.

Point Process Model for Spike Trains
We assume that the spiking of each neuron follows a point process and, following [19] (page 592), we write its conditional intensity function as λ(tjH t ,X t ), where H t represents the spike history (auto-history), and the covariate X t represents other external factors. In this work, we let X t include the stimulus and the LFP phase, denoted by X t = (S t ,F t ). We assume the conditional intensity takes a multiplicative form, which becomes additive on the log scale: where t Ã is the last spike time preceding t (see Materials and Methods). We use splines to capture stimulus and auto-history effects, and circular splines to capture LFP phase effects. Our point process model thus takes the form of a standard generalized linear model (GLM). We also ensure identifiability by imposing a set of restrictions (Eqs (20) and (21)), which are implemented within a maximum likelihood estimation (MLE) algorithm. The parametric bootstrap is used for acquiring 95% confidence bands.
To illustrate the ability of the MLE algorithm to recover the model in Eq (1), we simulated 100 spike trains ( Fig 1A) with known functions λ 1 (t), λ 2 (t − t Ã ) and λ 3 (ϕ). Using the simulated spike trains and phase of the oscillatory drive (representing a network-wide oscillation), the MLE algorithm accurately fit the underlying spike history (Fig 1B), stimulus ( Fig 1C) and phase modulation (Fig 1D) effects. Our approach can thus accurately recover the statistical relationships between firing rate and various external factors.
The model in Eq (1) is a "full" model including stimulus, auto-history, and an oscillatory factor. Importantly, we can remove selected factors from the full model (e.g., the LFP phase modulation) and still fit the spike trains using the same procedure. Indeed, in the following results, we also fit a simplified model lacking the oscillatory factor, log lðtjH t ; X t Þ ¼ log l 1 ðtÞ þ log l 2 ðt À t Ã Þ: ð2Þ

Estimation of LFP Phase Modulation
Many researchers have reported that firing rate is modulated by the phase of specific networkwide oscillations in different brain areas, such as monkey V1 [20], rat hippocampus [21], rat prefrontal cortex [22], mouse olfactory bulb [23], human pedunculopontine nucleus [24], lamprey reticulospinal neuron [25], and so on. Almost all of these results [20-22, 24, 25] used spike phase histograms to show how firing rate is modulated by the oscillation. The significance of phase locking can be evaluated using Rayleigh's Z statistic [22]. The model in Eq (1) offers an alternative method of computing LFP phase modulation. We simulated N spike trains and estimated LFP phase modulation using two different methods: 1) the classical spike phase histogram, and 2) by fitting λ 3 (ϕ) with the point process regression of model (1). The true LFP phase modulation function is defined as λ 3 (ϕ). The discrepancy between the estimatedl 3 ðÞ and λ 3 (ϕ) is measured by the integrated squared Using each method, we can derive point-by-point standard errors and 95% confidence bands for the LFP phase modulation. The mean integrated squared error (MISE) is then defined as: where n is the total number of data sets and i is the index of ith data set.l i 3 ðÞ is computed given N repeated trials of spike train in ith data set. We can decompose MISE in terms of the sample mean l 3 ðÞ in the form of: which provides an estimator of variance plus bias squared. The histogram method is highly dependent on the bin size for smoothing. We picked the optimal bin size that minimizes the MISE. Fig 2C illustrates how the MISEs of the two methods are dependent on number of trials N. Both methods achieve smaller MISEs when more data are used, but the spike phase histogram method consistently exhibits a much larger MISE than the GLM method. Indeed, the spike phase histogram MISE reaches an asymptote for high N that is much larger than the MISE of the GLM method. In Fig 2F we show the variance and bias separately for the two methods. These results show that the spike phase histogram method retains a large bias, explaining the MISE asymptote in Fig 2C. The LFP phase modulation estimated by the spike phase histogram method additionally exhibits significantly larger variance than the GLM method for small sample sizes (< 17 trials).
Two additional comments can be made about the results shown above (Fig 2). First, when few trials or samples are available, only the GLM method can provide an accurate estimation of the LFP phase modulation of a neuron's firing. Second, for moderately large samples, the error in the estimation of the LFP phase modulation by the spike phase histogram method arises primarily from estimation bias. We can explain this second point by considering the definitions of the two methods. The term λ 3 (ϕ) describes how an oscillation changes the firing rate and is independent of other factors (stimulus, auto-history, etc.). In contrast, the spike phase histogram method provides the distribution of phases when a spike occurs, denoted as P data (ϕ). Since the generation of a spike train is influenced by factors other than the oscillation, especially for a non-Poisson process, P data (ϕ) is conceptually different than λ 3 (ϕ). Below, we explore this conceptual difference further.
Bias in the estimation of the LFP phase modulation using the spike phase histogram method. To simplify our GLM, let us assume that the stimulus effect is a constant λ 1 (t) = C and f is the frequency of the oscillation. Now the firing probability model is: Suppose we observed a spike train {u k }. The spike phase histogram method provides an estimate of the distribution of {ϕ u k }, where P data (ϕ u k = ϕjϕ u k − 1 = ϕ 0 ) is the conditional probability of ϕ u k = ϕ given the phase of its previous spike is ϕ 0 . This conditional probability can be computed using the distribution of the waiting times [19] (page 602), If we have a spike at phase ϕ, while its previous spike is at phase ϕ 0 , then the waiting time should be within the set W ¼ fDu : Eq (3) shows that P data (ϕ) is an eigenfunction of P data (ϕ u k = ϕjϕ u k − 1 = ϕ u k − 1 ). This is very hard to compute analytically, but we can get its numerical solution instead. When we discretize ϕ and write P data (ϕ) as a vector P, Eq (3) can be rewritten as where P 2 R m × 1 , m is the number of bins to discretize ϕ 2 [−π,π), and A2R m × m is the transition probability Thus P is the eigenvector of A and its corresponding eigenvalue is 1. Numerical tools were used to compute P, which is the discrete version of P data (ϕ). In this way we can theoretically determine the LFP phase modulation given by the spike phase histogram. This theoretical prediction accurately predicts the LFP phase modulation estimated from simulated spike trains using the spike phase histogram method (Fig 3B, 3C, 3E and 3F). Using this theoretical prediction, we examined how bias emerges in spike phase histogram estimates of the LFP modulation. We considered Poisson and non-Poisson firing across different mean firing rates. When a neuron's firing approximates a Poisson process (i.e., λ 2 (t) = 1) (Fig 3A), the results of the spike phase histogram match λ 3 (ϕ) independent of the mean firing rate C (Fig 3B and 3C). Indeed, it can be shown that P data ðÞ ¼ l 3 ðÞ 2p is a solution to Eq (3). Specifically, Eq (4) inserted into Eq (3) with some basic substitutions yields: When λ 2 (t) = 1, we replace P data (ϕ) with l 3 ðÞ 2p in Eq (5). Then the integrand in the right side is , which is a probability density function and hence has the integral of 1, making the right side λ 3 (ϕ). Thus, λ 3 (ϕ) is one solution of Eq (3). These results show that the spike phase histogram estimate of the LFP phase modulation is accurate for Poisson firing. In contrast, for non-Poisson firing (in which the neuron's firing rate is influenced by its firing history) λ 2 (t) is no longer a constant ( Fig 3D). As a result, estimates of the LFP phase modulation curve diverge from λ 3 (ϕ) in a firing rate-dependent manner (Fig 3E and 3F). Thus, the GLM method for estimating a neuron's LFP phase modulation is more accurate than the spike phase histogram method for smaller sample sizes (Fig 2) and for non-Poisson firing (Fig 3).

Comparison with Spike Field Coherence
Spike field coherence (SFC) is commonly used to report interactions between spikes and specific oscillations in LFP. Lepage et al. [17,26] showed that SFC is dependent on the expected rate of spiking, and they proposed to use intensity field coherence, which is a rate-independent measure, for inference of spike field synchrony. They also used GLMs to estimate spike field association [17]. In their work, they assumed that the LFP phase modulation is a sinusoidal function with period of 2π, which might not be accurate enough in some cases [22,23]. In our model, to approximate this periodic function we use circular splines [18], which remain easy to fit while being more flexible than a sinusoidal function.
Here, we provide two examples showing that when estimating spike field relationships, the SFC can be misleading. First, we simulated spike trains with three different mean firing rates, then computed SFCs with GNU software Chronux [27]. Fig 4A shows that the three SFCs are different even though they were generated by the same λ 3 (ϕ) = 1 + 0.4cos(ϕ + π). On the other hand, when we use our model to fit the LFP phase modulation functions, Fig 4C shows that there is no difference in phase modulation strength in these three cases. Second, we show that two neurons exhibiting different LFP phase modulations can have the same SFC ( Fig 4B) because they have different firing rates. Again, we can use our model to distinguish these two conditions by their respective LFP phase modulation curves (Fig 4D).

Synchrony and Oscillatory Phase
We now use point process regression of our GLMs Eqs (1) and (2) to analyze the contribution of network-wide oscillations to the synchronous spiking of two neurons. We first present numerical simulation results where ground truth is known and then apply the same technique to experimental neural recordings. Simulation results. In the Introduction, we described how GLMs can be used to assess the role of some potentially relevant factors in modulating spike synchrony. We designed a scenario in which we tested the contribution of a network-wide oscillation (i.e., an oscillatory LFP) to the number of synchronous spikes observed. This scenario is illustrated schematically in Fig 5 for two neurons. The stimulus effects (i.e., the tuning) of the two neurons are different, and both neurons' spiking activities are influenced by their own recent spike histories. Critically, these two neurons also receive a common oscillatory signal with phase F t that modulates their firing rate, but their individual phase modulation curves are shifted (i.e., they have different preferred phases F pref ). Because the preferred phase modulates the average timing of each spike in one oscillatory cycle (in this example, * 10 ms), differences in preferred phase lead to a relative shift in spike timing between the two neurons. The larger this shift, the less synchronized are their spikes. As a result, the observed number of synchronized spikes is dependent on the difference of preferred phase ΔF pref .
This simple scenario was used to demonstrate the effectiveness of the procedure, in principle, and to investigate its statistical power. The assumption that two neurons have different phase modulation curves has been reported both experimentally [20,28] and theoretically [29]. Jia et al. [20] have shown that neurons in area V1 have various preferred phases and the distribution of the preferred phase can change in response to different stimuli. Richardson [29] computed analytically the modulation of the oscillatory signal for an exponential integrated-andfire neuron. He showed that the modulation is influenced by biophysical properties of the neuron. He also showed that there is a phase lag between the peak firing rate and the peak of the oscillatory signal, which corresponds to the preferred phase in λ 3 (ϕ), and this phase lag is dependent on properties of the neuron. Usually the oscillations near two neurons in a small area are very similar; thus the assumption that two neurons receive a common modulation is reasonable. For two neurons located far apart (e.g., two brain areas), this assumption should be useful as long as the two oscillations are coherent. This more general case is relevant to hypotheses about mechanisms of neural communication [4,20].
To demonstrate directly the relationship between an oscillatory LFP and spike synchrony, we simulated spike trains from two neurons, then fitted models (1) and (2). For each model we Observed number of synchronized spikes Predicted number of synchronized spikes of its theoretical counterpart z 12 defined in [16]. Under conditional independence, we have log z 12 = 0, while conditional dependence yields either excess synchrony (log z 12 > 0) or suppressed synchrony (log z 12 < 0). We tested H 0 : logz 12 = 0 using a parametric bootstrap (see Materials and Methods). Results are shown in Fig 6. Using model (2) (i.e., without the oscillatory factor) we found that logẑ 12 is dependent on ΔF pref . This is because the relative phase preference of the two neurons changed the observed number of synchronized spikes, while the predicted number is almost the same when the contribution of the oscillatory LFP is disregarded. In contrast, when we included the oscillatory factor according to Eq (1), we found logẑ 12 to be close to 0 and independent of ΔF pref . Thus, including the oscillatory factor in our model removes the apparent conditional dependence of the predicted spike synchrony on the relative phase preference of the two neurons, and we can conclude that spike synchrony is associated with the oscillatory phase.
We picked two different values of ΔF pref (purple and cyan arrows in the Fig 6A) to demonstrate the described hypothesis test. In the first example, we obtained evidence against the null hypothesis of log z 12 = 0 using the simplified model ( Fig 6B). That is, there is evidence that the two neurons are not conditionally independent given only the stimulus effects and spike history effects: they exhibit significant levels of excess spike synchrony. Fig 6C shows that including the oscillatory factor accounts for this excess synchrony. In other words, consideration of the oscillatory LFP can explain the higher than expected levels of spike synchrony predicted by the stimulus and spike history effects alone. In turn, lower than expected levels of spike synchrony predicted by the stimulus and spike history effects alone can be explained by consideration of the oscillatory LFP in the second example, in which the oscillatory LFP suppresses spike synchrony (Fig 6D and 6E).
We also investigated the amount of data needed to reliably detect excess synchrony by generating spike trains with varying numbers of trials, varying values of z, and two levels of firing rate, and then computing the probability of rejecting the null hypothesis (i.e., the statistical power). Fig 6F displays the power when we used the same simulation parameters (apart from z and number of trials) as in Fig 6A-6E. A standard target for power in the statistics literature is 0.8, and we have indicated this level of power with a red line in Fig 6F. Thus, to attain this high level of power when z = 1.125 we need 70 trials, but when z = 1.4 we need only 5 trials. This number is also highly dependent on the mean firing rate. When we change the firing rate from 25 Hz to 10 Hz, we need much more data to detect excess synchrony (Fig 6G). The simulation procedure is computationally slow, but a fast approximation is given by where T is the length of one trial, λ 1 and λ 2 are mean firing rates of two neurons, and δ is the bin size for detecting synchronized spikes. The approximate power from this formula is given by the green curves in Fig 6F and 6G. The formula is derived in Materials and Methods.

Applications to Experimental Neural Recordings
To further demonstrate the value of our approach, we next examined the relationship between an oscillatory signal and spike synchrony in experimental neural recordings from two distinct preparations: hippocampal CA1 pyramidal cells recorded in vitro and V4 neurons recorded in vivo.
Hippocampal CA1 pyramidal cells. We first designed an experiment to resemble the scenario proposed in Fig 5 using whole-cell patch clamp recordings in a controllable acute slice preparation. In this experiment, we recorded the spiking response of, and spike synchrony between, two CA1 pyramidal cells (Fig 7A and 7B) in response to an arbitrary stimulus with and without a shared oscillatory signal. Critically, to directly test the relationship between the oscillatory signal and the resulting spike synchrony, we limited potential confounding influences on spike synchrony (e.g., common neuromodulatory influences, coupling between two neurons) by recording these neurons sequentially in two separate slices. Each neuron was injected with 100 trials of a 2 s-long 150 pA step current overlaid with a slow sinusoidal current (2 Hz frequency, 25 pA amplitude) and white noise (σ = 10 pA) to evoke physiological spike trains with low trial-to-trial reliability. The slow 2 Hz component was the same for all trials, and it generated visible time-varying fluctuations that are visible in the raster plots in Fig 7A  and 7B, which led to a time-varying PSTH that was captured by the λ 1 (t) term in Eq (1). On 50 random trials ("Exp. 2"), an additional sinusoidal current (40 Hz frequency, 15 pA amplitude) with random initial phase (but identical between the two neurons) was also injected to simulate a gamma frequency network-wide oscillation. The 40 Hz component was not consistent over trials due to varying initial phases (S1A Fig), and its effect is, therefore, not captured by λ 1 (t). Instead, this 40 Hz modulatory effect was captured through the term λ 3 (F t ) in Eq (1).
Thus, in the 50 trials without the simulated network-wide oscillation ("Exp. 1"), each neuron fired according to the its own stimulus and auto-history effects, generating a certain level of largely spontaneous spike synchrony reflecting the neurons' fluctuating stimulus-driven firing rates. Using our simplified model (Eq (2)), we fit the spike trains from "Exp. 1" and predicted the number of synchronous spikes. As expected, the observed and predicted number of synchronous spikes closely matched (Fig 7D), consistent with the two neurons being conditionally independent given the arbitrary stimulus waveform and their own recent spiking histories. That is, no other factors were necessary to explain the observed number of synchronous spikes. However, using our simplified model to fit the spike trains from "Exp. 2" (Fig 7A and  7B), we observed a significantly greater number of synchronous spikes than could be explained by the stimulus and the neurons' spike histories alone (Fig 7E). This conditional dependence between the two neurons arose because the firing of the two neurons was modulated by the simulated network-wide oscillation (Fig 7C). Indeed, using our full model (Eq (1)) to fit the spike trains from "Exp. 2" (S1B, S1C and S1D Fig), the number of synchronous spikes observed closely matched the number of synchronous spikes predicted (Fig 7F). This experiment demonstrates that when two experimentally recorded neurons are not modulated by a shared oscillatory signal, then the simplified model (Eq (2)) can account for the observed number of synchronous spikes. However, when two neurons are modulated by a shared oscillatory signal (such as an oscillatory LFP, reflecting a network-wide oscillation), then a model including this oscillatory factor (Eq (1)) is necessary to account for the observed number of synchronous spikes. In contrast with our simulation above, the firing of these CA1 neurons is not described by the GLM in Eq (1) exactly. This model mismatch did not restrict the application of our method.
V4 neurons. In this experiment, spike trains from a pair of neural units in V4 were simultaneously recorded (Fig 8A and 8D) with a multi-electrode array during a fixation task in which spontaneous activity was measured. These data have been analyzed in another paper [28], which examined the relationship between individual neuron's activity and large-scale network state. Here we wanted to test whether network-wide oscillations contribute to the excess pairwise synchrony. For each neuron, we defined its surrounding LFP as the average of LFPs recorded at its adjacent electrodes (see Materials and Methods). The spike-triggered average of the LFP for two neurons showed that the two neurons are phase locked to their surrounding field potential (S2C and S2D Fig and [28]). We also found that the LFP showed a prominent slow oscillation (Fig 8B and 8E). The LFP is thought to be the integrated effect of synaptic and spiking activity [30] near the recording sites. We filtered the LFP on each electrode within the 4-25 Hz band and extracted its phase to fit our full model (Eq (1)). Using the same procedure as in the case of the hippocampal CA1 pyramidal cells, we found that a significantly larger number of synchronous spikes were observed than could be explained by the simplified model (Fig 8C), while the full model fully explained the spike synchronization observed between the two neurons ( Fig 8F). These results show that for these two neurons in vivo, spike synchronization is associated with the network-wide oscillation.

Discussion
In this paper, we have shown how the GLM methods of [15][16][17]26] may be combined in order to assess the potential contribution of network-wide oscillations to neural synchrony. The novel approach presented in this study complements existing alternatives [31][32][33] by: introducing models of single neuron firing based on stimulus-related fluctuations as well as a networkwide oscillatory signal; using those models to make predictions about spike synchronization; and quantifying departures from those predictions in the observed data. We demonstrated the advantages of this novel approach using both neural simulations and experimental neural recordings in vitro and in vivo.
In our analyses, we have utilized a repeated-trial structure, which allowed us to estimate the stimulus effects as a function of time, λ 1 (t). We note, however, that the same approach could be applied using a linear response filter [34][35][36] or analogous nonlinear methods. Previous work has shown the close relationship between GLM neurons and integrate-and-fire neurons [37][38][39]. We only considered one band of oscillation in simulation and experimental examples, but it is straightforward to extend this method to the case of multiple oscillations by including additional terms in the model of Eq (1). Sometimes the firing probability may be related to the amplitude of the oscillation A t , or the magnitude of an LFP B t (cf. [17]). If so, we can change Overall, the key step of this method is to build an approximately correct GLM. The specific form of GLM depends on the data and we can check model performance using time rescaling [40]. We have also included a simulation to show that even when the model is mis-specified, and therefore less sensitive, it can detect spike-LFP relationships (S3 Fig). We have also defined spike synchrony to involve the firing of two neurons within a few milliseconds of each other (i.e., with zero lag on average). In other contexts, however, interest may focus on two neurons firing in procession with a consistent positive or negative lag of many milliseconds. Our approach could be easily applied to such lagged-synchrony cases as well.
In this paper, we consider only pairwise synchrony. By combining our approach with the procedure proposed by [16], we can also test the role of oscillations in three-way synchrony. Briefly, we fit all single neuron firing probabilities and then compute the pairwise synchrony coefficientŝ z ij ; we can then use an iterative algorithm to estimate the three-way synchrony coefficientẑ ijk , and to test the null hypothesis of two-way interactions, instead of three-way interaction. In principle the same steps may be followed for more than three neurons, but simulations in [16] show that very large data sets would be needed in order to demonstrate higher-order interactions convincingly in the absence of stronger assumptions about the nature of those interactions.
It has been argued that synchronous firing resulting from network-wide oscillations could provide an essential mechanism of network information flow, and further serve as a a marker distinguishing normal from diseased states (e.g., see [41][42][43][44][45][46][47][48]). On the other hand, there has been considerable debate on this subject (see [49] and references therein). We remain agnostic on this, and importantly, the value of our methods does not depend on the ultimate outcome of this debate. Instead, we view synchrony, more descriptively, as a feature of spike train data that needs to be explained. To this end, the framework that we have introduced here is useful for quantifying the extent to which oscillations, as a feature of neural activity, are associated with synchronous spiking among neurons. Armed with this method, future experiments can measure oscillations and synchrony in a statistical framework in which their contributions to cognitive and behavioral processes can be accurately quantified.

Point-Process Framework
In a continuous time interval (0,T], a neuron can fire a spike at any discrete time point u i . A series of spikes {u i } for 1 i N forms the spike train, where 0 u 1 < Á Á Á < u N T. We take the spike train to be a point process, which is characterized by its conditional intensity function where N(t) is the total number of spikes prior to time t, H t is neuron's own spiking history prior to time t, and X t includes all other relevant covariates. When Δ is small, λ(tjH t ,X t ) Á Δ approximates to the firing probability in the time interval (t,t + Δ). To determine how different factors contribute to firing rate we write λ(tjH t ,X t ) as a function of (H t ,X t ) We can include different factors into this model and study their effects. Usually the stimulus S(t) is included when neurons show selectivity to stimuli. In this work, because we are interested in phase modulation by an oscillatory signal, the phase of the specific oscillation F(t) is also included.

Generalized Linear Model and Maximum Likelihood Estimation
To take advantage of the generalized linear model (GLM) framework we divide T into K equally spaced bins, thus taking the bin width to be Δ = T/K. Δ is small enough that no more than one spike event in each bin, e.g. Δ = 1 ms. Therefore the probability of observing one spike in kth bin is Using the vector Y & R K × 1 to represent the spike train {u i }, y k is the number of spikes in kth bin. Since we choose small bin width Δ, y k is not bigger than 1, i.e. y k 2 {0,1} and, from the Poisson approximation to the binomial for small p we take the probability of observing y k given H t k and X t k to be where t k = kΔ. The loglikelihood function is and this is maximized to determine the MLE fit.
We assume that log[λ(tjH t ,X t )] can be written as a sum of specific functions of each covariate. Here we are studying three factors, stimulus, recent post-spike auto-history, and oscillatory phase, and we write log ½lðtjH t ; X t Þ ¼ f 1 ðstimulusÞ þ f 2 ðauto À historyÞ þ f 3 ðoscillationÞ ð 11Þ Here, S(t) is a possibly time-varying stimulus and f 1 (stimulus) determines the trial-independent time-varying firing rate, i.e., the effect that is usually associated with the peri-stimulus time histogram (the PSTH), which may be estimated due to the repeated trial structure of the experiment. The recent post-spike auto-history effect is assumed here to be dominated by effects subsequent to the most recent spike t Ã prior to time t, as in [50], so we assume f 2 (autohistory) has the form f 2 (t − t Ã ). The oscillatory term f 3 (oscillation) is defined as f 3 (F t ), where F t is the phase of specific oscillation. In summary, our spike train model has the form and we will assume f 1 (Á), f 2 (Á)and f 3 (Á) are smooth functions.

Approximate Function with Spline Basis
To fit the smooth functions f 1 (Á), f 2 (Á)and f 3 (Á) we use cubic splines of the form where {a i (t)} is a B-spline basis set for f 1 (t) within the range t 2 (0,T], fb j ðt À u Ã t Þg is a B-spline basis set for f 2 (t − t Ã ), and {r k (ϕ)} is circular spline basis set for f 3 (F t ). Thus, we use maximum likelihood to fit the coefficients Θ = {α,β,γ}. We used open source software FDAfuns [51] to create each B-spline basis sets after manually selecting knots. For the circular spline we pick knots equally spaced in [−π,π]. Once we get all knots {ϕ i }, acquiring the related basis function is straightforward [18] using 4 cos ð2pmð À k ÞÞ: ð16Þ In numerical implementations, we usually cut the summation from m = 1 to m = 4 because amplitude of each term decreases quickly.

Maximum Likelihood via Iteratively Re-weighted Least Squares
Because L in Eq (10) is a concave function, we can use iteratively reweighted least squares (IRLS), as in typical GLM implementations. From Eqs (8), (10), and (15), we can rewrite loglikelihood in matrix form Here we have three parameter sets to fit {α,β,γ}. If we fit all three parameter sets together, the dimension space of this GLM model is relative large. To make model fitting efficient, we prefer back-fitting, i.e., fitting each parameter set separately, and iterating cyclically. For example, when we fit the parameters {α}, we hold the parameters {β,γ} constant and rewrite Eq (18) as where θ 2 {α,β,γ} and V is the corresponding covariate matrix. We fit {α,β,γ} in a sequence and then iterate the loop until convergence. We also must place identifiability restrictions on {β,γ} because both the auto-history and oscillatory effects modulate the spike trains and the parameters must be constrained to provide unique solutions. We use the constraints To avoid over-fitting of the model, we also add an l 2 penalty into the objective function. Now the problem becomes minimizing objective function Because the objective function Q is convex, we can iteratively maximize Θ by following the updating rule where H is the Hessian of Q and rQ is the gradient of the function, which are obtained as where W is a diagonal matrix The algorithm is summarized as Algorithm 1, shown below.

Spike Synchronization
For a pair of neurons labeled 1 and 2, we fit conditional firing rate for each of them to getl 1 ðt j H t ; X t Þ andl 2 ðt j H t ; X t Þ. Then we can predict the number of synchronized spikes given temporal bins with Δ = 5ms, as in [15], using Given spike trains from these two neurons, we can also get the observed number of synchronized spikes N obs by counting. If a pair of spikes from two neurons has an time interval less than 5 ms, then this pair is counted as a synchronized spike. Synchrony is measured by taking the ratio of these two numbersẑ or logẑẑ When two neurons are conditional independent, Eq (27) can make relative good predictions andẑ % 1 (or logẑ % 0).

Bootstrap Method
Once we have logẑ, we also need to determine its standard error and confidence interval. Furthermore, a p value is required to test the hypothesis logẑ ¼ 0. We use a parametric bootstrap method for these purposes, as in [16]. For example, givenl 1 ðt j H t ; X t Þ andl 2 ðt j H t ; X t Þ we can obtain the p value as follows: • For i in 1: G do 1. Simulate each of the two sets of spike trains, across the same number of trials as in the data, using the respective spike train models withl 1 andl 2 .
2. Compute log z i from the spike trains generated in step 1.
• Compute the number of values {log z i } (out of a total of G such values) for which j log z i f gj >j logẑ j and divide by G. This is the p value.

Power Analysis
Statistical power is the probability of correctly rejecting the null hypothesis when it is false. We used the GLM model in Eq (1) to study power as a function of z and N (N being the number of trials). We simulated N trials of spike train data for each of two neurons, independently, using Eq (1) with intensity functions λ (1) (tjH t ,X t ) for the first neuron and λ (2) (tjH t ,X t ) for the second. The synchronous spikes in the resulting spike trains occur with probability corresponding to z = 1 (independence). In order to obtain sets of spike trains for other values of z we removed all the synchronous spikes from the N simulated spike trains and replaced them with synchronous spikes generated from an intensity function z Á λ 1 (tjH t ,X t ) Á λ 2 (tjH t ,X t ), i.e., for each time bin of width δ, synchronous spikes occurred with probability z Á λ 1 (tjH t ,X t ) Á λ 2 (tjH t ,X t )δ 2 . However, while this is the desired probability of synchronous spikes, it leaves the wrong marginal probability of spiking for each neuron. To adjust these we consider the spike trains made up of only the non-synchronous spikes, and we thin these with probabilities p (j) (t) given by p ðjÞ ðtÞ ¼ l ðjÞ ðtjH t ; X t Þ À z Á l ð1Þ ðtjH t ; X t Þ Á l ð2Þ ðtjH t ; X t Þd l ðjÞ ðtjH t ; X t Þ À l ð1Þ ðtjH t ; X t Þ Á l ð2Þ ðtjH t ; X t Þd for j = 1,2. Note that when we multiply the numerator and denominator of this expression by δ we have the ratio of the desired probability of a non-synchronous spike to the probability of a non-synchronous spike under independence (the latter probability corresponding to the process we are thinning). After obtaining all N trials we then fitted the model to these simulated spike trains, found the estimateẑ, and applied the hypothesis test using the bootstrap method. This procedure was carried out for each z and N in our simulation. Because the simulation is computationally time-consuming, for the benefit of any future efforts along these lines, we also derived a formula to approximate the number of trials needed to get 0.8 power. Suppose we have N trials, each trial is T seconds, the bin size for synchrony detection is δ. Denote the instantaneous firing rates for two neurons on trial i by l ð1Þ t;i and l ð2Þ t;i . The number of synchronized spikes within the tth bin is y ð12Þ t;i and y ð12Þ t;i $ Poisson z Á l ð1Þ t;i l ð2Þ t;i Á d 2 , where z is the synchrony coefficient. The total number of observed synchronized spikes given l ð1Þ t;i and l ð2Þ t;i is N obs j l ð1Þ t;i ; l ð2Þ Assuming l ð1Þ t;i and l ð2Þ t;i are independent, we have E½l ð1Þ t;i E½l ð2Þ t;i d 2 ¼ NTl 1 l 2 d; where λ 1 and λ 2 are the mean firing rates of two neurons. Then we have We next assume that the distribution of logẑ is (approximately) normal, i.e., logẑ 2 N log z; 1 z 1 NTl 1 l 2 d , so that to get the power to equal 0.8 with type I error.05 we need where x is the threshold of rejecting null hypothesis log z = 0. We can then solve for N as the number of needed trials for detecting excess synchony: : Experiment Acute slice electrophysiology. Experiments were completed in compliance with the guidelines established by the Institutional Animal Care and Use Committee of Carnegie Mellon University. Whole-cell patch clamp recordings of hippocampal CA1 pyramidal cells were performed similar to previously described methods [52]. Briefly, a postnatal day 16 Thy1-YFP-G mouse [53] was anesthetized with isoflurane and decapitated into ice-cold oxygenated dissection solution containing (in mM): 125 NaCl, 25 glucose, 2.5 KCl, 25 NaHCO 3 , 1.25 NaH 2 PO 4 , 3 MgCl 2 and 1 CaCl 2 . Brains were rapidly isolated and sagittal slices (310 μm thick) containing the hippocampus were cut using a vibratome (5000 mz-2; Campden, Lafayette, IN, USA). Slices recovered for * 30 min in * 37C oxygenated Ringer solution that was identical to the dissection solution except for lower Mg 2+ concentrations (1 mM MgCl 2 ) and higher Ca 2+ concentrations (2 mM CaCl 2 ). Slices were then stored in room temperature oxygenated Ringer solution until recording. During recording, slices were continuously superfused with warmed oxygenated Ringer's solution (temperature measured in bath: 32°C). CA1 pyramidal cells were identified by morphology and laminar position using infrared differential interference contrast microscopy. Whole-cell recordings were made using electrodes (final electrode resistance: 5-7 MO) filled with (in mM): 120 potassium gluconate, 2 KCl, 10 Hepes, 10 sodium phosphocreatine, 4 Mg-ATP, 0.3 Na 3 GTP, 0.2 EGTA, 0.25 Alexa Fluor 594 (Life Technologies, Carlsbad, CA, USA) and 0.2% Neurobiotin (Vector Labs, Burlingame, CA, USA). The liquid junction potential was 12-14 mV and was not corrected for. Pipette capacitance was carefully neutralized and series resistance was compensated using the MultiClamp Bridge Balance operation. Data were low-pass filtered at 4 kHz and digitized at 10 kHz using a MultiClamp 700A amplifier (Molecular Devices, Sunnyvale, CA, USA) and an ITC-18 acquisition board (Instrutech, Mineola, NY, USA) controlled by custom software written in Igor Pro (WaveMetrics, Lake Oswego, OR, USA). Cell morphology was reconstructed under a 100X oil-immersion objective and analyzed with Neurolucida (MicroBrightField, Inc., Williston, VT, USA).
V4 neurons. Experimental procedures were approved by the Institutional Animal Care and Use Committee of the University of Pittsburgh. A separate analysis of these data has been previously reported ( [28,54]).
Subjects: We implanted one, 100-electrode "Utah" array (Blackrock Microsystems) in right V4 in one adult male rhesus macaque (Macaca mulatta). The basic surgical procedures have been described previously [55], and were conducted in aseptic conditions under isoflurane anesthesia. In addition to the microelectrode arrays, the animal was implanted with a titanium head post to immobilize the head during experiments. We recorded neurons with receptive fields centered * 4°from the fovea in the lower-left visual field.
Behavioral task: We trained the subject to maintain fixation on a 0.6°blue dot at the center of a flat-screen cathode ray tube monitor positioned 36 cm from its eyes. The background of the display was 50% gray. We measured the monitor luminance gamma functions using a photometer and linearized the relationship between input voltage and output luminance using lookup tables. The subject was trained to maintain fixation on the central dot for 2 seconds while no other visual stimulus was presented, at which time the fixation point was moved 11.6°in a random direction and the animal received a liquid reinforcement for making a saccade to the new location.
Microelectrode array recordings: Signals from the microelectrode arrays were band-pass filtered (0.3-7500 Hz), digitized at 30 kHz and amplified by a Grapevine system (Ripple). Signals crossing a threshold (periodically adjusted using a multiple of the root-mean-squared [RMS] noise for each channel) were stored for offline analysis. These waveform segments were sorted using an automated clustering algorithm [56] followed by manual refinement using custom MATLAB software [57] (available at http://www.smithlab.net/spikesort.html), taking into account the waveform shapes and interspike interval distributions. After sorting, we calculated the signal-to-noise (SNR) ratio of each candidate unit as the ratio of the average waveform amplitude to the standard deviation of the waveform noise [57]. Candidates with an SNR below 2.5 were discarded. Signals were also filtered from 0.3-250 Hz with a digital Butterworth filter and sampled at 1 kHz to provide LFPs.
LFP preprocessing: We assume that the oscillation modulating spiking activity is explicitly within the surrounding LFP. The naive way of selecting LFP is using the one recorded at the same electrode for each neuron. Since spike waveforms might contaminate the LFP spectrum [30,58], we computed LFP related to each neuron as the average of LFPs recorded on its neighboring electrodes. Another way of avoiding spike bleed-through is to choose the LFP on any electrodes adjacent to the neuron. In S2A and S2B Fig, we show that LFPs selected by all three methods are very similar. We also computed the spike-triggered average (STA) field potential using these three different methods. Their shapes are almost the same (S2C ad S2D Fig). We then bandpass filtered the LFP using Chebyshev type II filter design with passband 4-25 Hz. After we got the filtered oscillatory signal (Fig 8B-8E), we applied the Hilbert transform to estimate the instantaneous phase for further model fitting [20]. In this example, the firing rate is modulated by the magnitude of the oscillation B t = A t cos(F t ), where the amplitude A t is time varying and the modulation curve is 1 + B t = 1 + A t cos (F t ). We want to show that even though the phase-modulation assumption is violated, our method can still explain partly the role of oscillation in synchrony. (A) Amplitude and magnitude of the oscillatory signal; (B) Bootstrap-generated distribution of log z 12 values under the null hypothesis of log z 12 . Arrowhead shows the value of log z 12 predicted by the simplified model. A significantly larger number of synchronous spikes is observed than predicted by the model lacking an oscillatory factor. (C) Including an oscillatory factor in the model yields an accurate prediction of the observed number of synchronous spikes. (EPS) S1 Dataset. Experimental dataset used in this paper. Two datasets used for Fig 7 and Fig 8 were included in S1_Dataset. They were named as CA1_data and V4_data respectively. Details of data format were described in README file of each dataset. (ZIP)