|Home | About | Journals | Submit | Contact Us | Français|
The hypothesis that cortical networks employ the coordinated activity of groups of neurons, termed assemblies, to process information is debated. Results from multiple single-unit recordings are not conclusive because of the dramatic undersampling of the system. However, the local field potential (LFP) is a mesoscopic signal reflecting synchronized network activity. This raises the question whether the LFP can be employed to overcome the problem of undersampling. In a recent study in the motor cortex of the awake behaving monkey based on the locking of coincidences to the LFP we determined a lower bound for the fraction of spike coincidences originating from assembly activation. This quantity together with the locking of single spikes leads to a lower bound for the fraction of spikes originating from any assembly activity. Here we derive a statistical method to estimate the fraction of spike synchrony caused by assemblies—not its lower bound—from the spike data alone. A joint spike and LFP surrogate data model demonstrates consistency of results and the sensitivity of the method. Combining spike and LFP signals, we obtain an estimate of the fraction of spikes resulting from assemblies in the experimental data.
A common hypothesis concerning the processing of information by cortical networks involves the propagation of activity through synchronously firing groups of neurons, termed assemblies (Hebb 1949). A hallmark signature of an activated assembly is the repeated coincident firing of neurons in relation to a specific stimulus or behavioral demand resulting in coincidence counts that exceed the chance level estimated from the firing rates (Gerstein et al. 1989). Despite the inherent undersampling of state of the art multiple single-unit recordings, experimental studies indirectly substantiate the assembly idea with findings of such behavior-related significant synchronous spiking (e.g., Riehle et al. 1997; Kilavik et al. 2009). Independently thereof, a signal on the population level, like the mesoscopic local field potential (LFP), typically exhibits temporally structured oscillations commonly interpreted as correlated network activity. Synaptic transmembrane currents have been identified as the primary contributor to LFP generation (Mitzdorf 1985; Logothetis and Wandell 2004). Counterintuitively, although synchronized membrane potential oscillations of neurons in the vicinity (Katzner et al. 2009) of the recording electrode show strong correlations with the LFP (Poulet and Petersen 2008; Okun et al. 2010), these synchronized potentials do not induce the same degree of coincident spiking between the same neurons. Correlated spiking appears largely independent of synchrony on the level of membrane potentials (Tetzlaff et al. 2008; Poulet and Petersen 2008).
In a recent study (Denker et al. 2010) we were able to demonstrate in data of the motor cortex of the awake behaving monkey the missing link between significant spike synchrony and the LFP. A conceptual model enabled us to derive lower bounds for the fraction of spike coincidences β originating from observed assembly activity identified by periods containing excess spike synchrony, and the fraction of spikes of neurons γ caused by assembly activity whether observed or not. The results were obtained by comparing the locking of spike coincidences to the LFP in time periods with significant synchrony to the locking outside of these periods.
In the remainder of this section we first introduce the Unitary Events analysis method (Section 1.1) used to quantify significant excess spike synchrony. We then review our earlier findings (Denker et al. 2010) on the phase locking between spike-spike coincidences and the LFP (Section 1.2) and the model-based interpretation of the data (Section 1.3). In the present work we remove the limitation to lower bounds of β and γ by evaluating a stochastic model of the composition of spike coincidence counts. Section 2 introduces the model and derives an estimator for β only based on the configuration of spikes in the respective time interval. In Section 3 we explain that the probability of spikes γ to be part of an assembly—in contrast to β—cannot be extracted from the spike trains alone. However, equipped with the estimate of β we are in the position to compute γ using the previously found relationships between the phase distributions of the spike-LFP coupling (Section 1.3). The possibility to extract a parameter that relates to macroscopic features of the network dynamics illustrates the relevance of the LFP as a population signal in addressing the undersampling problem.
We demonstrate the consistency of the concept with the help of a joint spike-LFP toy model. In the physiological range, the parameters can reliably be determined and the lower bounds obtained from the phase locking are compatible with the parameter estimates yielded if the spike statistics is considered in addition. Finally, we utilize our new tool to reanalyze the experimental data (Section 4) for β and γ. Although highly simplified the toy model enables us to discuss (Section 5) various aspects of the distribution of parameters obtained from the experimental data set.
In the following we briefly summarize the main results obtained in Denker et al. (2010) which serve as a starting point for our analysis. In the first step we analyzed simultaneously recorded single units from monkey motor cortex (cf., e.g., Roux et al. 2006 and Section 4.1 for details) for excess spike synchrony by employing the Unitary Events analysis method (Grün et al. 2002b; Grün 2009). The method compares the empirically measured spike coincidences to the number expected by chance given the neurons’ firing rates. The expected number defines the mean of the Poisson distribution realizing the null-hypothesis of statistical independence. If the p-value of the empirically found number of coincidences evaluated by comparison to this distribution exceeds the significance level α, the coincidences are considered significant and are termed Unitary Events (UE). This analysis is performed in a sliding window fashion (windows of 100 ms) to account for the dynamics of correlations and the non-stationarities in the firing rates in time (cf. Riehle et al. 1997 for details). Per sliding window data are considered from the various trials, and typically only a small fraction of the spikes in a window is coincident (Grün et al. 2002b, 2010). Non-stationarity across trials is accounted for by calculating the expected number of coincidences as the sum of trial-by-trial expectancies (Grün et al. 2003). Coincident spike events with a temporal precision of up to 3 ms were collected by the multiple-shift method to avoid artifacts due to binning (Grün et al. 1999). As a result of this analysis we are able to detect the dynamics of spike synchrony and identify time intervals that exhibit UE, i.e. contain excess spike synchrony (compare Fig. 1(a), top). Non-UE time periods contain coincidence spike events at chance level only. Thus UEs indicate time windows where the neurons under consideration exhibit coordinated activity across trials. The (unknown) set of all neurons within the network that generate this correlated activity is termed the assembly.
In the next step we investigated the relation of spikes to the LFP. In particular, we were interested if spikes in synchronous events, i.e., chance coincidences (CC) or UE, have a different relation to the LFP than spikes that are not involved in a coincidence, i.e. isolated spikes (ISO). Therefore we classified each spike of a data set into one of these three classes and analyzed the classes separately (see Fig. 1(a) for an illustration; color codes for ISO, CC, and UE are retained throughout this article). Spike triggered averages (STA) of the LFP revealed for all classes an oscillatory structure at about 17 Hz with spikes occurring preferentially at the decaying amplitude. However, the amplitude of the STAs derived for the different spike classes were strikingly different: consistent across the whole data set the STA triggered on UE spikes exhibited the largest amplitude, for CC spikes a smaller amplitude and for ISO the smallest.
However, the STA analysis cannot uncover the reason for the differences in STA amplitude, since it may be due to differences in the LFP amplitudes or due to different degrees of phase locking between spikes and the LFP. In order to disentangle these aspects we performed a phase-amplitude analysis of the LFP. Using a Hilbert transform of the LFP we gained the instantaneous phase and amplitude as time dependent functions, and extracted the respective measures at spike times (the detailed procedures will be explained in Section 4.1). The phase distribution exhibits the phase preferences by non-uniformity of the phase histograms (see Fig. 1(b)), and we observe that UE spikes express the strongest degree of phase locking, CC less (at the predicted chance level) and even less for ISO. A resampling procedure ensures that histograms are compared with comparative numbers of spikes: The phase histogram of CC (UE) is compared to resampled versions of the ISO (CC) phase histogram, each containing the same number of data points as the CC (UE) histogram (bands in Fig. 1(b), middle and right graph). The differences in phase locking were found to be consistent across the whole population of recorded neurons, and 68% of the neurons (n=123 neuron pairs; see Section 4.1 for details) exhibit a stronger phase locking of UE than CC. In contrast we found much smaller or even negligible differences of the LFP amplitude measured as its envelope at spike times for the different spike classes.
Our analysis of the relation of spikes to the LFP revealed several unexpected results. Firstly, we found a difference in the locking degree for spikes involved in different categories of coincidences, i.e. UE vs CC. Following the hypothesis that active assemblies are expressed by coordinated spiking activity, UE coincidences are interpreted as a signature of such active assemblies. In the following, we will define the assembly as the set neurons exhibiting coordinated activity that is reflected in the excess synchronous events during UE periods. Our new finding on the phase locking of UE spikes to the LFP leads to the interpretation that assembly activity occurs in a pulsed fashion, locked to the LFP oscillation. The fact that spikes of chance coincidences also occur phase locked seems contradictory. However, their degree of locking is fully explained by the (weak) locking of single spikes, which trivially leads to an enhanced locking of such spikes if occurring coincidently (cf., predictor in Fig. 1(b)). Thus the question is rather, why isolated spikes outside UE periods exhibit phase locking at all.
To get a better understanding of these puzzling observations, we developed a conceptual model which consistently explains all our findings (Fig. 2(a)). A direct conclusion of the experimental data, which show that UE periods exhibit a better locking to the LFP than CC periods, is that assembly spikes must be better locked to the LFP than spikes not participating in an assembly. Hence, one basic assumption of the model is that spikes involved in assemblies exhibit locking to the LFP while non-assembly spikes do not. According to our findings that individual spike trains typically contain a mixture of all spike categories (ISO, CC, and UE), it is reasonable to assume that an individual spike train is composed of assembly and non-assembly spikes. Due to the severe undersampling of the system caused by the limited number of recording electrodes, it is highly likely that also ISO spikes contain assembly spikes, but the corresponding partner neurons could not be identified. Consequently, the phase histogram of ISO must contain a mixture of unlocked and locked spikes. Assuming a uniform phase distribution pn(ϕ) for non-assembly spikes and an unknown non-uniform phase distribution pa(ϕ) for assembly spikes, the phase distribution of the ISO spikes pISO(ϕ) can be expressed by a weighted combination of the two components:
Thus also ISO spikes show locking, and its strength is determined by the factor γ.
CCs result from pairs of independent spike trains that express coincidences by chance. Under the above assumption that spike trains are composed of assembly and non-assembly spikes, chance coincidences are the result of combinations of different spike types: unlocked–unlocked spikes, locked–unlocked spikes (two possible combinations), and locked–locked spikes. Thus, pCC(ϕ) is also expressed as a composition of pn(ϕ) and pa(ϕ):
In case of time periods that contain UEs, coincidences resulting from an active assembly are present, but likely intermixed with chance coincidences that are not related to this assembly. Thus, the phase distribution of UE is assumed to be composed of a contribution of chance locking of a percentage of (1−β), and to a percentage of β of assembly coincidences (see sketch in Fig. 2(a)):
As a result, the enhanced locking of UE as compared to CC is a consequence of the presence of assembly coincidences i.e. β>0.
In terms of spike coincidences, β is the ratio of the number of coincidences resulting from assembly activation and the total number of coincidences in a given time window
The task therefore is to construct an estimate of nc on the basis of the known properties of the observed spike trains. For two independently spiking neurons with n1 spikes in one spike train and n2 spikes in the other, the probability to observe k coincidences in a time window T is given by the hypergeometric distribution (Grün et al. 2003)
where we introduced the shorthand
and h is the resolution of the discretized time axis. We define outside of . The expected number of coincidences is
The expectation value (Eq. (7)) is just the probability to observe a coincidence in a given bin multiplied by the number of available time bins Th. In the presence of nc deterministic coincidences, however, only a reduced number of spikes of the two neurons is available to form chance coincidences
The construction of the respective spike trains is illustrated in Fig. 2(b). Let us now assume that we measure the number of coincidences averaged over many repetitions of the experiment with exactly the same values for nc, n1, n2. In this case we can equate the empirical average with the expectation value . In this expression nc is the only unknown variable. Hence, we can express nc in terms of n1, n2, and . Given only the triplet of a single realization we can still hope that the measured nemp is typical and write
shows that the quadratic terms in nc cancel
and collecting terms with nc gives
Using Eq. (9) we can already compute β for all experimental UE periods with their individual triplets and obtain the distribution of β as well as its mean βUE (see Section 4). Clearly Eq. (9) is an approximation, the expression is negative for small nemp, in particular when no coincidence has been observed (nemp=0).
In order to capture coincidences with a temporal jitter larger than the resolution of the data h, the multiple-shift method (Grün et al. 1999) sums the coincidences over a range of s=2b+1 displacements −b,...,+b of one spike train with respect to the other. At each shift only some of the nc coincidences are detected and a particular coincidence is by definition only detected at a single displacement. Under the assumption that in the s displacements the spike counts n1and n2 are conserved, the relation (Eq. (8)) between nemp and nc holds if we replace the number of spikes available for chance coincidences by sni−nc and the number of available bins by sTh−nc. The estimate of the number of coincidences originating from assembly activation (Eq. (9)) then reads .
Next we consider a slightly more realistic model with a variable number of assembly activations nc. The purpose of the model is to help us to understand the conditions under which we can reliably extract β. In the absence of any knowledge about the process generating the additional coincidences nc we assume that each value consistent with the observed number of coincidences is equally likely. Thus, each possible value of nc occurs with probability 1/(nemp+1). For a particular nc, however, the probability to be consistent with the nemp observed coincidences is again given by the hypergeometric distribution (Eq. (5))
The probability that a particular nc underlies the observation therefore is
and the total probability to observe nemp coincidences is
Consequently, the expected nc given the observed triplet is
which reduces to
Note that the denominator is not unity as the sum extends over probabilities from different hypergeometric distributions.
The intuitive approximation (Eq. (9)) is compatible with the more detailed model result (Eq. (10)). To see this, we first approximate the individual hypergeometric distributions in Eq. (10) by the one for the average nc and extend the range of the sums to the maximum value Th−nc:
This reduces the denominator to unity. With the substitution k=nemp−i this reads
The sum only has contributions for k≥0 and we again extend the range to the maximum value
This suggests a decomposition of the sum to
where the first sum is unity and the latter term is the expectation value of k. Hence,
Before we turn to the experimental data in Section 4 we need to assess the accuracy of our estimator of nc. To this end we construct a surrogate data set with parameters adapted to the experimental data. The data are organized into 32 blocks containing an increasing number of assembly activations from nc=0 to 31. A block is composed of M=2,700 time windows, consistent with the number of UE windows found in the experimental data. A time window has a duration (Eq. (6)) of Th=5,000 at a resolution of h=1 ms, corresponding to the 100 ms segments of UE analysis covering 50 trials. Each time window contains surrogate spike trains of two neurons with totals of exactly n1=n2=100 spikes (corresponding to a rate of 20 Hz). In each spike train the ni−nc non-assembly spikes are uniformly distributed over the time window.
Figure 3(a) shows the distribution of injected coincidences nc in the data set organized by the total number of coincidences in each window nemp. Here and in the following only windows containing a significant number of coincidences (nemp>nα) are analyzed (cf. Section 1.1). The estimator (Eq. (9)) of nc assigns a unique value to each value of nemp because the spike counts are identical for all windows. Therefore at a given nemp the estimate of nc is identical for each window and thereby identical to the average over all windows. In Fig. 3(b) we observe that the estimate well describes the mean (red curve) of the actual distribution of nc at a particular nemp. The panel further demonstrates that the approximative (Eq. (9)) and the exact (Eq. (10)) estimator only start to deviate for nemp below nα. Figure 3(c) uses the same representation as Fig. 3(a) to illustrate the excellent correspondence between the distribution of the actual β and the estimated values. The total distribution of β (Fig. 3(d)) in the model is asymmetrical because of the lack of a typical nc and the constant spike count. In conclusion, our measure is a faithful estimator of the average β in the data.
In the following section, we extend the simple spike model defined in Section 2 to include a representation of the LFP locking and verify numerically that the results obtained in the previous section yield a reliable estimate of our model parameters β and γ.
First, we choose a parameter γset as the assumed probability that any spike in the network is part of an assembly activation, in agreement with our conceptual model. In the following we simulate for each choice of nc the set of M=2,700 UE time windows as combinations of an injection process and a background process for two neurons as described in Section 2.2. In addition, CC time windows are generated only by a background process with equal spike counts n1=n2=100 (see Fig. 2(b) for an illustration). Again, only significant UE windows (nemp>nα) and non-significant CC windows (nemp≤nα) are retained. To model the experimental results, we assign a label ’a’ to all spikes that originate from an assembly activation (Fig. 2(b)). By definition of our assembly process, every spike that originates from the injection process receives a label. In addition, a random proportion γset of spikes from the background process is labeled (in both, UE and CC time windows). In our simulations, the overall probability for a spike to belong to an assembly is set to γset=0.1. Next, we define two distributions pn(ϕ) and pa(ϕ), where the latter has a larger modulation depth, which describe the locking of non-assembly and assembly spikes, respectively (Fig. 2(a)). Here, pn(ϕ) is modeled as a uniform distribution, whereas for pa(ϕ) the distribution is modeled as a Gaussian as an approximation of a von Mises distribution. The modulation of the Gaussian was chosen to mimic that of the experimentally observed distributions pCC(ϕ), and pUE(ϕ) (compare Figs. 1(b) and 2(a)). Classification of spikes into the groups ISO (taken here as the single spikes in CC windows), CC and UE allows us to calculate the simulated phase distributions pISO(ϕ), pCC(ϕ), and pUE(ϕ) as mixtures of pn(ϕ) and pa(ϕ).
The setup allows us to follow the same analysis steps that we will perform on the experimental data in the following section. From the conceptual model introduced in Section 1.3, which is formally expressed in Eqs. (1–3), we infer the lower bound of coincidences originating from assemblies during an observed UE period. Substituting the measured population phase distribution of chance coincidences pCC(ϕ) of Eq. (2) into Eq. (3) yields an expression relating the known phase distribution pUE(ϕ) of UE coincidences to the parameter β and the squared phase distribution of assembly spikes :
By systematic variation of the parameter β we thus obtain a corresponding phase distribution by solving the equation separately for each bin of the respective distributions. However, for small values of β the assembly distribution must exhibit a strong modulation to compensate for the large difference between pUE(ϕ) and (1−β)pCC(ϕ). In the extreme case, it can become necessary for a bin of to contain negative values, ruling out that particular choice of the parameter β. From this consideration, we define as the lowest value of β that leads to a phase distribution with non-negative entries. This lower bound depends on the choice of the assembly phase distribution pa(ϕ) we initially introduced into our model: the smaller the difference in locking between assembly spikes pa(ϕ) and non-assembly spikes pn(ϕ), the lower are the values that may assume. In Fig. 3(d) we show that in our spike-LFP model we obtain a value for well below the mean βUE extracted from the spike analysis (Section 2). We note that when generating the model data with a more strongly modulated Gaussian to model pa(ϕ), the minimum approaches βUE from below (not shown).
We are now prepared to extract the parameter γ from the simulated phase distributions using the estimate of βUE derived in Section 2. First, using β=βUE we extract from Eq. (3) by solving the equation separately for each bin of . Thus, pa(ϕ) is determined independently of the distribution pn(ϕ) describing the general entrainment of spikes to the LFP. Taking the square root pa(ϕ), we renormalize the distribution to unit area, and insert it in Eq.(1). By variation of the parameter γ, we find the value γϕ,UE that minimizes the sum of the absolute bin-by-bin differences between the distribution pISO(ϕ) and the right side (1−γ)·pn(ϕ)+γ·pa(ϕ) of Eq. (1).
Using this method, we first analyze separately each of the 32 blocks of data with a fixed number of injected coincidences nc. Let and be the estimated values of βUE and γϕ,UE within a single block. In Fig. 3(e) and are compared to their theoretical means (calculated as the mean of nc/nemp over all significant time windows from one set of simulations with fixed nc) and γset. The estimates are in good agreement for datasets where nc was set above the significance threshold nα. Below this threshold, the number of injected coincidences, and hence and , is overestimated. Clearly, in this regime contributing UE periods become significant only due to an unusually high number of coincidences in the background process. Therefore, the mean approximation (Eq. (9)) does no longer hold. In principle it is possible to correct for the bias due to the selection of significant periods, although it is not possible to arrive at an expression in closed form for nc. Nevertheless, in fact only a small fraction of windows that were created with low values of nc actually become significant due to such exceptionally high coincidences counts from the background (the shaded graph in Fig. 3(e) shows the probability for a window to become significant for given nc). Therefore, in practice, where the true nc is unknown, only few windows enter the analysis when nc is low. For this reason, the bias does not affect the estimate when re-sorting the data according to the observed nemp (Fig. 3(a)–(c)), and the original value of γset is well estimated as from the complete dataset spanning all nc (black cross in Fig. 3(e)). Consistently, the lower bounds and are a lower bound on the respective estimates and for all values of nc.
To quantify whether our procedure works for other values for the total spike counts n1 and n2 we introduce two error measures. First, we note that our prediction of β will be least accurate when the coincidence count nemp is low. Thus, we numerically quantify the absolute error of β for the smallest value of nemp to become significant (Fig. 3(f)). Moreover, we quantify the absolute error of our final estimate γ obtained from the complete spike-LFP model (Fig. 3(g)). Both measures indicate that the parameters are reliably extracted for a wide range of activity levels in both neurons, even in the case of large rate differences between the two. In summary, these calibrations using a simple combined spike-LFP model demonstrate that our method and the approximation (Eq. (9)) are well suited to estimate the parameter γ from the phase distributions in a dataset with a realistic choice of parameters.
In this section we estimate the percentage of coincident spike events reflecting assembly activity (model parameter β) and the percentage of spikes that are part of an assembly activation (model parameter γ) from neuronal data of primary motor cortex of monkeys.
Two rhesus monkeys were trained to perform arm movements in two different tasks involving an instructed delay. One of these tasks requires a self-initiated movement after one of two preset delays, while the other requires the discrimination between two such delays. In this study we exclusively analyze the delay activity during the preparatory period for the upcoming arm movement. Only correct trials were considered, in which the monkey responded within a predefined time window and in which movements were performed in the required movement direction. LFPs and spikes were recorded simultaneously in primary motor cortex using a multielectrode device of 2–4 electrodes. Spikes of single neurons were detected by an online sorting algorithm. The inter-electrode distance was on the order of 400 μm. LFPs were sampled at a resolution of 250–500 Hz and hardware filtered (band pass, 1–100 Hz). In total, we analyzed 53 recording sessions, which yielded 143 single neurons or 570 combinations of neurons and behavioral conditions. Each of these eight possible behavioral conditions is a combination of the task during which the neuron was recorded, the length of the delay (short or long), and the direction of the arm movement (left or right). This selection included only those neurons which exhibited an average firing rate of 5 Hz or more and exhibited at least 25 spikes in total. On average 33±11 (mean±standard deviation) trials were recorded per experimental condition. In analyses that combine spikes and LFP, each neuron enters only once, and we never combined LFP and spikes that were recorded on the same electrode to exclude the possibility of spike artifacts in the signal. We confirmed that simultaneously recorded LFPs are highly synchronous in the frequency regimes of interest. For experimental details see Roux et al. (2006).
The Unitary Events analysis was applied to all simultaneously recorded pairs of neurons recorded on different electrodes thus totaling 123 pairs of neurons. Data were analyzed separately for the possible behavioral conditions and each experimental session. Defining a neuron by the combination of its identity and the behavioral condition (see above) during which it was recorded, data from the same neuron may enter a population average up to eight times. Based on the results of the UE analyses each individual spike was marked as either ISO, CC or UE (see Section 1.1 and Denker et al. (2010) for details). We observed a dominant component of the β-band (in both monkeys around 17 Hz) in the LFP during the preparatory period (see, e.g., Murthy and Fetz 1996a) and therefore filtered the LFP accordingly before applying the Hilbert transform to extract its instantaneous phase and amplitude (i.e. the envelope). LFPs of both monkeys are filtered with a zero-phase 10–22 Hz band pass filter (Butterworth, 8-pole). Finally, we derived the phase histogram pooled over the whole population for each class of spikes (see Fig. 1(b)). The histogram was pooled after averaging, such that each spike enters the distribution unweighted (as is required for our analysis).
As described in Section 3.2, we derive the lower bound of the percentage of coincidences originating from assembly activity . Again, we make use of the Eq. (11) that relates the phase histogram of the UE to the model parameter β. This equation contains two components that we extract from the data, i.e. the phase distribution of the UE spikes pUE(ϕ) and the phase distribution of the CC spikes pCC(ϕ). We vary β in the interval [0, 1] and retrieve the corresponding . Four examples are shown in Fig. 4 for different values of β. The first two examples demonstrate how choosing a small value of β may lead to distributions that have negative entries and indicate a non-valid solution. From the systematic variation of β we derive the value that just leads to non-negative entries in all bins of (see Fig. 4, bottom). This result indicates that at least about a fourth of the coincidences during a UE period are reflections of the observed assembly. In addition, we performed an analogous analysis using von Mises distributions fitted to the experimental distributions pCC(ϕ) and pUE(ϕ) (dashed lines in Fig. 4). This approach is very conservative because the fits neglect the prominent systematic peaks present in the original distributions. In this analysis, the minimum value for β is determined as .
As introduced in Section 2 the parameter β can be estimated by the comparison of the expected number of coincidences nexp and the number of empirical coincidences nemp. Figure 5(a) shows nemp as a function of its respective nexp value for all analysis windows that were detected to contain UE for all pairs and sessions. Due to the selection of windows that contain UE, the empirical coincidence counts are larger than the minimal number of coincidences nα required to be significant given the significance level α. nα does not have a constant difference to nexp but increases non-linearly as a function of nexp since the Poisson distribution used for the evaluation of the significance becomes broader for larger expected number of coincidences (see for details in Grün et al. 2002a).
Based on Eq. (9) we estimate the number of coincidences resulting from the active assembly nc for each nexp, nemp combination. We use a version of the UE analysis that also allows to detect temporally imprecise coincidences by employing the multiple-shift method. This method detects coincidences of systematically shifted spike trains up to a predefined shift, and sums the coincidence counts from all shifts. For the UE evaluation, the expected number of coincidences then has to be adjusted by a factor of 2s+1 with s being the maximal shift (cf. Section 1.1; for details see Grün et al. 1999). Thus to estimate the number of coincidences nc that generically result from active assemblies, Eq. (9) was adjusted correspondingly (Section 2). Figure 5(b) shows nc versus its corresponding nemp and we find an average number of assembly coincidences of .
Next, we can derive for each pair of (nc, nemp) an estimate of the percentage of coincidences reflecting assembly activity β (Fig. 5(c)). We find these values peaked around a mean of βUE=0.51 with a standard deviation of 0.123. The obtained distribution of β is in good agreement with the minimum obtained in Section 4.2, as lies well below the mean βUE of the distribution of β as expected for an estimate of a lower bound. Although in fact seems to be a bound for the complete distribution of β, this must not necessarily be the case: As a measure derived from the population estimate of phase distributions, it can only have predictive power on the population mean βUE, but not the estimate β of single windows (compare Fig. 3(d)).
Finally, we estimate the fraction γ of spikes in the network that are part of an assembly activation (compare Section 3.3). Again, by systematic variation of the parameter γ in Eq. (1) we find the best fit between the measured pISO(ϕ) and the right side of the equation (see Fig. 6). Here, the assembly distribution pa(ϕ) is known from Fig. 5 by setting β=βUE, while the phase distribution of non-assembly spikes pn(ϕ) is taken as the uniform distribution (see Section 5). Taking all analysis steps together, the best fit is derived for γϕ,UE=0.22 indicating that on average 22% of the spikes of any neuron in the network originate from assembly activity (Fig. 6, bottom). To corroborate this result, we obtain a similar value of γϕ,UE=0.26 using the same procedure described above to solve (Eq. (2)).
Despite the complex mechanisms that contribute to the formation of the local field potential, it is well established that a primary contribution to the oscillatory LFP dynamics arises from the superposition of synchronized, slow transmembrane currents of cells close to the recording site (Mitzdorf 1985; Logothetis and Wandell 2004). Nevertheless, how rhythmicity in the LFP is linked to synchrony on the spiking level has remained an open question (Poulet and Petersen 2008; Tetzlaff et al. 2008) due to the unspecificity of the LFP signal and lack of a clear global oscillatory spiking activity. Several authors have interpreted the LFP as reflections of the specific synchronous synaptic activity responsible for the coactivation of neurons in the context of cell assemblies (Eckhorn et al. 1988; Murthy and Fetz 1996b) which, in the simplest case, time their activations to the network rhythm revealed by the LFP (Singer 1999). Our recent experimental findings (Denker et al. 2010) demonstrate that indeed only assembly activity, which is identified as transient periods of significant excess spike synchrony between two neurons (UE analysis, Grün et al. 2002a, b), shows an exceptional phase relationship to LFP oscillation that exceeds expectation (Fig. 1(b)), thus confirming the hypothesis.
This link provides a handle on characterizing the spike synchronization dynamics in the context of the network oscillations. Here, we introduced a simple model that captures the main experimental findings on the spike-LFP relationship in such a way that it includes a parameter γ that measures the overall participation probability of individual spikes in an assembly—independent of whether it is observed or not. However, the estimation of this network parameter from the data requires knowledge of a second model parameter β, the expected relative number of (excess) coincidences stemming from an active assembly (as opposed to chance coincidences) during UE periods. Estimating this number from the spike data alone and integrating it into the conceptual model yields an average participation of single spikes to assembly activity of γ=22%, a parameter that describes precise spike synchronization on the network level.
Our model assumes that the observed excess synchrony is the result of the specific activation of the observed neurons. An alternate hypothesis states that spike synchronization is solely caused by fluctuations of the spiking likelihoods due to the entrainment of neurons to a common LFP oscillation. In this case, no UEs would be observed due to the lack of a spike coincidence count exceeding the expectation based on the rate time course. Only if the analysis would unsuccessfully correct for non-stationarities of the firing rates false-positive detections of UE periods emerged. Consequently, these false-positives would trivially correlate with the LFP. To exclude this possibility, we repeated the analysis presented in this study by replacing the parametric distribution of coincidences used for evaluating the significance of observed coincidence counts by numerically derived distributions based on surrogate data. The employed surrogate method (spike train dithering, see Grün 2009) considers further statistical features of the experimental spike trains (in particular non-stationarity of rates on a short time scale and the inter-spike interval distributions) while destroying precise spike coincidences. Simulations showed that this method is more conservative (Louis et al. 2010) than the analysis used here, but despite the decreased sensitivity confirms the differences in the phase distributions for ISO, CC, and UE. These findings are the essential experimental foundation of our study.
The combination of measurements of synchrony on the local and mesoscopic scales enables us to access parameters of the network dynamics that remain hidden on the two individual levels of observation. The nature of the estimation process via the population phase distributions obtained from the LFP-spike locking requires a vast pooling of the experimental data in different sessions and two different monkeys in order to obtain a large sample of neurons as an appropriate representation of the network activity in motor cortex. The resulting value for γ must therefore be seen as a coarse estimate of the degree of assembly activity. Despite the finding that the parameter β shows a rather narrow distribution (Fig. 5(c)), the extent to which an individual neuron contributes to the assembly dynamics will fluctuate around the network mean γ. By quantifying the quantiles for 2 standard deviations of the experimental distribution of β we find corresponding γ values between γ=0.14 and γ=0.39 as an estimate of the range of γ across neurons.
In understanding the underlying dynamical structure, more interesting than the precise value of γ itself are two simple observations: first, γ<1. Therefore, not all spikes are part of an assembly activation, and some spikes must be attributed to a complementary mechanism. Intuitively, this is clear from the observation that in individual UE periods, it is always the falling phase of the LFP oscillation where increased locking is observed. Second, γ>0, specifically about one fourth of all spikes must originate from assembly activity (with a non-zero lower bound ). Thus, under the assumptions of our model, the activation of assemblies is an ubiquitous phenomenon in the network, providing compelling evidence for the presence of an assembly coding scheme in addition to the correlation of synchrony with behavior (Kilavik et al. 2009). Indeed, we typically observe a fraction of about 26% of neurons that show UEs during a given task. Therefore, combined with the large value of γ this suggests that typically any neuron is part of one or several assemblies. Nevertheless, information about the existence of higher-order synchrony between neurons is required for a further characterization of the assembly, such as an estimate of the number of participating neurons. Moreover, such an analysis could better disambiguate whether spikes that do not originate from assembly activation show an intrinsic degree of phase locking that deviates from our model assumption of a uniform pn(ϕ).
A crucial part of our analysis is the estimation of the average relative amount of excess synchrony β during UE periods directly from the synchrony analysis. By design, all currently available methods to detect the presence of an active assembly activation, such as the UE analysis, rely on a significance test (Grün 2009). Therefore, we must assume that the amount of excess synchrony β is influenced by the explicit choice of the significance level α used for the significance test: For a very restrictive significance criterion, only UE periods with very high excess synchrony nc will be detected, resulting in higher values of β. Nevertheless, by choosing the same α level for the phase locking analysis, we will in turn also obtain a correspondingly more modulated distribution for pUE(ϕ), reflecting that it contains a higher proportion of assembly spikes. Therefore, β yields a consistent phase distribution pa(ϕ) of the assembly spikes independent of α.
The approximate formula (Eq. (9)) does not make any assumptions on preselecting significant periods of synchrony in the first place. However, we apply this estimate specifically in periods that display a significant surplus of synchrony, i.e. UE periods. If the number of coincidences nc coming from the injection process is small, significance is only reached with an unexceptionally high amount of coincidences originating from the background. Therefore, in our stochastic spike-LFP model, in those very few trials that are significant for a given small nc, the average background rate will be underestimated and the average injection rate overestimated (Fig. 3(c)). However, restructuring the model data to the realistic case where only nemp is measured (pooling across the data sets with different nc) automatically incorporates the low frequency of significant windows with small nc. Therefore the overall estimate of β is not significantly affected.
Our spike-LFP model (Section 3) is created in the spirit of providing a simple abstraction of the experimental findings in order to test the analysis under controlled conditions (equal counts per spike train with fixed injections). Due to this intentional lack of experimental detail in the model, the distribution of β that is obtained in the spike-LFP model naturally deviates from the one found in the data (compare Figs. 3(d) and 5). The model assumes an equal probability for a large range of nc (number of injected coincidences). In real data, nc likely follows a much more narrow distribution that does not exploit this range, resulting in lower values of β (compare Fig. 3(e)). Moreover, we assumed fixed counts n1, n2 for all neurons in the model. In reality, rates vary considerably across neurons. Especially for higher rates, where a higher proportion of coincidences can be assumed to originate from the background, we would expect a tendency towards lower β values.
In generating the spike-LFP model in Section 3 we did not explicitly model an LFP oscillation to place spikes at specific points of the field potential. In extreme situations this might be an oversimplification, where the constraints placed on the spikes due to the LFP locking influence the probability to detect coincidences. In our model, however, non-assembly spikes are associated with a uniform phase distribution, such that pn(ϕ) does not impose any constraints on the spiking probability in time. The Poissonian spike interval statistics of non-assembly spikes thus remain unaffected. In contrast, the assembly spikes from the background process must in principle be adjusted with respect to a hypothetical LFP so that their phase distribution matches the assumed the distribution pa(ϕ) of assembly spikes. However, as pn(ϕ) is uniform, doing so will only influence the probability of finding a coincidence between two spikes that are both assembly spikes, i.e. that participate by chance in two different assemblies at the same time. Nevertheless, the probability (γset)2=0.01 of this to happen is negligibly small. Finally, assembly spikes from the injection process are synchronous by definition, and therefore remain unaffected by the choice of pa(ϕ). Due to their small number they can always be freely placed in the time window T in accordance with pa(ϕ). Taken together, the simple model introduced in Fig. 2 to represent the experimental findings in the context of the spike model includes sufficient detail without the need to explicitly model the actual positions of individual spikes with respect to an artificial LFP.
The analysis we provide in this manuscript focuses on estimating a global parameter that characterizes the network dynamics. Given the overall consistency of our conceptual assembly model, the results suggest that oscillations and assemblies are tightly coupled. However, one may still speculate whether LFPs are in fact the direct cause of the assembly process, or whether they act as a supporting mechanism that coordinates synchronized firing. Whichever the initial cause of the LFP oscillation, a promising next step is to return to the level of the single neuron and integrate our knowledge on the dependence of spike synchrony on the LFP in order to discriminate which of the spikes are likely candidates to originate from the assembly process, and which are not. Two pieces of information aid us in this task: first, the knowledge of the observed UE phase distributions in relation to the ISO and CC distributions determines the likelihood of a single spike to be part of an assembly activation. Second, the parameter γ is informative of the relative number of spikes to assign to the assembly dynamics. Indeed, a third parameter not discussed in this manuscript is the dependence of spike synchrony on the LFP amplitude (Denker et al. 2010). Evaluating these three factors on a spike-by-spike basis, there is reason to believe that analyzing the LFP in relation to the spiking activity may well provide an independent indicator to identify the probabilities with which recorded neurons become transiently synchronized as part of an assembly activation. Moreover, the pulsed activation of assemblies triggered on the LFP oscillations cycle suggested by the experimental data stimulates the idea that maybe even assemblies, where only one neuron is recorded, may be inferred from the phase locking statistics in a probabilistic manner. Thus, incorporating a mesoscopic signals may serve to overcome problems related to the undersampling in multiple single-neuron recordings. The success of the method presented here not only corroborates the assembly hypothesis of neuronal processing, but offers a promising vista on reinterpreting the dynamical implications of observed LFP signals.
Partially funded by the Helmholtz Alliance on Systems Biology, the French National Research Agency (ANR-05-NEUR-045-01), EU grant 15879 (FACETS), and the Next-Generation Supercomputer Project of MEXT, Japan. Part of the work was carried out while SG and MDi enjoyed a scientific stay at the Norwegian University of Life Sciences, Ås in January, 2009.
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.