Synchronous spike events are widespread in multi-electrode systems

To demonstrate how artifacts appear in massively parallel spike train data, we start inspecting two macaque data sets from different laboratories and animals (Fig. 1a). In the raster display of the simultaneously recorded spike trains (Fig. 1b, top), we notice events where many neurons exhibit spikes simultaneously (highlighted in gray). To examine the temporal precision of these spike coincidences, we compute the population spike time histograms with different bin sizes (Fig. 1b, bottom). At the resolutions of 1 and  ms, we find some bins with outstandingly high spike counts, not visible at the resolution of 10 ms. We term these extremely precise spike coincidences within  ms as HSEs.

We further quantify the HSEs by their complexity (i.e., number of involved neurons; Grün et al., 2008) and compare the obtained complexities to those expected by chance using a surrogate method (Stella et al., 2022; Fig. 1c; see “Materials and Methods” for details). Figure 1d shows the distributions of complexity values from the original and the surrogate data, for two different bin sizes (1 and  ms). For both bin sizes, the original data contains more high complexity HSEs than expected by random chance. Considering that the coordination of spike times on the timescale of  ms is very unlikely, we conclude that many of the HSEs in the experimental data are artifacts, which we name synchrofacts.

We find artifacts in two separate data sets originating from two different laboratories, recorded in two different cortical areas of different macaques. Other arrays in these recordings, but not presented here, also show similar excess HSEs (see Supplementary Material: Fig. S7). Therefore we suspect that artifacts are common across recordings of the same type. While both data sets were recorded with Blackrock recording systems, artifacts also occur in other recording systems, e.g., in Neuropixels recordings (see section “High-density probe recordings also contain artifacts”). Taken together, our observations highlight the ubiquity of artifacts.

One might think that the data can be cleaned by simply removing all HSEs. However, that is not an appropriate solution to get rid of artifacts, since this also removes many genuine neuronal spikes and, thus, genuine neuronal synchrony. One can demonstrate this by removing all spikes in the bins with complexity of two or larger at resolution, and then calculating the complexity distribution, to be compared to the complexity distribution of its surrogates. The surrogate distribution obtained in this way contains complexities of two or higher that occur by chance. However, the original distribution after the removal totally lacks complexities higher than 2, and hence it differs vastly from the surrogate distribution (see Supplemental Material: Fig. S1). Thus, removing all HSEs significantly alters the statistical properties of the spike trains by selectively eliminating higher-order chance synchrony, which can lead to false results regarding fine temporal correlations in the spike trains.

To summarize, HSEs in massively parallel spike trains commonly contain synchrofacts, which are artificial spike synchrony events on a timescale of the data sampling frequency. Synchrofacts can only be spotted when analyzing the data on very fine timescale and show themselves as excess HSEs. Since HSEs contain not only artifacts but also genuine neuronal synchrony (as well as chance coincidences), removing all HSEs for cleaning the data is not a proper approach. Hence, we aim at discriminating artifacts from chance coincidences of real spikes, for selectively removing artifacts from the data.

Synchronous spike events with nearly identical waveforms are not of neural origin

Following the observations in the previous section, we now take a deeper look into the HSEs for features discriminating artifacts from chance spike coincidences. Whereas an HSE with a very high complexity is likely to be an artifact, we cannot rule out the possibility that it is a real spike synchrony. In fact, there is no way of discriminating between them on the description level of spikes. Thus, seeking for additional information to spot artifacts, we inspect the raw signals from which the spike events are extracted.

Here for the macaque data sets, the spike extraction was performed by thresholding the band-pass filtered raw signals in the frequency range of 500–7,500 Hz (see “Materials and Methods, Signal processing” for details). For simplicity, we henceforth refer to the band-pass filtered raw signal merely as band-passed signal. Figure 2a shows the band-passed signals of multiple channels around the timing of one HSE. We observe highly similar spike-like waveforms in multiple channels simultaneously, but with different amplitudes. Furthermore, the ongoing signal fluctuations before and after the HSE appear to be correlated across these channels. The electrodes of these channels are not necessarily in close spatial proximity, as shown in Figure 2b. To quantify this observation, we measure the distances between the electrodes of channels participating in each and every HSEs with high complexities. The distribution of the obtained distances (Fig. 2c) shows that the high complexity HSEs often involve electrodes that are far apart, up to several millimeters. In particular, in the case of macaque L, HSEs appear even across the two arrays located on opposite banks of the lunate sulcus.

To summarize, we have observed (i) sub-millisecond synchronization of spikes across channels, (ii) nearly identical waveforms of those spikes, (iii) visible correlation between band-passed signals, and (iv) partly large distances between the participating electrodes. Taken together, these observations strongly suggest that many of these events are artifacts due to electric cross-talk between channels, rather than multiple electrodes recording the same neuron or some novel form of fast neuronal synchronization.

Cross-correlation and HSE can detect artifacts

An electric cross-talk between a channel pair is expected to increase the Pearson correlation coefficient between the band-passed signals of these channels. Figure 3a shows the correlation coefficient c_i,j between the band-passed signals of channels i and j (see “Materials and Methods, Cross-correlation of band-passed signals” for details), for all channel pairs in the two data sets in a matrix form (only for i ≠ j). The distribution of the c_i,j values in the matrix is shown in Figure 3b. The coefficients are mostly larger than zero, i.e., the band-passed signals are generally positively correlated. The coefficient values vary strongly across channel pairs, with some pairs showing quite strong correlation, in both recording setups from different laboratories.

There may be multiple possible origins for such strong correlations, ranging from damages on the electrodes to electric shortcuts between channels due to various reasons, e.g., crossing of the implanted cables below the skull, dirt in the headstage connector, interference between the cables connected to the amplifier, etc. (Yu et al., 2009). As shown in Figure 2b,c, the artifacts are not limited to spatially proximal channels but extend over the whole array or even across arrays, indicating that a deficit in a local group of channels cannot be the primary reason for these correlations. Another possibility is a local deficit on the headstage connector, which might be the case for Macaque Y but not for Macaque L (see Supplemental Material: Fig. S2). Thus, no single reason can explain the observed correlations, but rather there would be multiple causes residing on different stages of the setup. Eliminating all possible causes from the setup is practically not feasible, and hence we cannot avoid having such correlations in the data.

For channel pairs with a large c_i,j value, their band-passed signals must be similar to each other. Hence, if one channel of such a pair had spike-like waveforms, the other channel should also have those at the same time. This leads to the expectation that a channel pair with a larger c_i,j should show more artifacts. To confirm this, we first count, for each channel pair, the number of HSEs containing spikes from that pair. We then define the HSE index I_i,j of channel i and j as the ratio of the obtained HSE count to the spike counts of these two channels (see “Materials and Methods, Measure of spike synchrony” for the formula). Figure 3c shows the I_i,j values obtained from the two data sets in a matrix form (only for i ≠ j). The distribution of those values in the matrix is shown in Figure 3d. We find I_i,j values ranging from 0 to 0.5, while the value expected from an independent pair of spike trains is on the order of 10⁻³ (see Supplementary Fig. S3).

We then plot I_i,j against the correlation coefficient c_ij to look for a systematic relation between them (Fig. 3e). We find that a majority of channel pairs have HSE index values close to zero, meaning that hardly any spikes are shared by these pairs. Especially the channel pairs with low correlation coefficients c_i,j, such as <0.4, have very low HSE index values. As the correlation c_i,j gets larger than 0.4, the HSE index value rapidly increases, indicating a critical correlation beyond which spikes becomes very likely to be detected in both channels. Note that an HSE index larger than zero does not immediately indicate artifacts. As shown in Figure 1d, a certain amount of HSEs can be explained by chance. However, as we mentioned before, the HSE index expected from independent spike trains is very small, on the order of 10⁻³. Hence, an HSE index excessively greater than this strongly indicates presence of artifacts.

Cross-talk model explains the emergence of artifacts

To understand the origin of the relation between the band-passed signal correlation c_i,j and the HSE index I_i,j, here we present a simple model of cross-talk between channels (see “Materials and Methods, Cross-talk model” for details). This model enables us to derive an analytic expression of the relation, which is in agreement with our experimental observations (Fig. 3e).

Our model assumes two Gaussian white noise (mean 0, variance σ²) time series x₁ and x₂ as the cross-talk-free “ground truth” signals. We model the cross-talk as a linear mixing of these two signals (Fig. 4a), such that the mixed signals s₁ and s₂ measured on channel 1 and 2, respectively, are: where 0 ≤ α ≤ 1 is the strength of the cross-talk: s₁ and s₂ are identical with α = 1, and independent with α = 0. The Pearson correlation coefficient c_1,2 between s₁ and s₂ can be written in terms of α as (see “Materials and Methods, Cross-talk model” for derivation): Thus the correlation c_1,2 increases nonlinearly with the cross-talk strength α, as shown in Figure 4b.

For the derivation of the HSE index, we need to introduce spikes to the signals. We model spikes by subtracting a value A_i, representing the spike amplitude, from the ground truth signal x_i ( ) at random time points. We parameterize A_i as , where is the signal-to-noise ratio of the spike waveform in x_i. Spikes are then extracted from the mixed signals s₁ and s₂ by thresholding, as commonly done for experimental data. We set the threshold for s_i at , where is the standard deviation of s_i and m is a multiplier (Quiroga et al., 2004). The expected numbers n₁ and n₂ of spikes detected in s₁ and s₂, respectively, are: where N₁ and N₂ are the number of spikes introduced in x₁ and x₂, respectively, and p_i,j is the probabilities that a spike from x_i is detected in s_j as a function of the parameters α, , , and m (see “Materials and Methods, Cross-talk model” for more details).

The spikes transferred across channels by the cross-talk are detected in both channels, so the sum of the numbers of such spikes, i.e., n_1,2 = p_1,2N₂ + p_2,1N₁ represents the number of artifacts in this pair of channels. We can now calculate the HSE index as I_1,2 = n_1,2/min(n₁, n₂) (see “Materials and Methods, Cross-talk model” for detailed derivations). While this model could be extended to include more realistic details such as different types of noise, spike waveforms, and so on, the simple model as described here suffices for the current purpose of deriving the HSE index analytically.

Through the dependence of p_i,j on the cross-talk strength α, the HSE index I_1,2 is dependent on α and thereby the correlation coefficient c_1,2 (Fig. 4c). Interestingly, hardly any artifacts are observed in the model for α < 0.2 or c_1,2 < 0.4, suggesting that some degree of cross-talk might be tolerable in recording systems. However, for higher α or c_1,2, the HSE index rapidly increases in the form of a sigmoid, depending on the signal-to-noise ratio r^SN. This is in good agreement with what we observed in our macaque datasets (Fig. 3e).

In the real data, the signal of one channel is typically correlated not only with another channel but with multiple other channels at once, as seen in Figure 3. The model can be extended to represent such data by introducing multiple channels (let N_ch be the number of channels) and expanding the cross-talk strength parameter α to a cross-talk strength matrix α_ij (1 ≤ i ≤ N_ch, 1 ≤ j ≤ N_ch), which represents the cross-talk strength between channel i and j (see “Materials and Methods, Cross-talk model” for more details).

Here we demonstrate that the extended multi-channel model can simulate artifacts and correlated signal fluctuations as seen in our experimental data (Fig. 2a), given a proper cross-talk strength matrix. We consider a model of 100 channels and a cross-talk strength matrix as shown in Figure 4d. The diagonal elements of the matrix are all set to 1 (to be consistent with the channel-pair model shown above), and the remaining elements are defined such that the data contains three groups of 20 channels (i.e., channels 0–19, channels 40–59, and channels 80–99) with strong cross-talk within each group. The ground truth signal of each of the 100 channels is generated independently as a Gaussian white noise time series with a Poisson spike train (20 Hz firing rate, 0.3 ms spike width) injected (see “Materials and Methods, Cross-talk model” for more details of the simulation), as shown in Figure 4e, left. Since the injected spike trains are independent between channels, the complexity distribution of the spikes in the ground truth signals is just as expected from chance spike coincidences between channels (Fig. 4f, blue). From these ground truth signals, mixed signals are generated according to the cross-talk strength matrix (Fig. 4e, right). Strong correlations in the signals among the channels within each cross-talking channel group (black traces) are evident, compared to the correlations among the channels not participating in those groups (gray traces). As a consequence, when there is a chance spike coincidence in the ground truth signals of a cross-talking channel group (e.g., at the middle of Fig. 4e, left, for the bottom channels), the strong cross-talk within the group replicates those spikes into the other channels of the group at the same timing, resulting in a formation of an artifact (Fig. 4e, right, for the bottom channels). Accordingly, the complexity distribution of the mixed signals exhibits a heavy tail of high order complexities (Fig. 4f, gray), which resembles the complexity distribution of the real data (Fig. 1d).

ZCA whitening effectively removes artifacts

From this model we now see that the correlations between the signals lead to emergence of artifacts. Conversely, the model poses a possibility that if one could properly decorrelate the measured signals, the ground truth signals could be recovered and the artifacts would be dissolved. Since our model assumes a linear mixture of the ground truth signals via symmetric cross-talks, the measured signals can be decorrelated by whitening them based on their covariance matrix, a method called ZCA whitening (Bell and Sejnowski, 1997; Kessy et al., 2018).

In Fig. 5a we illustrate an example application of ZCA whitening to simulated data of the model. Similarly to the example shown in the previous figure (Fig. 4e), a chance spike coincidence in the ground truth signal (Fig. 5a, left) develops to an artifact in the mixed signal (Fig. 5a, center) due to cross-talk. After an application of ZCA whitening to the mixed data, the ground truth signals are perfectly recovered for all channels (Fig. 5a, right), and the spikes that compose the artifact are removed except for the two that are the original chance coincidence spikes that caused the artifact. All the other spikes outside the artifact are also kept intact.

To quantify the performance of ZCA whitening on linearly mixed signals, we statistically compared the ground truth and the ZCA whitened signals obtained from simulations of our cross-talk model with a simple setting: 10 channels with all-to-all, uniform cross-talk and a common, homogeneous firing rate. We systematically varied the uniform cross-talk strength α, the duration T of the signals, and the common firing rate λ, and examined how the Pearson correlation coefficient between the whitened signals and their respective ground truth signal varies depending on these parameters. As summarized in Figure 5b, the correlation coefficient is almost 1 for the signal durations longer than 10 s, while it drops toward shorter durations, but never falls below 0.9995 (median over 10 channels times 10 simulations per case). The firing rate and the cross-talk strength do not seem to affect the correlation. Thus, for experimentally meaningful data durations (>10 s), the ZCA whitened signals are practically identical to the ground truth signals. This means that spikes can be detected in the whitened signals with the same waveforms and timings as in the ground truth signals. Furthermore, even for shorter durations, we confirmed that all spikes in the ground truth signals are detected in the whitened signal with exactly the same spike times.

In fact, ZCA whitening is exactly the inverse operation of the linear mixture performed by the proposed cross-talk model, as long as the covariance matrix of the mixed signals is correctly estimated. The declined performance of the ZCA whitening seen for shorter data durations is due to inaccurate estimation of the covariance matrix because of a smaller sample size for the estimation. Thus, given enough data, ZCA whitening should perform almost perfectly to recover the ground truth signals for any mixed signals where the correlations are generated according to the cross-talk model.

Here we note that, if temporal correlations across the spike trains, e.g., excess spike coincidences reflecting some physiological processes taking place in the recorded neuronal population (Grün et al., 2002a, 2002b), exist in the ground truth signals, they are not removed by the ZCA whitening and hence they correctly remain in the whitened signals. To illustrate this, we show in Figure 5c an application of the ZCA whitening to a simulation of our cross-talk model where spike coincidences mimicking physiological correlations in spike timing are manually introduced. The “physiological” coincidences in the ground truth signals are exaggerated by the cross-talk, generating a lot of copied artificial spikes in the mixed signals. The ZCA whitening removes all these copied spikes, while all the original spikes, including the manually introduced coincident spikes, remain the same in the whitened signals as in the ground truth signals. Only if the rate of the manually introduced spike coincidences is un-physiologically high (for an estimate of physiological coincidence rate, e.g., Maldonado et al., 2008), some of the coincident spikes are missed in the whitened signals (Fig. 5d). Thus, the ZCA whitening does not destroy physiological correlations in spike timing that may exist in the ground truth signals.

Now we apply ZCA whitening to our experimental data and see if this method can remove the artifacts. Figure 6b shows the complexity distributions for the macaque Y (top right) and macaque L (bottom right) data with application of ZCA whitening [and for comparison, the corresponding distributions without ZCA whitening (Fig. 6b, left), which are the same data as in Fig. 1d]. With ZCA whitening, the complexity distribution matches quite tightly to the distribution from surrogate data, i.e., independent data, with a considerable reduction of high order complexities seen in the complexity distribution without ZCA whitening.

To quantify the effects of ZCA whitening on single channels, we introduce a measure evaluating the participation of each channel in HSEs. We term this as the global HSE index I_i for channel i, defined as the ratio of the number of HSEs including channel i to the number of all spikes on channel i (see “Materials and Methods, Measure of spike synchrony” for the formula). The HSE index I_i,j that we have introduced before for a pair of channels i and j is henceforth referred to as pairwise HSE index. Furthermore, we also introduce a channel-wise measure of signal correlations, i.e., the highest correlation coefficient for each channel i to all other channels, denoted as C_i for channel i and defined as , where c_i,j denotes the correlation coefficient between the band-passed signals of channel i and j. Figure 6c shows, for data with (right) or without (left) ZCA whitening, the global HSE indices I_i plotted against their respective highest correlation coefficients C_i for all channels in macaque Y and L data sets. In the case without ZCA whitening (Fig. 6c, left), as expected from the model, the channels with high C_i values are indeed also the channels with high I_i values, while the channels with C_i below 0.4 show rather small I_i. On the other hand, in the case with ZCA whitening (Fig. 6c, right), all channels show almost zero correlations, which is a direct consequence of the whitening, and also show considerably lower I_i values, indicating an overall suppression of HSEs by ZCA whitening.

Figure 6d shows the statistics of the two data sets with (y-axes) and without (x-axes) ZCA whitening. Figure 6d (top) compares the global HSE index I_i with and without ZCA whitening. All channels show a reduction in I_i with ZCA whitening. This reduction is achieved without reduction in the firing rates of the spike trains in a majority of channels, as evidenced by the comparison of channel-wise firing rate with and without ZCA whitening (Fig. 6d, middle). Most channels show similar firing rates with and without ZCA whitening, indicating that ZCA whitening removes only a tiny fraction of spikes, which are most likely the ones copied by the cross-talk between channels. Similarly the spike waveform SNR does not change from after the whitening: the signal is largely unchanged (see Supplemental Material: Fig. S5).

High-density probe recordings also contain artifacts

So far we have focused on the occurrence and removal of artifacts in recordings with Utah electrode arrays, where the inter-electrode distance is at least . Nowadays more and more laboratories start to use high-density electrodes, such as Neuropixels probes (Jun et al., 2017), with inter-electrode distances as small as at minimum. With such a high electrode density, a spike of one neuron can and is intended to be recorded at multiple neighboring electrodes, which causes an apparent synchronous spike event among the respective channels. Such synchronous spikes (of different shapes) are used to sort out the contributions of single units by use of a waveform template considering signals from multiple neighboring channels, as is done in, e.g., the Kilosort sorting algorithm (Pachitariu et al., 2024). However, this does not necessarily mean that spike trains sorted in such a way are free from artifacts. They can still contain artifacts due to cross-talks between distant channels pairs, or any other causes that introduce strong correlations between signals at distant channels. In this section we focus on spike train data from a Neuropixels recording sorted by Kilosort to show that these data also contain artifacts and our ZCA whitening-based approach presented above can be used to remove them.

The data we consider here were recorded with a Neuropixels 1.0 probe (see “Materials and Methods”) in an awake mouse during visual stimulation. The probe was inserted through the primary visual cortex (V1; Fig. 7a), with the tip of the probe reaching the SC. The band-pass filtered raw signals appear largely similar across all channels (Fig. 7b, left). This is primarily because all recorded channels are relative to the external reference electrode (in this recording placed on the dura above the cerebellum) and are therefore equally affected by signal fluctuations at the reference. Consequently, the signals can show very strong correlations between channels, even at distant channel pairs (Fig. 7c, top).

A standard signal pre-processing step is therefore to first remove signals from electrodes outside of the brain and then to calculate the median of signals across all remaining channels at each time point. By subtracting this median signal from each channel—commonly known as median removal or common average referencing—common noise patterns that are shared across all channels are effectively removed (for more details see “Materials and Methods”). The signals after this median removal show more variability across channels (Fig. 7b, center) and generally contain much less correlation between channels (Fig. 7c, middle). However, there are still strong correlations in channel pairs with very short distances, as expected from the conductivity of the brain tissue. Besides that, there are clear negative correlations between channels separated by more than , which can be recognized also in the signal traces in Figure 7c (middle panel) for top and bottom channels showing inverted waveforms. These negative correlations seem not to be of physiological origin, because the bottom electrodes are in another brain region than the top electrodes (namely, SC and V1, respectively). One may additionally notice oscillatory modulations of the correlation coefficient as a function of the inter-electrode distance. In the scope of the present study, we do not discuss these oscillations further, but we suspect that they most likely stem from the hardware or get introduced by the correction for the time shifts introduced by the multiplexing during the recording (see “Materials and Methods”). The time-shift correction is a necessary pre-processing procedure to obtain properly time-aligned signals from the raw recordings.

ZCA whitening on the median-removed signals almost totally gets rid of the positive and negative correlations (Fig. 7c, bottom). Kilosort sorting applied to the median-removed and ZCA whitened signals results in a considerable difference in terms of the complexity distribution of the obtained spike trains (Fig. 7d): the distribution for the median-removed signal (Fig. 7d, “BP + MR”) shows a heavy tail of high order complexities up to 20, while the distribution for the ZCA whitened signal (Fig. 7d, “BP + MR + ZW”) shows a faster decay with a maximum complexity of 7.

We need to note that ZCA whitening is already included in the typical and recommended workflow of the Kilosort sorting algorithm as a signal pre-processing step, which we intentionally disabled when we compute the complexity distribution for the median-removed signal in Figure 7d. The ZCA whitening in Kilosort is by default applied to the 32 channels closest to the target channel. In this manner, the signals in the channels separated by more than 32 channels are not decorrelated to each other, i.e., the negative correlations between distant channel pairs we found in the median-removed signals are not removed. This could lead to artifacts between sorted single units in the top and the bottom channels, and in fact, the complexity distribution for the median-removed signals sorted with the Kilosort default settings (Fig. 7d, “kilosort default”) shows more high order complexities than our ZCA whitened signals, where ZCA was applied to all channels simultaneously (Fig. 7d, “BP + MR + ZW”). To further examine the degree and location of artifacts in these data, we compute the global HSE index I_i for each sorted single unit and plot it against the position of that unit on the probe (Fig. 7e). In the case without ZCA whitening (Fig. 7e, left), the units around the top and the bottom of the probe stand out with particularly high I_i, as we expected from the negative correlations between top and bottom channels. These correlations cause Kilosort to extract spikes in these channels at the same time point just with opposite signs. In the case with ZCA whitening with Kilosort default settings (Fig. 7e, center), much less number of spikes are detected and hence units naturally show lower I_i values, but the units around the top and the bottom of the probe still exhibit relatively high I_i values compared to the other units. This is a reflection of the negative correlations that still remain in the signals by the spatially restricted application of ZCA (32 channels in parallel) by Kilosort. In the case with ZCA whitening applied to all channels altogether (Fig. 7e, right), a similar number of spikes and units are detected as in the previous case, whereas the units around the top and the bottom of the probe do not show particularly high I_i values, indicating that artifacts including these units are removed. Although whitening across all channels slightly reduced the overall waveform amplitude, the SNR of the detected units remained essentially unchanged (Supplemental Material: Fig. S6).

With these results, we conclude that, although the internal cause of artifacts is probably different from the case using Utah array(s), Neuropixel recordings can contain artifacts (although of different appearance) and ZCA whitening applied to all channels together effectively removes those artifacts from the data. We need to note that the use of ZCA whitening here cannot be considered as resolving the linear mixture of the signals according to our cross-talk model, since the background fluctuations in neighboring electrodes cannot be considered as independent as the model assumes, due to the close proximity of the electrodes in the Neuropixels probe. What we want to stress here is that HSEs do exist in high-density probe recordings even after spike sorting, and that the application of ZCA whitening to 32 neighboring channels only, which is the default setting of the Kilosort spike sorter and widely used in many studies, is not capable of eliminating these HSEs. We have shown that the ZCA whitening to all channels at once removes HSEs in the Neuropixels recording, with approximately the same number of spikes and sorted units as obtained with the default setting. Hence, we recommend to apply ZCA whitening always to all channels at once to obtain spike data free from the synchronous artifacts, even though the application to high-density probe recordings is not validated by the cross-talk model as in the case for Utah array recordings.