One of the key innovations of the SIMNETS framework consists of using single neuron SSIM matrices as a “computational fingerprint” to represent the intrinsic structure of the output space of individual neurons. This is accomplished by performing within-neuron measures of pairwise spike train similarity across all spike trains for each neuron, resulting in an SxS single neuron SSIM matrix for each of the N neurons (Figure 1 step 2; Figure 2B). Multiple approaches to estimate the similarity between pairs of spike trains have been proposed (Paiva et al., 2010) and could be used for this step. Here, we use the spike train metric proposed by Victor and Purpura (Houghton and Victor, 2010; Victor and Purpura, 1997; Victor and Purpura, 1996; See Box 1 for description).

In the current work, each single neuron SSIM matrix represents the relationship (i.e., similarity) between all spike trains for a single neuron for the specific time windows selected for all simultaneously recorded neurons. Figure 2. presents three examples of SSIM matrices derived from the spike trains of simulated neuron responses (Linear Non-Linear Poisson functional model; Meyer et al., 2016). Entries equal to zero (black) in the matrices represent identical spike trains (zero VP cost), while higher values represent increasing differences between the spike trains. For ease of viewing, the rows/columns of the matrices have been sorted according to simulated experimental conditions. Note that a specific ordering of the spike trains is not required for the subsequent comparisons (and is only included for illustrative purposes). It can be useful to graphically illustrate the relationship between individual spike trains across the experiment captured by the single neuron SSIM matrices using dimensionality reduction (Figure 2B). Here, we use PCA-initialized tSNE (Maaten and Hinton, 2008): a combined unsupervised dimensionality reduction technique that aims to preserve the local neighborhood structure, via t-SNE, as well as information related to the global shape, via PCA-initialization (Kobak and Linderman, 2021; Lee et al., 2015). The non-random PCA-initialization also ensures reproducibility across iterations (Kobak and Linderman, 2021). This mapping step is not necessary for the following steps of the SIMNETS analysis pipeline, but these maps can be useful to visualize and compare the features coded in the population. For example, in Figure 2B, the single neuron SSIM maps illustrate the separation of orientation selectivity in simulated neuron n1, the spatial frequency selectivity in n2, and the combined spatial and orientation selectivity of n3 (See Supplementary Figure S1; supplementary material for additional discussion and examples). Changes in the shape of the SSIM output space map may also indicate the influence of other variables which are not explicitly manipulated.

Each SSIM matrix characterizes the intrinsic geometry of the output space for a specific neuron. The next step of the algorithm compares these “computational fingerprints” across the full set of recorded neurons. This second-order comparison is made by correlating each pair of SSIM matrices (Figure 1, step 3; Figure 3). Here we use Pearson’s correlation (Figure 3A). Each correlation value is entered into a symmetrical (N x N) matrix (Figure 3B) representing the relationships across neurons, which we term the Computational Similarity (CS) matrix. Applying dimensionality reduction (here, PCA-initialized t-SNE; Maaten and Hinton, 2008) to the CS matrix creates a CS map (Figure 1, step 4; Figure 3C). Note that, while each colored point in a single neuron SSIM map represents an individual spike train (Figure 2C), each point in the CS map represents a neuron (Figure 3C). In the population CS maps, the distance between neurons is a relative measure of the similarity of computational fingerprints (the SSIM matrix). Therefore, neurons performing similar computations will tend to aggregate in clusters within the population CS map of neurons.

After neurons are embedded in the CS map, the task of finding groups of neurons representing potential neuronal subnetworks can be addressed using a wide variety of clustering algorithms. Here, we use the k-means clustering algorithm, an efficient centroid-based clustering method, however, other density-based, distribution-based, or hierarchical-clustering methods would also be appropriate. The k-means algorithm works by iteratively assigning neurons into the pre-specified number of k clusters until the optimal solution is reached. However, the challenge is to ensure that the clustering structures identified by the k-means algorithm reflect genuine clusters within the CS neuron map and not false discoveries. To address this issue, we combined the k-means algorithm with two different cluster validation techniques, a silhouette graphical analysis and our novel shuffle-based statistical test, which together help avoid false cluster discovery (Figure 4).

Silhouette analysis is a validation technique used to determine the most likely number of clusters in a dataset (Figure 4; Rousseeuw, 1987). A silhouette value represents the ratio of the between to within-cluster distances: a point within an ideal cluster will be close to members of the same cluster and far from points assigned to different clusters, resulting in a high silhouette value. Choosing the cluster number with the highest average silhouette value—the Silhouette Coefficient (SC) score—maximizes the separation between potential clusters (Figure 4C, red point). We use a novel shuffled-based resampling procedure to assess the statistical significance of the SC score to avoid false cluster discovery. The test relies on creating 1000’s of surrogate datasets through a shuffling or randomization procedure that is applied to the rows/columns of the single neuron SSIM matrices We use the matrix permutation technique used in the Mantel test, which was specifically designed to randomize pairwise distance matrices. This is accomplished by permuting both rows and columns together in order to preserve matrix symmetry. In this way, it is possible to produce a new distance matrix for each neuron where the relationship between the trials (i.e., which trials are most similar) is randomized. This procedure effectively destroys any dependencies among the neuron’s single neuron SSIM matrices and, consequently, the associated clusters with the CS neuron map (Figure 4B). The empirical SC score is only considered statistically meaningful if it falls outside a distribution of surrogate SC values calculated from the 1000’s of surrogate datasets created in this manner (Figure 4C, see blue dots; for further insights on conditional resampling; see Amarasingham et al., 2012; Schemper, 1984). This unsupervised cluster detection and validation procedure supports the identification of statistically meaningful clusters of computationally similar neurons within the CS neuron maps.

We evaluated the performance of the SIMNETS algorithm using simulated neural data with known ground truth. This allowed us to verify that the algorithm clusters neurons with similar informational content even when they use different encoding schemes. We simulated a population of 180 neurons (N) composed of three computationally distinct subnetworks, SN1, SN2, and SN3 (Figure 5A). We generated an output structure for these subnetworks across three simulated experimental conditions (C1, C2, and C3). Each subnetwork emitted a similar “baseline pattern” for two of the conditions and a different pattern for a unique “preferred” condition (C1 for SN1, C2 for SN2, and C3 for SN3). Neurons from a common subnetwork were designed to encode similar information (i.e., respond differentially to one specific condition) using different encoding formats; for example, 20 neurons in SN1 encoded condition A through a change in spike rates (rate-coding), another set of 20 neurons in the same subnetwork responded with a change in the timing of their spikes (temporal-coding neurons), while the remaining neurons responded with changes in both spike rates and the spike time variations (mixed coding; e.g., Figure 5B). The other two subnetworks were also generated to include examples of rate coding, temporal coding, and mixed coding artificial neurons in the same proportions. We simulated 30, one-second-long spike trains for each of the 180 artificial neurons, which included 10 repetitions of each condition (S = 30 spike trains per neuron; See Supplementary Table S1 for analysis summary).

SIMNETS can be used to quantify spike train similarity and explore single unit encoding properties across an experiment without assuming what coding schemes are relevant or what variables are generating spike train differences. It is also possible to examine the temporal accuracy of single unit spiking. To illustrate this point, the SIMNETS algorithm was applied using three different temporal accuracy values for the VP spike train similarity metric: a pure rate code (q = 0, i.e., no added cost for shifting spikes in time), 100 ms (q = 10), 5 ms (q = 200). The q parameter determines how sensitive the similarity metric is to differences in the timing of individual spikes, by setting an upper limit on how far spikes can be shifted in time (See Figure 5C, for three different single neuron SSIM matrices calculated using different q values). With a setting of q = 0, the simulated neurons operating with a rate-based encoding scheme and mixed coding neurons are grouped into three functionally distinct clusters in the CS map (Figure 5D, 2nd column), while temporal coding neurons form a single cluster at the center of the map (Figure 5D, 3^rdcolumn). As the value of q increases, the VP algorithm “cost” function becomes sensitive to differences in spike timing in addition to the total number of spikes (Victor and Purpura, 1996). SIMNETS correctly groups all neurons into three distinct clusters that reflect the ground-truth functional ensemble assignments (Figure 5F). At the highest q value (Figure 5H), the CS map shows sub-groupings within each of the three ground-truth subnetworks that correspond to neuron’s encoding formats (Figure 5H, 3^rd column); however, the optimal number of clusters remains in agreement with the ground-truth functional subnetwork assignments (K_sc = 3). By specifying a higher partition value for the k-means clustering step of the algorithm (e.g., k = 9), the sub-groupings within the detected clusters are defined by the coding properties of the neurons. For a demonstration of the interaction between the cluster number and the SIMNETS hyperparameters, perplexity and q, see Supplementary Figure S2.

To evaluate whether SIMNETS could successfully identify subnetworks, we compared the distribution of Euclidean distance between individual neurons in the SIMNETS CS map for neuron pairs within and between the artificially generated (‘ground truth’) subnetworks (Figures 5E,G,I, histogram). The within-subnetwork similarity was significantly higher than between-subnetwork values in all cases (Mann–Whitney p < 0.001). For q values > 0, there was no overlap between the two distributions, indicating the complete separation of the artificially generated subnetworks. Our results demonstrate that the SIMNETS algorithm can accurately separate neurons according to their computational properties (the latent variables upon which they operate), even if they employ different coding schemes to represent information using spike rates or temporal patterns.

We next applied SIMNETS to three distinct datasets of simultaneously-recorded, multiple single-unit recordings from three different areas in two species: (1) nonhuman primate (NHP) primary visual cortex (V1) during visual stimulation (Kohn and Smith, 2016; Smith and Kohn, 2008), (2) NHP motor cortex (M1) neural recordings during a center-out reaching task (Rao and Donoghue, 2014), and (3) rat CA1 hippocampal region recordings during a left–right alternation maze task (Pastalkova et al., 2015; Pastalkova et al., 2008). These datasets were chosen because the functional roles of the neurons in these areas have been reasonably well-established and the information content of the single neuron and ensemble spiking patterns extensively described (Gur and Snodderly, 2007; Livingstone and Hubel, 1984; Rao and Donoghue, 2014; Smith and Kohn, 2008; Truccolo et al., 2008).

We note that these results are meant to illustrate the capabilities and utilities of the SIMNETS analysis pipeline and are not intended to represent a full characterization of the role of any identified putative subnetworks. For a full analysis and interpretation, many more datasets (sessions/subjects) would be required, as would various parameter sweeps to evaluate computational processes.

We analyzed a dataset of 112 simultaneously recorded V1 neurons (recorded using a 96-channel “Utah” electrode array) during the presentation of drifting sinusoidal gratings in an anesthetized Macaque monkey (Kohn and Smith, 2016; Figure 6A). We extracted 1 s of spiking data from the first 30 repetitions of each stimulus (S = 360), starting 0.28 s after stimulus onset (see Figure 6B for raster plots and single neuron SSIM matrices for two example neurons) to capture the neural stimulus response period (eliminating response delay). SIMNETS was applied to these NxS spike trains generated with a fixed temporal accuracy setting of 1/q = 50 ms. The resulting NxN CS neuron matrix is projected into a lower-dimensional CS space upon which further statistical and clustering analyses are performed (Figures 6D–F). An unsupervised k-means cluster analysis identified 3 clusters of related single neurons in the V1 CS map (Figures 6D,G; SC = 0.78, peak max average silhouette value), indicating that computational features did not fall along a functional continuum but instead were organized into three functionally distinct putative subnetworks.

To validate the performance of the SIMNETS algorithm, we examined the relationship between the structure of the CS map and the known tuning properties of V1 neurons. We used classic parametric tuning models to estimate the receptive field properties of all neurons and labeled the neurons in the CS neuron map according to their “preferred” stimulus orientation (Figure 6E) and tuning depth (Figure 6F). Each neuron’s preferred orientation was calculated by fitting a Gaussian distribution to the stimulus-dependent firing rates and finding the orientation that maximizes the function over the range of orientations angles (θ = [0,180]) across both left and right drift directions. Tuning depth was calculated as the normalized difference between the peak and trough of the tuning function (See Supplementary methods). Again, we emphasize that the SIMNETS procedure is an entirely unsupervised approach and does not rely on any knowledge of the experimentally-induced neuron tuning functions to generate the CS map or to identify putative CS neuron clusters. We simply labeled each neuron with values obtained from its relationship to an experimental variable (here, grating direction) after the neurons were plotted in the SIMNETS CS map. A pairwise circular-linear correlation (r_cl) statistical analysis revealed a strong positive correlation between neurons’ distance in the low-dimensional CS neuron map and the differences in their preferred orientation parameters (Pearson, r_cl = 0.89). This relationship confirmed that neurons with similar preferred orientation were generally mapped near one another in the CS map (Figure 6E). The exception to this trend was observed for neurons with weaker tuning strengths (Figure 6F, brown points), which aggregated at the center of the CS neuron map. Overall, the distribution of the neurons’ preferred orientations corresponded well with the three observed clusters in the CS neuron space (Figure 6H). Collectively, these results confirm that the neurons’ sensitivity to the different stimulus orientations is a major contributor to the structure of V1 computational architecture. The CS map also reveals the presence of three computationally distinct groups of neurons within the context of this experimental paradigm.

The overall agreement between known V1 tuning features and the CS map generated using SIMNETS suggests that the algorithm will be able to identify computational architectures for neural systems where the tuning functions are not known a priori.

The DR techniques used to visualize the computational geometry representations of the single neuron SSIM maps (introduced in Figure 2) can be expanded to represent the combined output space of multiple neurons simultaneously, such as the CS neuron clusters (i.e., putative subnetworks) or the full neuron population. For multi-neuron SSIM maps, each point corresponds to the collective spiking pattern generated by all neurons within an ensemble (e.g., CS neuron cluster or full population) on a single trial, and the distance between points corresponds to the difference (spike train distance) in the neurons collective spike pattern across time. Here, we use this approach to visualize the relationship between activity patterns for the full neuron ensemble (Figure 7B), the three CS neuron clusters (Figure 7C), and example single neurons from each of the three CS neuron clusters (Figure 7D). Standard single neuron plots (parametric tuning functions and rasters) were also included for comparison.

As a second example, we applied the SIMNETS algorithm to a dataset of 103 M1 neurons recorded using a 96-channel electrode array in a macaque performing a planar, 8-direction arm reaching task (see Methods for more details). We extracted 1-s long spike trains (S = 114) from each neuron during all trials where the monkey successfully reached the cued target, starting 0.1 s before movement onset (Figure 8A; see Supplementary Table S1 for algorithm inputs and parameters). SIMNETS was applied to these NxS spike trains generated with a fixed temporal accuracy setting of 1/q = 50 ms. SIMNETS analysis revealed that the neurons were organized into three neuron clusters in the M1 CS map (Figures 8D,G, SC score = 0.71, M1 peak average silhouette).

To examine the relationship between these clusters and functional properties of the neurons, we characterized each neuron’s preferred reach direction (PD) angle and tuning depth (Figure 8F) in a similar manner to the V1 dataset. Each neuron’s PD was estimated by fitting a von Mises distribution to the firing rates as a function of direction (Supplementary Figure S3; Mardia and Zemroch, 1975) A circular-linear correlation (rcl) analysis between the difference in PD and CS map distance between pairs of neurons revealed a significant positive relationship (Pearson, rcl = 0.92; p = 0.001), confirming that neurons with similar preferred directions tended to cluster together in the CS map, while neurons with dissimilar preferred directions were spatially distant (Figures 8D,E, right plot). After plotting the distribution of preferred reach direction, we again observed that the neuron tuning properties were similar across each of the three detected CS neuron clusters (Figure 8H).

SIMNETS analysis suggests that the computational properties of M1 neurons are not uniformly distributed along a functional continuum, as evidenced by the distinct CS neuron clusters identified in the population CS space. The separation between the CS neuron clusters is most apparent for CS neuron cluster 2 and cluster 3 (Figure 8D, orange and yellow CS clusters), which mirrors the discontinuity in the distribution of preferred reach angles (Figure 8E). These results agree with previous findings, which support the hypothesis that the biomechanical constraints of the limb are reflected in an uneven distribution of preferred directions among motor cortical neurons (Lillicrap and Scott, 2013).

As with the V1 dataset, we used a multi-scale latent space analysis to visualize the relationship between spiking patterns across different scales within the CS map (Figure 9A), from the full ensemble (Figure 9B) to CS neuron clusters indicative of potential subnetworks (Figure 9C) and individual neurons functional properties (Figure 9D). Standard single neuron plots (parametric tuning functions and raster plots) are again included for comparison (Figure 9C). Each M1 neuron had a unique single neuron SSIM map that captured several computationally relevant features of the neurons’ spike train outputs (Figure 9C; SSIM matrices/maps for two example neurons). In contrast to the regular conical and circular parabolic shapes of the V1 single neuron SSIM spaces (Figure 7C), the M1 SSIM maps displayed more irregular and heterogenous global geometries (Figure 9D).

Lastly, we assessed the CS structure of a population of hippocampal neurons using SIMNETS, a region that is both functionally and structurally quite different from the primary visual and motor neocortical areas examined above. We applied SIMNETS to a publicly available dataset of rat CA1 hippocampal neurons recorded using multi-site silicon probes in a rat performing a left/right-alternation navigation task in a “figure-8” maze (Pastalkova et al., 2015; Figure 10A). All neurons with firing rates greater than 5 Hz were included in the analysis (N = 80/106 recorded neurons). The rat performed 17 correct trials (T = 17, trials), taking an average of 4.3 s to reach the reward location at the end of the left or right arm of the maze. Spike trains used for the SIMNETS algorithm were obtained by extracting 1-s of spiking activity for each 80 mm increment of distance traveled by the rat (Figure 10B, left, example single neuron raster plots with S spike trains; see Supplementary methods). This resulted in a set of S = 442 spike trains (17 trials x 26 time-epochs) per neuron.

SIMNETS was applied to the set of SxN neurons, using the VP temporal sensitivity setting of q = 35 (Figures 10B,C; see Supplementary Table S1 for all algorithm inputs/outputs). After generating the low-dimensional CA1 population CS neuron map, we characterized the organizing properties of the neuronal relationships using the unsupervised k-mean cluster detection/validation procedure. This procedure identified Ksc = 5 statistically significant CS neuron clusters within the CA1 CS neuron map (SC = 0.70, peak max average silhouette value; Figure 10D), suggesting that the neurons are organized into five computationally distinct clusters of neurons. To verify this, we used two different classic functional characterization approaches to assign functional ids to the neurons within the SIMNETS CS map and assess their functional organization within the context of the CS neuron map (Supplementary Figure 10F). Specifically, we categorized the CA1 neurons as having place cell-like activity (n = 60, place-cells, PCs) or lacking place cell-like activity (n = 20, non-place-cells, non-PCs) using an information-theoretic measure (Skaggs et al., 1992) and a spatial firing rate tuning functions (Figure 10B, example single neuron raster plots; see Supplementary methods for analysis details). As is reported in other studies, some of the neurons categorized as place-cells exhibited complex firing responses, such as multiple peak place-dependent firing responses (Supplementary Figure S3D, show example single neuron response functions). An examination of the physiological properties of the identified k-means neuron clusters revealed that one cluster was almost entirely (95%) composed of non-PCs (Figures 10E,F, K-means cluster c3, purple), while the other four k-means CS neuron clusters were entirely or almost entirely (>92%) composed of PC’s (Figures 10E,F). This indicated one of at least three possibilities: (1) the “non-PC CS cluster,” cluster c3 corresponded to a cluster of untuned, “uninformative” neurons, (2) that the individual neurons computed a common task-related signal that was unrelated to place-coding, or (3) that place-coding signals were only apparent at the level of the neuronal cluster rather than the individual neurons (Figure 10E, cluster c3). While a full exploration of these different possibilities is beyond the scope of this work, we do provide a demonstration of how the SIMNETS CS neuron maps may be used to visualize the relationship between spiking patterns across scales for hypothesis formation and testing (Figure 11).

Once again, we used a multi-scale latent space analysis to visualize the relationship between spiking patterns across different scales within the CS map (Figure 11A), from the full ensemble (Figure 11B) to CS neuron clusters indicative of potential subnetworks (Figure 11C) and individual neurons functional properties (Figure 11D). As expected, the spatial layout of the task-space was reflected in the topology of the population-level SSIMS map (Figure 11B), while each cluster-level SSIMS maps captured the place-dependent modulation of the clusters spiking patterns as the rat traversed a specific region of the maze (Figure 11C). Surprisingly, the “non-PC neuron cluster,” (Figure 11C, cluster c3) did not appear to reflect a cluster of non-responsive or “uninformative” neurons; instead, it had a cylindrical shape that resulted from at least two dominant patterns of spike pattern modulation: a within-trial periodic spiking pattern in the x-y plane (Figure 11B, cluster c3, x-y viewing angle) and a location or distance-dependent variance along the z-axis (See Supplementary Figure S3C, for x-z viewing angle). This cylindrical shape was also observed in the single neuron SSIM maps of the neurons within cluster c3 (see Figure 11C, purple box/stars, neurons c3-n7 and c3-n20), however, the distance-dependent x-z variance was less apparent at the single neuron level. Overall, this suggests that the activity of the non-PC neuron cluster, c3, displays dynamics that may reflect a different task variable that is in some way correlated with the rat’s performance on the task (e.g., head direction, speed profile along maze arms) or, alternatively, the intrinsic circuit dynamics (Pastalkova et al., 2008; Richard et al., 2013; Wang et al., 2015). Latent spaces displaying 3-dimensional ring topology reflecting animal position and running direction encoding have been recently described in CA1 (Esparza et al., 2025). Further investigation of this phenomenon is outside of the scope of the present work but is an example of how SIMNETS may be used to quickly uncover different types of functional relationships among neural representations that might not be otherwise evident when using classic population and single unit analyses. These results highlight the advantages of applying SIMNETS to neural recordings where tuning properties or other coded variables are not readily apparent or unknown a priori.

The SIMNETS framework can rapidly analyze large numbers of neurons using a standard personal computer (e.g., 6-Core Intel Core i7, 64GB of RAM). With the current Matlab® implementation (see included SIMNETS software package), the simulated neuron population (N = 180, S = 30) and Macaque M1 neuron population (N = 103, S = 360) were processed in less than 2 s. Larger datasets, such as the Macaque V1 neuron population (N = 112; S = 360) and rat hippocampal CA1 neuron population (N = 80, S = 442), take approximately 10 s to process. This is significantly faster than a typical cross-correlation-based functional connectivity analysis. For example, a dataset of 1,000 neurons can take an hour to processes with a standard cross-correlation analysis (Figure 12, blue line; utilzing Matlab® xcorr function), whereas our implementation of the SIMNETS algorithm takes only about a minute to process the same dataset (Figure 12, black line). The SIMNETS analysis is faster because it involves fewer spike train comparison operations (NxSxS) than the cross-correlations analysis (N2xSxS, across all time-lags), and the number of spike train comparisons per added neuron has a favorable linear scaling. If parallel or distributed computing resources are used to calculate the single neuron SSIM matrices, the execution time can be reduced even further to seconds or milliseconds. This would be especially useful for datasets with larger trial numbers (e.g., S > 1,000), which are more computationally intensive than datasets with larger neuron numbers.

We tested and validated SIMNETS using simulated data with known ground truth as well as three different experimental neural recording test datasets to demonstrate the general applicability of the analysis framework across a range of coding schemes, neural systems (visual, motor, cognitive), behaviors (awake or anesthetized), and different species (rat and macaque). This included a demonstration of a novel shuffle-based statistical procedure based on the Mantel test for identifying non-spurious CS neuron clusters within a low-dimensional population CS neuron map.

Our analysis of simulated data with ground-truth subnetworks and known single neuron properties (i.e., encoded content and encoding timescale) demonstrates the current implementation of the SIMNETS algorithm can capture a neuron’s local and global computational relationships, even under conditions where computationally similar neurons relied on very different encoding schemes (e.g., rate coding vs. fine temporal coding schemes). This agnosticism can be beneficial when analyzing neurons that generate highly heterogeneous spike train outputs because of differences in their internal machinery (e.g., inhibitory vs. excitatory neurons) yet still reflect task-relevant information (Churchland and Shenoy, 2007; Diba et al., 2014; Martina et al., 1998). We demonstrated how the VP metric temporal accuracy parameter (q) could be used to interrogate the effect of temporal resolution on the computational configuration of the neurons.

The application of SIMNETS to three publicly available multi-neuron recording datasets validated the capabilities of the method in revealing known dominant feature representation in V1, M1, and CA1 (orientation, direction, place, respectively) without imposing stimulus or movement-driven parametric tuning models a priori. Of note, the results show that the computational properties of neurons are not uniformly distributed, but instead are organized into apparent subnetworks (i.e., CS clusters in the low-dimensional CS map) that emphasize certain collections of feature representations. Our results also suggest that SIMNETS may be able to detect distinct subnetwork structures hypothesized to support ensemble place-coding or complex feature conjunctions (Buzsáki, 2019; Diba et al., 2014). Although it was beyond the scope of this report to demonstrate the functional or computational significance of the detected putative subnetworks, our results suggest that the detected clusters are physiologically meaningful; our subsequent unsupervised cluster validation statistical procedure further provided support to these findings. Finally, our choice of datasets allowed us to demonstrate that this method generalizes well to neural recordings from a variety of brain regions, including sensory areas (Kohn and Smith, 2016), motor areas (Rao and Donoghue, 2014), and memory/cognitive areas (Pastalkova et al., 2008); different species, including rat and non-human primate; and across different recording technologies including laminar silicone probes (Buzsáki, 2004) and fixed-dept multi-electrode arrays (Maynard et al., 1997).

Collectively, our results demonstrate that single neuron SSIM matrices provide a simple, scalable, and powerful format for characterizing the computational properties of neurons, as well as identifying computational similarities and differences between them. In addition, examining the arrangement of neurons in the CS maps may suggest novel coding schemes that are yet unexplored by experimenters, generating new hypotheses for future experiments.

The SIMNETS analysis framework builds upon rich theoretical and mathematical literature spanning multiple domains and disciplines. Geometric and metric space representational models of similarity have a long history of application in the field of psychology where they have been used to model the perceptual relationships between sensory stimuli as latent or low-dimensional perceptual metric-space (Forstmann, 2011; Goldstone et al., 1991; Hyman, 1974; Shepard and Chipman, 1970; Shepard, 1987; Shepard, 1964; Thurstone, 2017; Tversky, 1977; Zenker and Gärdenfors, 2015). This approach has recently been adapted to study the intrinsic structure of neural representational spaces in visual, motor, and cognitive brain regions (for reviews, see; Edelman et al., 1998; Haxby, 2012; Kriegeskorte et al., 2008; Kriegeskorte and Kievit, 2013; Lehky et al., 2013). Several mathematical variations of this framework have been developed to study the information content of macroscale fMRI BOLD signals (Diedrichsen and Kriegeskorte, 2017; Haxby et al., 2001; Kriegeskorte and Bandettini, 2007), and both population-level and neuron-level spiking patterns (Houghton and Victor, 2010; Kiani et al., 2007; Vargas-Irwin et al., 2015a, 2015b; Victor and Purpura, 1996). To the best of our knowledge, this is the first work that uses a measure of 2nd-order point-process similarity analysis to generate latent space embedding of a collection of individual neurons.

The concept of a low-dimensional embedding that captures the functional relationship between spiking neurons was introduced in the seminal papers by (Baker and Gerstein, 2000; Gerstein and Michalski, 1981) describing the use of “Gravitational Clustering” (Wright, 1977): a neuron clustering and visualization tool for identifying groups of neurons with synchronous spiking patterns. This method is based on an analogy of the physics of the gravitational forces governing the dynamics and interactions of macroscopic particles. It treats the N neurons as N particles moving within an N-dimensional space, where charges that influence the attractive and repulsive interactions between particles are dictated by the temporal dynamics of pairwise synchronous spiking activity between neurons. The result is a visualization of particle clusters (and their trajectories) that represent dynamically evolving assemblies of synchronously active neurons. More recent work by Kiani et al. (2015) involves the application of the dimensionality reduction techniques (i.e., multidimensional scaling) to pairwise measures of between-neuron spike rate covariations to detect natural functional modules in a population of pre-arcuate and motor cortex neurons (Kiani et al., 2015). Although the goal of this method is like that of SIMNETS—to the extent that it makes use of single-trial information to group neurons in an unsupervised manner—this approach is applied to binned spike rates and appears to capture the covariation information carried in the absolute firing rates of the neurons’ responses across trials (Kiani et al., 2015), rather than the information carried in the neurons’ intrinsic spike train relational geometry. However, our analysis of the synthetic neurons using temporal codes, as well as rat CA1 neurons demonstrate that SIMNETS can detect a broader range of meaningful information processing motifs that reflect both condition-independent and condition-dependent information carried in the fine temporal structure of spike train outputs.

Yang et al. (2019) identified functional groupings of artificial neurons in a neural network by applying t-SNE to activation-based measures of single neuron selectivity across 20 different tasks (i.e., task variance; Yang et al., 2019). The method is similar to the SIMNETS method in that it projects neurons onto a CS map according to the similarity of measures of output variance. This method requires knowledge of trial labels within tasks to group neurons according to the similarities of their changes in selectivity across tasks. By contrast, SIMNETS does not require a priori knowledge about the similarities or differences of the information encoded on certain trials (i.e., trial labels), only that the trials are recorded simultaneously. This feature of SIMNETS is particularly valuable when trying to identify groups of functionally similar neurons in experiments with awake and freely moving animals, where the assumption of repeatable perceptual, cognitive, or behavioral states across trials is not always possible.

Several previous studies have used pairwise comparisons between the spike trains of different neurons to identify putative subnetworks or cell assemblies (Aertsen et al., 1989; De Blasi et al., 2019; Gerstein et al., 1985; Grün et al., 2002a; Humphries, 2011; Roudi et al., 2015; Singer and Gray, 1995). As with the gravitational clustering method, these studies have operated under the working hypothesis that the detection of time-series spike time or rate correlations is a signature of a potential functional link (Cole et al., 2016; Roudi et al., 2015; Singer and Gray, 1995; Ts’o et al., 1986). By contrast, SIMNETS identifies neurons with similar informational content even if they exhibit heterogeneous spike train outputs or utilize different encoding timescales (e.g., rate vs. precise spike timing). For the simulated neuron population, SIMNETS successfully identified clusters of the simulated neurons according to their ground-truth functional subnetworks and by increasing the sensitivity of the VP metric, the CS map was able to reveal further sub-groupings within the three identified functional CS clusters that highlighted subtle differences their output space geometries that corresponded to the different encoding timescale. This feature of SIMNETS could be particularly useful for determining if neurons encoded information across heterogenous timescales, e.g., fast vs. slow-timescale inhibitory neurons, or if they compute across a hierarchy of timescales (Cusinato et al., 2022; Fortenberry et al., 2012; Gorchetchnikov and Grossberg, 2007; Honey et al., 2012). Additionally, this feature of SIMNETS enabled us to cluster physiologically and functionally distinct groups of neurons in the CA1 dataset (i.e., place-cells vs. non-place cells), which were not obvious when VP metric was insensitive to spike times (i.e., VP metric q = 0; data not shown).

Capturing the computational similarity between pairs of neurons in a more general way could be accomplished using Information theoretic approaches (Shannon, 1997), which could include estimating the shared or mutual information between the output spike trains of different neurons. Other approaches focus on the asymmetry of the predictive power between variables at different lags, resulting in “directed” estimates of functional connectivity such as Granger Causality, Transfer Entropy, or the Directed Transfer Function (Granger, 1983; Ito et al., 2011; Ursino et al., 2020). These strategies are based on estimating joint probability distributions across the activity patterns of pairs of neurons. A similar approach can be applied to relationships between multiple neurons using Generalized Linear Models (Chen et al., 2009; Truccolo et al., 2005). The number of possible activity patterns is very large, so this type of calculation can be challenging even when relatively large amounts of data are available. Using large amounts of data presents additional problems for this approach because it assumes that relationships between neurons remain relatively constant across the time span used to fit the models. Hence, every additional neuron requires exponentially more data and computing power and analyzing hundreds to thousands of simultaneously recorded neurons becomes intractable.

A considerable amount of research in systems neuroscience has focused on identifying new classes of neurons based on their information-processing properties. The standard approach for many of these experiments involves recording single unit activity while a certain experimental variable of interest is manipulated (for example, providing systematic changes in stimulus features, or eliciting different behavioral responses; Meyer et al., 2016). Standard statistical tests (ANOVA, etc.) are then used to determine if each neuron displays significant changes in firing rate across the experimental conditions. The percentage of significant neurons is usually reported as a functionally distinct “class” of neurons sensitive to the variable of interest. It is common to exclude neurons that do not reach statistical significance or cannot be fit using a predetermined tuning model from further analysis. This approach is prone to both selection and confirmation bias, and ultimately produces “classes” of neurons identified based on arbitrary statistical thresholds imposed on what are likely continuous distributions of properties. The SIMNETS analysis framework is an unsupervised approach to determine if neurons are organized across a functional continuum or are organized into statistically separable functional classes, thereby mitigating the experimenter’s bias inherent in parametric neural discovery methods. The SIMNETS approach is different in that it seeks to organize the neurons into a low dimensional computational similarity (CS) map solely based on the intrinsic structure of their firing rates, without making any a priori assumptions about their tuning properties. This makes it possible to detect clusters of neurons with similar informational content and can also reveal gradients of functional properties if they exist. This is not to say that experimental design does not play an important role: the informational content of neuronal spiking is still evaluated within the specific context of the recording. However, SIMNETS does not target specific dependent variables, but instead tries to map the full population of neurons onto a single unified CS map. It is therefore possible to identify ‘orphan’ clusters of neurons that display variations in their spiking related to uncontrolled variables. While assigning a function to an orphan cluster related to an uncontrolled variable may be challenging, it highlights the potential impact of additional variables in the neuronal population, and presents the opportunity to design additional experiments and potentially identify new classes of neurons. In addition to providing a principled way to determine if a consistent organization of information processing modules can be found across sessions and subjects, we believe that the ability to intuitively visualize relationships within networks of neurons will provide a unique perspective leading to new data-driven hypotheses and experimental refinement.

Several important limitations of SIMNETS are worth noting. First, estimates of similarity using spike train metrics require that the time windows of interest be of equal length, making it difficult to compare neural responses with different time courses. This weakness is common to all trial-averaging models commonly used in the literature that we are aware of. It is important to note that the optimal length of the time window will depend on the specific experimental context and should be adjusted on a case-by-case basis, either based on previous work (as done here) or empirically. Using spike train similarity metrics makes it possible to use relatively long time windows without losing information encoded at fine temporal resolution. This feature makes it easier to adopt time windows that fully encompass behavioral events of interest.

Second, although the SIMNETS framework does not require a priori assumptions about the variables potentially encoded by neural activity, experimental design and data selection will still have a direct effect on the results obtained. For example, a set of neurons identified as a functional subnetwork could separate into smaller groups with different computational properties when additional task conditions are added to the analysis. Thus, the functional properties identified using SIMNETS are only valid within the context of the data examined and may not necessarily extrapolate to different experimental conditions. Ultimately, the computational space of the population is only as rich as the experimental design allows (Gao and Ganguli, 2015).

Third, SIMNETS does not provide a means to assess sources of variance. For example, subsets of neurons may exhibit systematic changes in spike train outputs through an experimental session, reward experience, timing, attention, or many other variables which can be revealed by the neural data; however, it will be up to the experimenter to identify their source. Nevertheless, SIMNETS can guide hypothesis generation for experiments that identify these influences.

Fourth, although peak silhouette values are indicative of the optimal number of clusters according to clustering statistics, they should be taken as a guideline and not an absolute measure of the number of subnetworks. The present study focused on a set of diverse recordings from various sources to highlight the performance of the SIMNETS algorithm for different types of data. An in-depth study of potential subnetwork structure using SIMNETS should include many similar datasets, ideally sampling multiple subjects and sessions to properly assess the number of clusters present across many CS maps before drawing conclusions regarding network architecture.

Fifth, it is possible that neuronal subnetworks are constantly re-arranged depending on changing external or internal signals, e.g., ethological demands, task demands, or task stages, or attentional levels. For example, a neuron could potentially exhibit rapid changes in its computational/functional interrelationships with other neurons during processing epochs that are shorter than the selected analysis time window. The current version of the SIMNETS algorithm was not designed to distinguish between rapidly changing sub-network memberships. However, applying the SIMNETS algorithm multiple times over different epochs can help with determining if different network configurations are engaged at different times.