Multi electrode recordings of neuronal activity provide an overwhelming amount of data that is often difficult to analyze and interpret. Although various methods exist for treating multielectrode datasets quantitatively, there is a particularly prominent lack of techniques that enable a quick visual exploration of such datasets. Here, by using Kohonen self-organizing maps, we propose a simple technique that allows for the representation of multiple spike trains through a sequence of color-coded population activity vectors. When multiple color sequences are grouped according to a certain criterion, e.g., by stimulation condition or recording time, one can inspect an entire dataset visually and extract quickly information about the identity, stimulus-locking and temporal distribution of multi-neuron activity patterns. Color sequences can be computed on various time scales revealing different aspects of the temporal dynamics and can emphasize high-order correlation patterns that are not detectable with pairwise techniques. Furthermore, this technique is useful for determining the stability of neuronal responses during a recording session. Due to its simplicity and reliance on perceptual grouping, the method is useful for both quick on-line visualization of incoming data and for more detailed post hoc analyses.
The highly distributed architecture of the cerebral cortex enables neurons to constantly process patterns of spikes originating from multiple presynaptic sources (Douglas and Martin 2004). Therefore a proper characterization of brain dynamics has to rely on parallel recordings that enable the simultaneous sampling of the activity of many neurons (Brown et al. 2004; Buzsáki 2004; Singer 1999). Such recordings, however, generate large amounts of data that are multidimensional in nature and are thus difficult to analyze. To handle such datasets, advanced quantitative methods are being developed (Abeles and Gerstein 1988; Chapin and Nicolelis 1999; Grün et al. 2002; Laubach et al. 1999; Perkel and Bullock 1968; Pipa et al. 2008; Puchalla et al. 2005; Schneidman et al. 2006; Yu et al. 2008).
Progress has also been made in developing techniques of data visualization that enable the experimenter to obtain a quick overview of multidimensional signals. These methods surpass the visualization of response properties by classical spike rastergrams or histograms, which, due to the arbitrary ordering of spike trains, provide little help in detecting multidimensional spike patterns (Toups and Tiesinga 2006). Techniques for visualization of multineuron activity include the “gravitational” method by Gerstein et al. (1985), which was designed for detecting correlations between multiple neurons, and a class of methods that visualize neuronal activity through a trajectory plotted in a reduced two- or three-dimensional state space (Bathellier et al. 2008; Brown et al. 2005; Friedrich and Laurent 2001; Galán et al. 2004). Two other methods enable the detection of groups of correlated neurons: the “correlation grid” (Stuart et al. 2004) and plots using “parallel coordinates” (Inselberg and Dimsdale 1990; Stuart et al. 2002). However, these latter ones do not provide means for visualizing the temporal dynamics of responses but offer only static representations of correlations within neuronal groups. Also most of the existing visualization methods are limited to a relatively small number of neurons and thus cannot deal with high-dimensional firing patterns.
Here we propose a method that uses colors to represent the identity of multineuronal activity patterns and enables the quick visualization of large multielectrode datasets, while preserving information about the structure (firing combination) of neuronal activity patterns. The method also allows for the inspection of activity patterns on various time scales and enables the estimation of the stationarity of cortical responses over time. Importantly, we will focus here only on the aspect of data visualization but will also discuss later some aspects of data classification, which is a different but related issue.
Experimental procedures and recording
Data were recorded from area 17 of four anesthetized and paralyzed adult cats bred in the facilities of the Max-Planck Institute for Brain Research. Experiments were carried out in accordance with the European Communities Council Directive of 24 November 1986 (86/609/EEC). Anesthesia was induced with ketamine (Ketanest, Parke-Davis, 10 mg/kg im) and xylazine (Rompun, Bayer, 2 mg/kg im) and maintained with a mixture of 70% N2O-30% O2 supplemented with halothane (0.5–1.0%). After tracheotomy, the animals were placed in a stereotactic frame. A craniotomy was performed, and the skull was cemented to a metal rod. After completion of all surgical procedures, the ear and eye bars were removed, and the halothane level was reduced to 0.4–0.6%. After assuring that the level of anesthesia was stable and sufficiently deep to prevent any vegetative reactions to somatic stimulation, the animals were paralyzed with pancuronium bromide (Pancuronium, Organon, 0.15 mg·kg−1·h−1). Glucose and electrolytes were supplemented intravenously and through a gastric catheter. The end-tidal CO2 and rectal temperature were kept in the range of 3–4% and 37–38°C, respectively. Stimuli were presented binocularly on a 21-in computer screen (HITACHI CM813ET) with 100-Hz refresh rate. To obtain binocular fusion, the optical axes of the two eyes were first determined by mapping the borders of the respective receptive fields and then aligned on the computer screen with adjustable prisms placed in front of one eye. All the experiments were conducted according to the guidelines of the Society for Neuroscience and German law for the protection of animals, approved by the local government's ethical committee, and overseen by a veterinarian.
Data were recorded by inserting multiple silicon-based multi-electrode probes (16 channels per electrode) supplied by the Center for Neural Communication Technology at the University of Michigan (Michigan probes). Each probe consisted of four 3-mm-long shanks that were separated by 200 μm and contained four electrode contacts each (1,250 μm2 area, 0.3–0.5 MΩ impedance at 1,000 Hz, intercontact distance 200 μm). Signals were amplified 10,000× and filtered between 500 Hz and 3.5 kHz and between 1 and 100 Hz for extracting multi-unit (MU) activity and local field potentials (LFPs), respectively. The waveforms of detected spikes were recorded for a 1.2-ms duration, which allowed the later application of off-line spike-sorting techniques to extract single units (SUs). Spike sorting was carried out with a custom-made software relying on principal component analysis of the waveforms. The software for visual stimulation was a combination of custom-made programs and a stimulation tool, ActiveSTIM (www.ActiveSTIM.com).
DATASET WITH DRIFTING SINUSOIDAL GRATING STIMULI (col05-e08).
Sinusoidal gratings moving in 12 directions in steps of 30° were presented in trials of 4,800-ms duration (1,000-ms spontaneous activity, 3,500-ms stimulus, 300-ms off response). Gratings spanned 12° of visual angle, had a spatial frequency of 2.4° per cycle, and drifted with a speed of 2°/s. Each direction was presented 20 times in a randomized order leading to the total of 240 presentations (trials). Dataset col05-e08 consisted of 26 simultaneously recorded SUs.
DATASET WITH PLAID AND BAR STIMULI (col08-e27).
Eight stimuli were used including: two simultaneously moving bars, in one case with a dot at the intersection, two single bars moving at different speeds, two plaids with pattern or component motion, a grating, and a moving dot (see Fig. 3). The plaids and the grating had a rectangular change in luminance. Bars spanned 7° and gratings and plaids 12° of visual angle. Trials were 6,500 ms long with stimuli presented between 1,000 and 6,000 ms, and 20 times each. Analyses were performed on 32 simultaneously recorded, unsorted MUs.
DATASET WITH NATURAL STIMULI (cer01-a50).
Three movies with natural images were presented to the cat (1 recorded by the authors, 2 extracted from a documentary “The Greatest Places” by the Science Museum of Minnesota). The movies contained indoor and outdoor scenes with various image statistics (slow moving, fast moving, dark, light, etc). The movies were presented across the entire screen (41 × 31° of visual angle) with a resolution of 800 × 600 pixels. Each movie was 28 s long and was presented 20 times. Analyses were performed on 34 simultaneously recorded SUs.
DATASETS WITH CENTER-SURROUND STIMULI (col11-b44 and col11-b68).
Sinusoidal gratings of three different sizes (small, medium, and large) and two orientations (horizontal and vertical) were presented individually or superimposed. Gratings spanned visual angles of 7, 14, and 21°, had a spatial frequency of 1° per grating cycle and were drifted at a speed of 1.5°/s, orthogonal to their orientation and in one direction only. Stimuli included six individual gratings, four superimposed gratings with a small central grating surrounded by an orthogonal medium or large grating, and four superimposed gratings consisting of a small grating separated by a gray ring (3.5° wide) from a surrounding large grating of identical or orthogonal orientation. The resulting 14 stimuli were randomly presented 20 times each, leading to a total number of 280 trials. Trials were 6,000 ms long with the stimulus presented between 1,000 and 5,000 ms. After spike-sorting, from a total of 66 SUs only 12 SUs could be confidently tracked (matched) in both datasets and further used in the analyses.
DATASET WITH HIGH SPATIAL-FREQUENCY SINUSOIDAL GRATING STIMULI (cer01-a45).
These stimuli were identical to those used in col08-e05 with the exception that gratings had a higher spatial frequency of 1.1° per cycle.
The visualization technique includes three steps: low-pass filtering spike trains, Kohonen mapping, and, finally, building color sequences (Fig. 1).
STEP 1: LOW-PASS FILTERING SPIKE TRAINS.
To handle multiple spike trains, one must first find a robust multidimensional representation. Most previous studies have applied binarization in combination with binning to define multidimensional activation vectors (Grün et al. 2002). While being convenient for defining patterns as binary vectors (Schneidman et al. 2006), binning prevents a flexible manipulation of the time scale because, for large bin sizes, multiple spikes fall into the same bin. Replacing binarization with instantaneous firing rate (spike count per bin) does not solve the problem because, especially for large bin sizes, one cannot distinguish different firing patterns within the bin. Instead of binning, spike trains can be first convolved with specialized kernels such as Gaussian (Eskandar et al. 1992; Heller et al. 1995; Richmond et al. 1990), sinusoidal (Sameshima and Baccalá 1999), or exponentially decaying kernels (Baker and Gerstein 2000; Häusler and Maass 2007). We used the latter because such kernels are essentially causal low-pass filters (i.e., the values of the signal depend only on the preceding spikes). The convolution preserves information about all spikes and allows one to control the time scale of interest by manipulating the decay (integration) time constant (τ). For each neuron i, a continuous signal, called activation, ai(t) was obtained using the formula (1) where ai(t) is the activation corresponding to neuron i at time t, and τ is the decay (integration) time constant. For small time constants (τ = 1–5 ms), one extracts synchronous patterns of spikes or joint-spike events. For large time constants (τ >100 ms), one extracts firing-rate patterns. Unless specified otherwise, an integration time constant of 20 ms was used because it is in the physiological range of membrane time constants of cortical neurons (Kasper et al. 1994; Lefort et al. 2009; Spruston and Johnston 1992).
Low-pass filtering of spike trains (Fig. 1, steps 1 and 2) yielded multiple continuous traces which were then sampled (Sameshima and Baccalá 1999) with a frequency of 1 kHz, and for each time step t (with millisecond resolution), an activity profile, called activity vector, was defined as (2) where n is the number of analyzed neurons. A trial was then described as a succession of activity vectors. Note that an activity vector also contains information about the activity in the recent past of the corresponding neurons (Fig. 1, step 2).
STEP 2: KOHONEN MAPPING.
After integrating the original spike trains using exponentially-decaying kernels, the activity vectors sampled at each millisecond were presented, in random order, to a three-dimensional (3D) Kohonen map (3DKM) for learning (Fig. 1, step 3). 3DKM is an extension of the classical 2D Kohonen self-organizing map (Kohonen 2001) simply by adding one dimension (3DKM of size N is defined on an N × N × N lattice). Each element in the lattice contains a vector of dimension equal to the dimensionality of the input space, termed model vector. At each step k of the learning algorithm, the map learned an activity vector (AVk) as follows: first, the most similar model vector to AVk was found in the map and called best-matching unit (BMU); second, the BMU and its neighbors (Kohonen 2001) were altered to make them more similar to AVk. The amount of change and the radius of the neighborhood were given by two monotonically decreasing functions: L(k) and R(k), respectively (3) where L(k) is the learning rate, modulating the degree to which model vectors were changed at each training step, k. L0 and LM are initial and final learning rates. We used L0 = 1 and LM = 0.01. The total number of training steps is denoted by M (4) where round denotes the rounding to the nearest integer, R(k) specifies the neighborhood size around the BMU within which elements were allowed to learn at step k. R0 is the initial radius of the neighborhood. g is the percentage of M after which R becomes 0 (only the BMU is modified for R = 0). We used R0 = N/2 and g = 66 (66% of steps were used to establish the topology of the map and the remaining 34% of the steps to fine tune the representation of activity vectors in the map, i.e., only the BMU was changed).
Within the above-defined neighborhood, model vectors further away from the BMU changed less than the ones closer to it by multiplying the learning rate with a 3D Gaussian envelope with a SD of R(k)/3 (5) where MVk[x,y,z] is a model vector, at step k of the training, located within the neighborhood of the BMU [distance from BMU ≤ R(k)] at position (x,y,z) in the 3D lattice. (xBMU, yBMU, zBMU) is the position of the BMU in the 3D lattice. AVk is the activity vector that is learned at step k, L(k) and R(k) are, respectively, the learning rate and the size of the neighborhood at step k.
The map was initialized with values of 0. The convergence of the learning algorithm depends on the total number of steps (M). Typically ∼100,000 steps are recommended, but in some cases, where “fast learning” is required, even 10,000 steps may be sufficient (Kohonen 2001). In the data investigated here, depending on the number and length of trials, the number of activity vectors varied across recording sessions (ranging from 1,040,000 to 1,680,000). The number of steps used in our analyses was typically three times the number of activity vectors, i.e., each activity vector was learned three times by the map. Activity vectors were presented in random order to the 3DKM to avoid learning temporal dependencies and to maximize the fidelity with which the clustering (model vectors) represented the activity vector space. In most cases, we used maps that included 1,000 points (clusters) in the 3D lattice (N = 10), ensuring a fine grain representation of the activity vector space. Thus for N = 10 there were 1,000 patterns (model vectors) that described a given dataset.
All activity vectors of a given dataset were learned using the 3DKM (Fig. 1, step 3). The Kohonen learning is essentially an ordered clustering, i.e., clustering and ordered mapping on the lattice (Kohonen 2001). After learning, the lattice holds model vectors that are in fact clusters approximating the input space. They represent stereotypically appearing spatiotemporal activity patterns spanning a finite time window determined by the integration time constant. For simplicity, model vectors will be called patterns.
During training, the approximation error, i.e., the error between the activity vectors of the dataset and their corresponding model vectors, was computed as (6) where E is the approximation error, is the Euclidean distance between the activity vector and its corresponding BMU from the map, and NA is the total number of activity vectors for a dataset. The approximation error drops dramatically in the first steps and then converges to a low value toward the end of the training (Supplemental Fig. S1A).1
It can be checked whether model vectors represent faithfully the activity vectors also by computing the distribution of Pearson coefficients of correlation between the two. In the example provided in Supplemental Fig. S1B, 82% of the coefficients were >0.8, i.e., model vectors approximated activity vectors with a high degree of precision.
The computation time required for training the Kohonen map is fast to moderate, depending on the amount of data and the size of the 3DKM. For example, the map used to produce Fig. 2B was computed in a few minutes (<15 min) on a personal computer with average performance. With modern parallel computing, the time to compute the map can be considerably reduced.
STEP 3: BUILDING COLOR SEQUENCES.
For visualization, each dimension of the 3DKM lattice was considered to be one color dimension in the RGB space (red, green, and blue; Fig. 1, step 3). A model vector was then assigned a corresponding color depending on its position in the lattice. The coloring strategy took advantage of the fact that, during learning, the Kohonen algorithm maps similar model vectors to nearby locations (i.e., ordered mapping). Due to the spatial (Euclidean) proximity of similar model vectors in the lattice (Kohonen 2001), these were assigned similar colors. In addition, the matching cubic topologies of the 3DKM lattice and the RGB space ensured that, for each analysis, the entire spectrum of colors was used to represent the model vectors.
A trial could then be visualized as a sequence of colors, whereby each activity vector was labeled by the color of its corresponding model vector. This was done as follows: for each time bin, the activity vector was used to find its closest model vector in the 3DKM cube (its stereotypical model). At the respective time bin in the trial the color of the model vector was marked through a vertical color band (a painted rectangle with thickness = 1 ms; Fig. 1, step 4). A color sequence represents the activity of all neurons along a single recorded trial.
Collections of color sequences can be presented in an adjacent and time-aligned manner—usually one color sequence is shown below the other (e.g., Fig. 2B). For the detection of regularities, our presentation method relies on the Gestalt grouping principles of the human visual system. The regions of similar or near-similar colors blend into an easily detectable larger one-color region (e.g., a blob; e.g., Figs. 2B and 3, A and C). Thus it is easy to detect the regions in which similar patterns (colors) occur consecutively along the trial and/or consistently at the same time points across different trials. The choice of the order in which the trials will be aligned (i.e., the grouping order) determines to a large degree the type of information about neuronal dynamics that can be visually extracted.
Importantly, one can also chose to inspect individual color sequences carefully, tracking given color bands corresponding to some specific patterns and thus detecting the latter even without relying on Gestalt grouping. In such cases, however, a global picture of the whole dataset is not possible to attain.
In the training process of the Kohonen map, activity vectors are presented in random order, and therefore different runs of the algorithm yield maps with different structure. The algorithm converges such that learning is not affected by the random presentation of the samples, e.g., across 10 different runs for a drifting grating dataset, the coefficient of variation for the learning error was 0.000625 (Supplemental Fig. S1A). However, because in different runs maps will have different structure, the assignment of colors to individual patterns is likely different. Therefore colors obtained with one map cannot be directly matched with those obtained with other maps. A color assignment could be more salient than the other. However, saliency can be improved by rotations of the Kohonen cube (i.e., changing the assignment of the red, green, and blue dimensions relative to the 3DKM).
The mapping on the Kohonen cube ensures that all colors are used to represent a dataset. However, the color similarity between model vectors cannot be appreciated with respect to the absolute distance between them. Depending on the structure of a dataset, a map might assign more similar colors to some vectors, whereas for another dataset, it might assign more different colors to the same vectors. In other words, the spread of the data might be relevant to understanding how color similarity should be interpreted with respect to absolute distance between vectors. A possibility to estimate the spread of vectors is to compute the distribution of all pairwise Euclidean distances between model vectors in the Kohonen map. In addition, the distances have to be normalized to the dimensionality of the vectors (i.e., neuron number) to have results comparable across datasets with different number of neurons. Example distributions for two different datasets are shown in Supplemental Fig. S1C. For a dataset recorded with drifting sinusoidal gratings (see Datasets), the distribution is shifted toward lower values, indicating that model vectors are rather close to each other. On the contrary, for a dataset recorded with plaid and bar stimuli (see Datasets), the distribution is shifted toward larger values, indicating that model vectors are farther apart from each other (Supplemental Fig. S1C).
To study the ability of the method to detect patterns evolving on different time scales, we generated two artificial datasets of random nonstimulus related spiking and embedded/inserted within them additional spikes that represented stimulus-specific activity. The datasets consisted of 25 spike trains, simulated as homogenous Poisson processes (Cox and Isham 1980) with a rate of 10 Hz. These were segmented in 3,000-ms-long trials, which were distributed among 10 hypothetical stimulation conditions, 20 trials per condition. In one dataset, we inserted stimulus specific joint-spike events (JSEs) (Pipa et al. 2008), consisting of synchronized spikes across five neurons. In the other dataset we added stimulus-specific rate covariations (RCs), consisting of simultaneous rate increases, lasting 100 ms, in three neurons. In this later case, firing rates increased within the first 25 ms from 10 to 50 Hz and stayed then constant for the next 50 ms, falling back to 10 Hz during the last 25 ms of the RC (Lerchner et al. 2006). For each stimulus condition, we randomly assigned different combinations of neurons participating in JSEs and RCs. JSE and RC events occurred two to three times per trial and always at the same time point for a given stimulation condition.
Pattern coefficient of variation
We defined a coefficient of variation for patterns [PCVj(t)], which gives a normalized, time resolved measure of the variability of all patterns across trials of stimulus j. This measure can be used to quantify how variable is the occurrence of patterns across different trials of the same stimulus. A high PCV indicates unreliable responses across trials, whereas a low PCV indicates that the same patterns are produced on repeated presentations of the same stimulus. For a single pattern p (model vector MVp), the classic coefficient of variation (CVp,j) in a sliding window [t − h,t + h], across trials of stimulus j, was defined as the ratio between the SD of the pattern counts and the average pattern counts in the window, across trials (7) (8) (9) where SDp,j is the SD of counts for pattern p across trials of stimulus j, t represents the time in the trial, h is half the size of the sliding window (we used h = 50 ms), rl,p,j(t − h,t + h) is the count of pattern p in window [t − h,t + h] of trial l belonging to stimulus j, and Tj is the number of trials for stimulus j.
Because multiple patterns can occur in a given window, we defined the average coefficient of variation for a stimulus j, in the sliding window centered at t, as (10) where 2h + 1 is the size of the window, r̅p̅,̅j̅(̅t̅)̅ is the average count of pattern p in window [t − h,t + h] (Eq. 8) for stimulus j, and Pj is the number of distinct patterns that occur in the window across trials of stimulus j.
Patterns within a window can have different degrees of similarity to each other. To consider this, we defined an average pair wise distance between patterns, as follows (11) where ∥MVpu − MVpv∥ is the Euclidean distance between patterns pu and pv, and are average pattern counts of patterns pu and pv in window [t − h,t + h] across trials of stimulus j, and Pj is the number of distinct patterns that occur in the window. If D̄j is small, these patterns are highly similar.
Finally, the coefficient of variation for patterns (PCVj) corresponding to stimulus j, was (12) PCVj(t) is 0 when only one pattern occurs in the window at t across all trials of stimulus j, and it can increase to different values, depending on the structure of the dataset. Note that PCV cannot be compared with the classic CV due to the weighing with the distance, which can take arbitrary values. In the analysis shown in Fig. 5, we averaged PCVj(t) across all stimuli j. It is important to emphasize that each PCVj was computed per stimulus and only then taken into the average. Thus a small average PCV does not mean that patterns corresponding to different stimuli were similar but rather that patterns were similar across multiple trials belonging to the same stimulus.
Pattern specificity index
To determine how specifically a pattern p occurs for a stimulus s, we defined the pattern specificity index (PSI) as the fraction of occurrences of pattern p for stimulus s with respect to the total number of occurrences of pattern p across all stimuli (13) The PSI of pattern p for stimulus s has a value of 1 when the pattern occurs exclusively for that stimulus and a value of 0 when the pattern never occurs for that stimulus.
Visualization of large datasets
We first analyzed a dataset recorded from cat area 17, evoked with drifting sinusoidal gratings and consisting of 26 simultaneously recorded SUs (see methods, Dataset with drifting sinusoidal grating stimuli). When the activity of all 26 neurons is represented as spike rastergrams grouped by trials belonging to four different stimuli (Fig. 2A), one can detect only the onset (on) and offset (off) responses and a certain degree of rate modulations between these two events. Spike rastergrams do not easily reveal the specific firing patterns (spiking combinations) that are evoked by different stimuli. By contrast, when color sequences are computed on the same data and are grouped according to the corresponding stimulus, specific activity patterns and their occurrence in time become visible (Fig. 2B). Color sequences reveal spontaneous activity (first 1,000 ms), onset and offset of stimuli (1,000 and 4,500 ms, respectively), and modulation of responses as the gratings pass the receptive fields (RFs) of neurons. Unlike spike rastergrams, color sequences show stimulus-specific patterns as certain colors predominate in the responses to certain stimuli. The expression of these colors is modulated by the passage of the grating over the receptive fields, reminiscent of phase-sensitive “simple” cells that participate in the multidimensional firing patterns (Hubel and Wiesel 1959). Such modulatorly effects are not visible in spike rastergrams when phase-insensitive “complex” cells (Hubel and Wiesel 1962) dominate the activity (Fig. 2A, 1st stimulus) because in such cases spikes corresponding to “simple” cell activity are embedded in the sustained activity of “complex” cells. However, because activity patterns represented by color sequences consist of both simple and complex cells, these patterns will be always modulated due to the contribution of simple cells. The color sequences can therefore always reveal the modulation of multineuronal firing patterns caused by simple cell activations (Fig. 2B, 1st stimulus). For the complete set of color sequences, across all 12 directions of the drifting grating, see Supplemental Fig. S2.
Different colors in the color sequences represent different firing patterns of neurons. In Fig. 2C, left, two portions of color sequences are shown, corresponding to two different stimuli from Fig. 2B, for which two distinct classes of colors were expressed. These colors correspond to different classes of activity patterns (model vectors) that can be traced down in the Kohonen map (Fig. 2C, 3DKM and Pattern). We further computed pattern-triggered spike-raster histograms (PTSRHs), i.e., sum of spike-rasters before each occurrence of a given pattern (Fig. 2C, Spike Histogram). When the same pattern occurred multiple times consecutively, due to the decay of the integration kernel, only the first occurrence of the pattern was taken as a reference. In Fig. 2C, two PTSRHs are shown for a pattern with three and a pattern with two active neurons (neurons that define the pattern by firing a spike are referred to as “active”), corresponding to colors green and red, respectively. The time of pattern occurrence is taken as a reference at 0 ms in the PTSRH. The bright stripes in the bottom PTSRH shown in Fig. 2C indicate that active neurons of the pattern emit spikes frequently participating in bursts. That is, spikes in the pattern were often preceded by other spikes of the same neurons, at interspike intervals <8 ms (DeBusk et al. 1997). The model vectors (Fig. 2C, Pattern), represented as vectors with grayscale-coded entries (black = active, white = inactive), were very similar to the spiking patterns in the PTSRHs. Thus colors from the color sequences have corresponding model vectors in the 3DKM, which faithfully represent the actual firing patterns of the recorded neurons. These firing patterns can be traced down in the 3DKM and subjected to analyses also with quantitative methods (e.g., pattern occurrence rate, stimulus specificity, etc).
Color sequences become richer when a larger variety of stimuli is used in the experiment. In Fig. 3A, color sequences are shown for multiunit activity simultaneously recorded from a set of 32 electrodes under eight conditions that varied considerably in the spatial extent and complexity of the used stimuli (see methods, Plaid and bar stimuli). It is evident that similar stimuli elicit similar responses (respectively, stimuli 1 and 3, 2 and 4, 5 and 6), the spatiotemporal evolution of stimuli is reflected in the temporal expression of patterns (stimuli 2 and 4: the speed of the moving bar determines the temporal activation profile of the “green” and “purple” patterns), some stimuli produce clear onset (on) (stimuli 1, 3, and 5–7, at 1,050–1,060 ms) or offset (off) responses (stimuli 5–7).
Another advantage of the present analysis is that the succession of patterns in time can be conceived as a trajectory in a multidimensional phase space (Bathellier et al. 2008; Brown et al. 2005; Friedrich and Laurent 2001; Galán et al. 2004). Given that model vectors were mapped on the 3DKM lattice, these trajectories could also be visualized in a 3D space. The successions of patterns in Fig. 3A are shown as trajectories in 3D space in Fig. 3B. To this end, we first computed the average pattern within nonoverlapping, 50-ms-long sliding windows, across all the trials of a given stimulus. For each position of the sliding window, we marked the location of the average pattern in the 3DKM lattice and then connected successive points with lines. Consistently with the similarity of the color sequences in Fig. 3A, the 3D phase spaces in B were also similar for stimuli 1 and 3, 2 and 4, and 5 and 6, respectively. Thus one can also construct 3D phase-space representations of multineuron activity using the 3DKM. However, the 3D trajectory representation is less informative than color sequences and represents only an added value to the method, not its main purpose.
For more realistic stimuli, comprising natural scenes, neuronal responses can become more complex (Fiser et al. 2004). In Fig. 3C, color sequences are depicted for a dataset consisting of 34 SUs with responses evoked by three movies with natural scenes, each in length of 28 s (methods, Dataset with natural stimuli). The grouping by stimulus reveals preferential expression of particular patterns at different moments in time. In some cases, patterns were reliably stimulus-locked across different trials (e.g., Fig. 3C, stimulus 1, around 6,500 ms), whereas in other cases they were not. Moreover, the color sequences associated with natural movies contain sometimes fewer color blobs than for simple stimuli, and this varies across different movies. This property is expected given the more sparse responses of neurons to natural scenes (Vinje and Gallant 2000). Also an important factor determining the structure of color sequences could be the degree to which the elements of the movies stimulate the receptive fields of the investigated neurons. Finally, the color sequence representation permits a quick and comprehensive inspection of the entire dataset and selection of periods of interest for further quantitative analyses. One can, for example, subsequently determine whether some patterns appear preferentially for given stimuli, and evaluate their stimulus-locking properties. Color sequences for the entire duration of the movies are shown in Supplementary Fig. S3.
Multiple time scales and higher-order correlations
We next investigated the ability of the method to detect firing patterns that evolve on different time scales. To have a precise control over the time scale of patterns, we generated two artificial datasets containing patterns that evolved on either fast, i.e., synchrony, or slow time scales, i.e., firing rate (see methods, Artificial datasets). By manipulating the integration time constant (τ) when computing the activity vectors one biases the analysis to various time scales. In the color sequences from Fig. 4A, corresponding to the dataset containing synchronous events (joint-spike events: JSEs), JSEs could be clearly detected when a small time constant (τ = 5 ms) was used but could not be detected with a large time constant (τ = 50 ms). When a too large time constant is used, the activity vectors corresponding to JSEs are contaminated by the spikes preceding the JSEs. Conversely, for the dataset containing slow rate covariation patterns, using a small time constant hampers the detection of the specific firing rate combination because not enough information is integrated. Therefore activity vectors fail to properly represent the firing rate profile (Fig. 4B, top). A larger τ on the other hand permits the proper integration of the firing rate profile across different neurons and thus enables the representation of the specific rate pattern by the activity vectors and, for that matter, by color sequences (Fig. 4B, bottom). Thus by manipulating the integration time constant and plotting color sequences at different time scales, one can identify the expression of multineuronal firing patterns across multiple time scales. Moreover, one can also infer the characteristic time scale on which a given firing pattern occurs and this is expected to be especially useful in studies that investigate dynamical stimulus encoding over multiple time scales (Butts et al. 2007).
To illustrate the ability of the method to detect high-order spike correlations (e.g., JSEs) not extractable with pairwise methods, we generated an additional dataset with neurons firing homogenous Poisson spike trains at 10 Hz. We chose 20 trials of 30 s length and inserted in each trial a single JSE with five active neurons. Supplemental Fig. S4 shows that although cross-correlation histograms between all the pairs of active neurons in the JSE fail to reveal a clear central peak, the visualization method can successfully identify the JSE in the data. Importantly, given the color richness of the color sequences, one can easily detect the JSEs if they are time-locked to the same position in the trial. When the JSEs are not clearly stimulus-locked, the method probably provides little help in visually detecting them and more advanced quantitative techniques may be required. A possibility to solve this problem is presented in Identifying meaningful patterns.
Detection of nonstationary responses and cortical state-changes
Nonstationary responses can occur during the experiment due to unstable recording setup, cortical state changes (Arieli et al. 1996), or other factors (Gur and Snodderly 2006; Sannita 2006). Such nonstationary responses can be quickly identified using color sequences. An example is depicted in Fig. 5A, where color sequences corresponding to 12 SUs recorded in response to 14 stimuli (see methods, Dataset with center-surround stimuli) are shown in the order the trials were recorded (trials of different stimuli were presented in random order). Responses to the same stimuli were recorded from the same neurons in two different sessions, ∼26 h apart. Session 1 (Fig. 5A, left) exhibited marked nonstationary responses for successive trials along the recording. on and off responses were poorly defined. Periods with relatively stable responses were interrupted by periods of weak and unreliable responses. By contrast, in session 2 (Fig. 5A, right), cortical responses were more locked to input and stable over time, on and off responses were reliable and precise. We further quantified these observations by computing the coefficient of variation for the patterns (PCV; see methods), in a time-resolved fashion, and averaging it across stimuli. A large PCV in a time window denotes a large variability in the identity of patterns expressed across the trials corresponding to the stimulus at that moment in time. In other words, the larger the PCV the more unreliable the multidimensional response pattern is across repeated presentations of the same stimulus. Average PCVs were clearly larger for session 1 than session 2 (Fig. 5B), indicating a marked difference in the reliability of pattern expression over time.
When computing autocorrelation histograms for the responses in session 2, we found oscillatory activity of ∼27 Hz, which is in the beta-high band (20–30 Hz; Fig. 5C, top right) but did not find oscillations in session 1 (Fig. 5C, top left). We quantified oscillatory response properties of all neurons by computing their oscillation score (Mureşan et al. 2008) in the beta-high/gamma-low bands (20–40 Hz). Significantly higher oscillation scores were found for session 2 than for session 1 (P < 0.001; Mann-Whitney U test; Fig. 5C, bottom), confirming that neurons in session 2 exhibited robust oscillations, whereas the ones in the session 1 did not. This result shows that differences observed in the color sequences and computed quantitatively on patterns, were associated with a change in the oscillatory behavior of neurons for the investigated case. The result is also consistent with reports that firing rate responses are less variable across repeated trials in brain states with strong gamma oscillations (Rodriguez et al. 2004).
In another example, responses evoked by high spatial-frequency sinusoidal grating stimuli (methods, Dataset with high spatial-frequency sinusoidal grating stimuli), exhibited also marked nonstationarities during the experiment (Supplemental Fig. S5).
Identifying meaningful patterns
Color sequences reveal all patterns that evolve on a given time scale. However, not all patterns are evoked by the presentation of the stimulus and sometimes similar patterns are evoked by multiple stimuli. As a consequence, it is useful to determine how specifically given patterns occur for given stimuli. To this end, we developed a measure, PSI, that reflects the proportion of occurrences of a pattern for a given stimulus out of the total occurrences of the pattern throughout all stimuli (see methods). The PSI can be used to threshold color sequences: for a given stimulus, color bands are painted only if their corresponding patterns have a PSI larger than the threshold. In Fig. 6A, we show color sequences for a dataset recorded with sinusoidal grating stimuli and two thresholded color sequences with thresholds of 0.25 (B) and 0.5 (C). For higher thresholds, fewer but more specific patterns are revealed.
The thresholding may come in handy when specific patterns are not simulus-locked and therefore their discrimination is hindered by the color richness of color sequences. A high threshold eliminates most colors and keeps only those corresponding to patterns with high specificity. Therefore these become detectable regardless of their stimulus-locking properties.
The Kohonen self-organizing map provides an ordered clustering (Kohonen 2001) that creates a discrete approximation of the activity vector space. This approximation allows for a flexible representation of the occurring multidimensional activity patterns through a limited set of representative model patterns. As we have shown, by introducing a color-based representation of these patterns, one can visually inspect large datasets and draw conclusions on the temporal distribution of classes of firing patterns. Importantly, the ordered mapping provided by Kohonen maps enables the consistency of the color representation. Thus two different segments of color sequences sharing similar colors represent similar activity patterns of the recorded neurons. The visualization technique we have introduced is relatively straightforward to implement and conceptually simple.
The present method represents advancement over previous PCA-based visualization techniques (Bathellier et al. 2008; Friedrich and Laurent 2001; Galán et al. 2004). First, the Kohonen mapping is nonlinear and can map the structure of the input space onto the 3D lattice more flexibly while preserving local similarity relations (Lee and Verleysen 2007). This is not the case with linear techniques such as PCA, which can introduce large approximation errors if the structure of the input space cannot be linearly decomposed. Second, instead of emphasizing a state-space trajectory representation (Bathellier et al. 2008; Brown et al. 2005; Friedrich and Laurent 2001; Galán et al. 2004), albeit possible with our method as well, we introduced a representation based on color bands. This technique offers an important benefit to the visualization especially when nonsimilar activity patterns appear in a fast succession. In that case, a state-space trajectory will appear noisy and difficult to interpret. By contrast, the color bands in the color sequence can be inspected in detail, one by one, being perceptually associated through their color with the actual activation pattern of the underlying neurons. Moreover, the precise occurrence time of patterns within the trial is also straightforward to estimate in color sequences. Therefore the latter provide a more salient and flexible visual description of neuronal dynamics than trajectories in a reduced state space.
Unlike PCA, Kohonen mapping is not a dimensionality reduction technique per se but a topological mapping technique. In fact, the 3DKM is an ordered clustering (Kohonen 2001) and for that reason, it has the important advantage of ensuring both an approximation of activity vectors (clustering) and a topological mapping on the 3D lattice that enables the coloring technique we have described. The 3DKM preserves the original dimensionality of the data through model vectors (cluster centers), and therefore one has immediate access to classes of firing patterns that occur in the data and are segregated through clustering during the training of the Kohonen map. In the model vectors, the identity of individual neurons from the original recording is preserved. In addition, both silent and active neurons determine the identity of the pattern (model vector), i.e., silent neurons also matter. On the contrary, PCA reduces the dimensionality of the data and the univocal relation to the original, high-dimensional firing patterns is therefore lost, the contribution of individual units being difficult to recover.
Our method is related to the gravitational technique of Gerstein et al. (1985) in that it can flexibly manipulate the time scale to detect multidimensional activation patterns evolving over different time scales. However, in Gerstein et al., although computations are performed in a multidimensional space (multiple particles coalesce when neurons fire in a correlated fashion), the visualization provided therein can only represent pair wise distances between particles corresponding to neurons (Gerstein and Aertsen 1985; Gerstein et al. 1985; Lindsey et al. 1997) and therefore that method differs from the one introduced here with respect to visualization. Nevertheless some parallels may be drawn to the gravitational technique. For example, a high-dimensional correlated firing pattern could be represented by a color in our method and by a decrease in the pair-wise distances between the corresponding particles in the method of Gerstein et al. In general, the two methods are complementary and could be combined to explore multineuronal correlation structures in more detail.
Color sequences are especially useful when multiple recorded trials are available for the same stimulus. When the evoked activity patterns are stimulus-locked and color sequences are grouped according to their corresponding stimulus, the visualization technique can provide especially salient information. Similar colors would be grouped in larger portions of color sequences and will pop-out from the background. Relatively small temporal jitter of patterns across trials is also detectable, as we have shown in the examples with natural movie stimuli. Furthermore, one can estimate the reliability of pattern occurrence over multiple trials. When patterns are not locked to stimulus onset but exhibit a large range of time shifts relative to the stimulus, their detection is visually less salient, but nevertheless possible by careful inspection or by thresholding the color sequences according to the specificity of patterns. The same applies when only one single recorded trial is available and therefore there is a single color sequence.
Because activity vectors integrate information over a controllable time scale, color sequences also reflect the snapshot of the activity patterns in the data at a given time scale. The temporal dynamics of cortical networks and the characteristic time scale on which information is encoded has been thoroughly debated, especially for the visual system (Shadlen and Newsome 1998; Singer 1999). Color sequences can reveal specific activity patterns that evolve on various time scales, and hence, they may become handy in detecting the characteristic time scale of stimulus-specific multineuronal firing patterns (Butts et al. 2007).
The multidimensional framework relying on activity vectors is a generalization of population coding. In the limit, for very large integration time constants, one obtains firing rate vectors as patterns. These firing rate vectors can represent a population vector similar to that described by Georgopoulos et al. (1986). As the time constant is decreased, the emphasis is shifted toward synchronous rate covariations and finally to synchronous patterns with millisecond precision. In fact, the multidimensional vectors can be interpreted as correlated fluctuations of firing on the time scale corresponding to the integration time constant. Therefore this multidimensional representation of the spiking signals is general as has already been described in other studies that have combined exponentially decaying kernels with multidimensional analysis (Gerstein et al. 1985; Nikolić et al. 2007).
Activity patterns detected by the method reflect high-dimensional correlations between recorded neurons, and therefore as we have shown here, one can extract information that is not accessible to pair wise techniques such as the cross-correlation (Perkel et al. 1967). Activity vectors are integrated with realistic, exponentially decaying kernels that mimic synaptic currents (Häusler and Maass 2007; Mureşan and Savin 2007) and exhibit fading memory (exponential decay). Hence the activity vectors can also represent serial (temporal) correlation in addition to the high-dimensional spatial correlation (multiple neurons). This transcends the representational power of the frequently employed binning and binarization (Puchalla et al. 2005; Schneidman et al. 2006; Yu et al. 2008), which is inaccurate for large bin sizes and fails to consider serial (temporal) correlation (Roudi et al. 2009).
Importantly, the method described here does not distinguish between multineuronal activation patterns that arise either by chance or are due to some internal coordination process (Singer 1999). The purpose of the method is merely the representation of multiple spike trains in a format appropriate for visualization. Nevertheless, activity vectors can also be used in quantitative analyses and issues regarding their statistical properties, stimulus specificity, stimulus locking, etc, can be further investigated. Activity vectors represent a snapshot of the input currents that a hypothetical postsynaptic neuron would receive from the recorded neurons at a given moment in time (Nikolić et al. 2007). The visualization through colors helps identifying these patterns quickly, albeit with some degree of loss in precision.
One aspect that needs to be considered is that color perception is dependent on context (Purves and Lotto 2003). The distribution of colors in the sequences can either facilitate or hinder the discrimination of colors similar to each other. From this point of view, the visualization method is limited by the fact that given a certain context, two similar but distinct patterns may not be distinguished in the color sequences because their similar colors are visually indistinguishable. Therefore color sequences alone cannot provide accurate quantitative details and may need to be complemented by a further inspection of the model vectors in the 3DKM. Color sequences compress a large amount of multidimensional information onto a limited color space such that there is a tradeoff between the global visualization capability and the limited precision that the visual representation can offer. Nevertheless, the ordered mapping provided by the 3DKM ensures that confounded patterns are similar to each other and one can always go back to model vectors to make a precise quantitative estimate of the underlying neuronal activity patterns.
Kohonen maps have also been employed before for data classification (Nicolelis et al. 1999). However, data classification, although related to data visualization, has a notably different purpose. While in visualization the emphasis is on finding an intuitive representation for a human observer, in data classification, an artificial classifier is used to discover stimulus-specific structure in the recorded data. In both cases, the data may be transformed in nontrivial ways, such that it provides to humans or artificial classifiers features that are easy to learn and discriminate. Here we have shown how Kohonen mapping can be successfully used to enable such a transformation, facilitating visualization by a human observer.
It is important to mention that the visualization technique can be implemented using methods other than Kohonen maps. In principle, any topological mapping method that ensures topology preservation of the input space, combined with a way to assign colors based on topology, can substitute the 3DKM (e.g., multidimensional scaling, generative topographic mapping, locally linear embedding, etc). We chose Kohonen maps because they are able to preserve the topology of the input space (activity vector space), cope well with large numbers of input samples (no need to compute distance matrices), and their predefined topology in the low-dimensional space can be mapped directly onto a 3D color space. The nonlinearity of the Kohonen map is another advantage. For example, between two patterns located at the extremes of the 3DKM cube one does not find a smooth transition from the first pattern to the second but can recover notably different patterns. Only the local similarity is ensured in the Kohonen map and this enables an efficient, less redundant use of the 3D lattice for representing the stereotypical firing patterns.
High-dimensional datasets are difficult to analyze and especially difficult to understand and explore graphically. The visualization technique we have introduced here can provide graphical/visual information about the activation patterns of all neurons across all trials of a given experimental session and can help guide subsequent quantitative analyses. In the future, this visualization technique may be applied to other types of multidimensional signals such as electroencephalography, magnetoencephalography, or local field potential data.
This work was supported by two Romanian grants (MECT/UEFISCSU: RP-5/2007 Contract 1/01.10.2007, ID-48/2007 Contract 204/01.10.2007), a Max Planck–Coneural Partner Group grant, Deutsche Forschungsgemeinschaft Grant NI 708/2-1, European Commission Grant GABA-FP6-NEST/043309, and the Hertie Foundation.
We thank the Science Museum of Minnesota for providing free of charge use of “The Greatest Places” movie and also V. V. Moca, R. V. Florian, S. P. Paşca, C. Wolff, and N.-H. Chen for useful comments and assistance. Free source code and libraries are available at: http://www.ovidiu.jurjut.ro/sources/visualization/.
↵1 The online version of this article contains supplemental data.
- Copyright © 2009 the American Physiological Society