Spike Train Coactivity Encodes Learned Natural Stimulus Invariances in Songbird Auditory Cortex

Stimuli

Starlings (both male and female) produce long, spectrotemporally rich, and individualized songs composed of repeating, shorter acoustic segments known as “motifs” that are learned over the bird's lifespan. Motifs are on the order of 0.5–1.5 s in duration. In these experiments, a library of 16 6-s-long “pseudo-songs” were used as stimuli. Each pseudo-song was constructed from six 1-s-long motifs manually extracted from a larger library of natural starling song produced by several singers. The first motif of each pseudo-song was an introductory whistle motif, as commonly occurs in natural starling songs. The same whistle was used for all 16 pseudo-songs. After the introductory whistle, all motifs were unique to one pseudo-song (i.e., none of the other motifs occurred in more than one song). The amplitude of each pseudo-song was normalized to 65 dB SPL.

Behavioral training

We trained 4 birds on a two-alternative forced-choice task to recognize four naturalistic pseudo-song stimuli, following previously described operant conditioning procedures (Knudsen and Gentner, 2013). Briefly, the operant apparatus afforded access to three separate response ports. Subjects learned to peck at the center response port to trigger the playback of one of four pseudo-songs associated with two different responses: a peck to the left response port or a peck to the right response port. Pecking at the left response port following the presentation of two songs, or at the right response port following the other two songs, allowed for brief (1–2 s) access to a food hopper. Incorrect responses, pecking the response port not associated with a given song, resulted in a short timeout during which the house light was extinguished and the bird could not feed from the hopper. The 12 remaining pseudo-songs were never presented during operant training for that bird. Assignment of songs for behavioral training was counterbalanced across subjects. The accuracy of the birds over the last 500 trials before neural recording was as follows: B1083, 86.2%; B1056, 95.4%; B1235, 97.6%; B1075, 91.2%. All birds exceeded 90% accuracy during training. Acquisition of the behavior was fast; the number of 500 trial blocks required to first exceed 90% accuracy was as follows: B1083, 30; B1056, 14; B1235, 5; B1075, 13.

Electrophysiology

Once trained, we prepared the birds for physiological extracellular recording from caudo-medial nidopallium (NCM). We anesthetized birds (20% urethane, 7 ml/kg, i.m.), affixed a metal pin to their skull with dental cement, removed an ∼1 mm² window from the top layer of skull, and placed a small craniotomy dorsal to NCM. We inserted a 16- or 32-channel silicon microelectrode (NeuroNexus Technologies) through the craniotomy into the NCM of the head-fixed bird positioned on a foam couch within a sound attenuation chamber (Acoustic Systems). Auditory stimuli were presented from a speaker mounted ∼20 cm above the center of the bird's head, at a mean level of 60 dB SPL measured at 20 cm. The microelectrode was lowered until all recording sites were within NCM, and then allowed to sit for ∼15 min to allow any tissue compression that might have occurred during insertion to relax and stabilize before stimulus presentations began. Auditory activity was confirmed by exposing the birds to non–task-related sounds and observing sound-evoked responses in the raw voltage signal on one or more channels. The 16 pseudo-song stimuli in the block were presented pseudo-randomly until a total of 20 repetitions for each stimulus. For B1083, after the completion of the first block, the electrode was driven deeper into NCM by a distance greater than the electrode length, to obtain another block from a disjoint population. Raw voltages recorded from each microelectrode site were buffered and amplified (20× gain) through a headstage (TBSI), bandpass filtered between 300 Hz and 5 kHz, and amplified (AM Systems, model 3600) and sampled digitally at 20 kHz (16 bit; CED model 1401) via Spike2 software (CED). We saved the full waveforms for offline analysis.

Spike sorting

Extracellular waveforms were spike-sorted offline using the Phy spike sorting framework (Rossant et al., 2016). After automatic clustering, the sort was completed manually using the KwikGUI interface. Clusters with large signal-to-noise ratio and <1% refractory period violations (taken to be 1 ms) were labeled as “Good” clusters and included in the analysis. Because the analyses are meant to quantify population-level structure, sorting priority was given to extracting as many neural signals as possible. Clusters identified automatically were mostly kept separate unless obvious duplicates were observed. Duplicate clusters were defined by overlapping distributions in PCA space, similar waveform shape, and cross-correlations that respected refractory periods. Despite these measures, the set of clusters labeled “Good” are likely to contain multiunit activity as well as single units.

Data processing

Data processing began by converting the neural population activity on each trial into an n × m matrix comprising all the spikes across all the isolated clusters and binned into 10 ms windows with a 5 ms overlap between windows. We used the Perseus persistent homology software for computing Betti numbers from simplicial complexes (Mischaikow and Nanda, 2013). For computing samples from the simplicial configuration model, we used the code described by Young et al. (2017). All other computations were performed using custom-written Python and C code available at http://github.com/theilmbh/NeuralTDA.

Mathematical background

Outside of neuroscience, invariance has been the subject of intense mathematical investigation. In particular, the mathematical field of topology is dedicated to studying the properties of arbitrary spaces that are invariant to different classes of transformations. This research has produced a wealth of techniques and methods for defining and quantifying invariant structures. Recently, these techniques have been increasingly applied to the problem of finding invariant structure in large datasets. Topology takes a global viewpoint, and the “quantities” of interest are entire spaces, rather than single numerical measures. By converting neural population activity into a topological space instead of reducing it to a series of numerical measures, we gain access to the library of techniques topology offers for characterizing invariant structure.

Here, we introduce the basics of the mathematics used for our analysis. More thorough introductions are available elsewhere (Hatcher, 2002), including for the many applications of algebraic topology (Ghrist, 2014). Our analyses work by constructing mathematical objects known as simplicial complexes. Simplicial complexes are topological spaces that are built from discrete building blocks known as simplexes. For each dimension, there is only one prototype simplex. For example, 0-dimensional simplexes are points, 1-dimensional simplexes are line segments, and 2-dimensional simplexes are triangles. Simplexes exist in every dimension, including dimensions >3 where visualization is not possible. For every simplex of dimension d, there is a collection of d-1 simplexes called faces. Consider the tetrahedron, which is the prototype 3-dimensional simplex (or 3-simplex). The faces of the tetrahedron are the four triangles, which themselves are 2-simplexes. The “faces” of the triangles are the 6 edges, which are 1-simplexes, and so on.

Simplicial complexes are built by “joining” simplexes along common faces. That information is carried through the boundary operators, which are linear maps encoding which d-simplexes are attached to which (d-1) simplexes. Given a simplicial complex, we construct a series of vector spaces, one for each dimension of simplex in the complex. These vector spaces are called “chain groups.” The basis vectors of these vector spaces are taken to be, formally, the simplexes in the various dimensions. For example, a graph has a chain group in dimension 0 consisting of all formal linear combinations of vertices, and a chain group in dimension 1 consisting of all formal linear combinations of edges. The boundary maps are linear maps (matrices) between chain groups of adjacent dimensions, mapping the dimension d chain group to the d-1 chain group. These matrices contain all the structure of the simplicial complex, which allows for all the computations described in this work to be performed using the standard, highly optimized numerical linear algebra routines available in scientific programming packages, such as SciPy.

Constructing simplicial complexes for neural data

To construct the simplicial complex associated to a population spike train, we follow methods similar to those developed by Curto and Itskov (2008). For each presentation of a stimulus, the neural response was divided into time bins of width dt and a percentage overlap between adjacent bins, both free parameters for the analysis. Unless otherwise noted, for all analyses, dt was 10 ms with 50% overlap, to preserve continuity of the spike trains. Each spike in the time period was assigned to its associated time bins. The result is an N_cells × N_bins × N_trials multidimensional array D, representing the population activity, where each element is a spike count in the time window dt. The spike counts were then divided by the window width to determine a “firing rate” for that cell in that time bin. The time average was taken to yield an N_cells × N_trials matrix of average instantaneous firing rates for each cell for each trial. Then, the multidimensional array D was thresholded by determining which bins had a firing rate greater than some threshold value times the average firing rate. This yielded an N_cells × N_bins × N_trials binary array B that represents significantly active cells. Unless otherwise noted, the threshold value of 4 times the average firing rate of the unit was used. Since the average firing rate of the units was low, the threshold value did not make a significant impact on the results.

For each trial and for each time bin in B, the N_cell-long binary vector V was used to create a list of cells that were active in that time bin. This was repeated across all time bins to create a list of “cell groups” active during the trial, much like the cell groups defined by Curto and Itskov (2008). This list of cell groups was taken to define the maximal simplexes for the simplicial complex. Using algorithms derived from published work (Kaczynski et al., 2006), this list of cell groups was turned into a list of lists of basis vectors for the chain complex. The n-th sublist represents the chain complex basis vectors corresponding to the n-dimensional simplexes in the simplicial complex. We also computed the matrices that represent the boundary operators in the chain complex.

Betti curves

Homology is one of the most basic topological invariants one can use to describe a topological space. Loosely speaking, homology counts the number of n-dimensional holes in the space. These holes are described purely algebraically by computing how n-dimensional cycles (e.g., a circle) are not the boundary of an n + 1-dimensional simplex. For each dimension n, the number of such holes is called the n-th Betti number. The reason that Betti numbers are a useful characterization of the topological spaces we construct is that any two topological spaces that differ in their Betti numbers cannot be topologically equivalent. Furthermore, topological holes interpreted neurally represent “gaps” in the coactivation pattern (which neurons do not fire together) and thus carry information about how the population activity is structured. We summarized the topology of neurally derived simplicial complexes by computing these Betti numbers across time to yield Betti curves.

Figure 1d illustrates the temporal filtration. From the start of the stimulus, cell groups from each time bin are added to the growing simplicial complex. The Betti numbers for the entire complex are computed to yield the value of the Betti curve for that bin. Then, the next bin is added, and the Betti numbers are recomputed. Simplicial complexes were fed into the Perseus persistent homology program (Mischaikow and Nanda, 2013; Nanda, 2013), which computed the sequence of Betti numbers in the evolving simplicial complex. In the language of applied topology, we computed the homology of the filtered simplicial complex using time as the filtration parameter (Giusti et al., 2016). Betti curves were computed for each dimension and for each trial of a given stimulus. Because the Betti numbers are discrete, we interpolated the output from Perseus using step functions. This allowed us to average the Betti curves across trials. We varied the time bin width dt between 5 and 250 ms and observed that the shape of the Betti curves was consistent for values of dt between 5 and 50 ms. For values of dt > 50 ms, the Betti curves looked like low-pass-filtered versions of the small-dt curves. Varying the percentage overlap from 0% to 50% did not appear to significantly alter the Betti curves.

The Betti curves reflect the accumulated effects of moment-to-moment spatiotemporal structure in the population activity. To serve as a control, we used various shuffling procedures to break this structure. We computed “fully shuffled” Betti curves by taking the neural response matrix for each trial and shuffling each row independently across columns. Since each row corresponds to a cell's response, this approach shuffles each cell's response independently across trials and across cells. The effect is to break the temporal correlations between cells while preserving the total number of spikes from a given cell in a given trial. We computed trial-shuffled Betti curves by using a similar shuffling procedure; but instead of shuffling a cell's response across time independently for each trial and each cell, we shuffled each response across trials independently for each cell and each time bin.

Spatiotemporal structure in population responses to stimuli reflects two kinds of correlation. The first is correlation between individual cell spike trains and individual stimuli. We call this the stimulus specificity of single-unit responses. The other kind of correlation is the set of n-wise correlations between neurons on a single trial across time (sometimes referred to as the noise correlation). The fully shuffled condition destroys both these forms of correlation, by assigning different single-unit spike trains, and therefore a random n-wise correlation structure, to each neuron on each presentation of each stimulus. To dissociate these two correlation sources, we created an additional shuffle of the population data, termed the “within-stimulus mask shuffle,” in which each stimulus is associated with a unique time bin-to-time bin permutation mask. We then shuffle the spiking response of each neuron in time individually and independently within a single trial for each stimulus, applying the same bin-to-bin permutation to all trials of that same stimulus. This creates a new, independent, stimulus-specific response for each neuron, and replaces the empirical n-wise population correlations with a random pattern of coactivity that repeats across repetitions of each stimulus.

Stimulus specificity

We analyzed the stimulus specificity of the empirical and shuffled Betti curves by first computing the distribution of the correlations between the Betti curves of pairs of trials. Then we compared the distributions of correlations from a pair of trials within a single stimulus to correlations from pairs of trials between different stimuli. The significance of these differences was computed using a linear mixed-effects model (Lindstrom and Bates, 1988). The model predicted the correlation between Betti curves on pairs of trials treating the within-stimulus versus between-stimulus distinction as a random effect and treating the average firing rate and an interaction of firing rate and pair type as fixed effects. We included firing rate effects since the shuffles preserve average firing rate but destroy precise temporal correlations.

Simplicial configuration model

To test whether the topology of the neurally derived simplicial complexes is likely to be seen in a random simplicial complex, we used the simplicial configuration model (Young et al., 2017). Because of the computational constraints of the configuration model, for each stimulus, we averaged the neural activity over trials to yield an average population spike train. The simplicial complex for this spike train was computed and seeded the simplicial configuration model software (Young et al., 2017). This algorithm produces a Markov chain that yields samples from the uniform distribution over all simplicial complexes with local connectivity structure identical to the seed complex derived from neural activity. For each sample, the Betti numbers were computed using Perseus. The distribution of these Betti numbers serves as a null distribution against which Betti numbers derived from the empirical data can be compared. We interpret data having Betti numbers that fall well outside the null distribution as containing nonrandom topological structure. To quantify how far outside the null distribution an observed Betti number falls, we calculated empirical p values using the formula as follows (Davison and Hinkley, 1997): (1)

Where b_i represents the Betti number for the i-th sample from the simplicial configuration model, N is the number of samples from the simplicial configuration model, and b_e is the empirically observed Betti number.

Simplicial Laplacian

We generalized the methods of De Domenico and Biamonte (2016) to define information-theoretic quantities related to simplicial complexes. Equipped with the boundary operator matrices , the simplicial Laplacians (Horak and Jost, 2013) were computed as follows: (2) where * indicates the adjoint, which in these analyses means the matrix transpose. Given a simplicial Laplacian , the related density matrix was computed as follows : (3) where β is a free parameter. Given a density matrix expressed in its eigen-basis, the eigenvalues represent the probability distribution over eigenstates with maximum entropy under the Hamiltonian given by the Laplacian. We can define the Kullback-Leibler (KL) divergence between any two such density matrices, ρand by first diagonalizing both ρand σ and then sorting each vector of eigenvalues by magnitude. The KL divergence is then defined by the following: (4)

Where is the i-th eigenvalue. The Jensen-Shannon (JS) divergence was defined using the following: (5) (6)

The advantage of sorting the eigenvalues is that the divergence becomes invariant to the labels assigned to individual neurons. This is desirable because these measures should depend only the neural activity and not the specific labels assigned to neurons, analogous to how different maps of the world describe the same geography with different coordinates. This definition of JS and KL divergence ignores the eigenvectors of ρ and σ. Indeed, ρ and σ may not be simultaneously diagonalizable (or even of the same dimension). Rectifying these issues requires a choice in the definition of the metrics, which like any other metric is justified by its ultimate use. We explain our choices in the following paragraph. Regardless, the eigenspectra of ρ and σ do form discrete probability distributions, and the definitions of KL and JS divergence here agree with the usual definitions of these divergences on discrete probability distributions.

One limitation in the naive generalization of the foregoing metrics to neural data are that they rely on the density matrices being square and of the same dimension. For real data, however, different stimuli or repetitions of the same stimulus often evoke different numbers of active neurons, giving rise to chain groups, boundary operators, and density matrices of different dimension. To address this, given two Laplacians of different dimensions, we expanded the dimension of the smaller Laplacian by padding with zeros. An alternative approach is to collapse the two simplicial complexes together along their common simplexes and derive “masked” boundary operators that only operate on the simplexes within the original, individual complexes, and then define the Laplacians using these masked boundary operators. Both the “zero-padding” and the “masking” approaches yielded similar results for our datasets. We report results from the computationally simpler method of zero-padding. In all analyses, the simplicial Laplacian was computed from the final simplicial complex constructed from the population response to an entire trial.

Beta parameter

The free parameter, β, that appears in the expression for the density matrix originates in statistical physics, where it is interpreted as “inverse temperature,” i.e., β is proportional to 1/temperature. The significance of β in the context of network theory is actively under study (Nicolini et al., 2018). Here we interpret β by noting that it acts to normalize the Laplacian matrix, and by extension, the eigenvalues of the Laplacian. In this sense, β acts as a “scale parameter” that sets the scale of the spectral features of interest. The eigenvalues of the Laplacian are always non-negative; and if β is large, then large eigenvalues are significantly damped by the negative exponential. Heuristically, increasing β decreases the proportion of larger eigenvalues that contribute to the quantities of interest. In our case, if β is too large, too little of the spectrum will be relevant for the JS or KL divergences, and we reasoned the results would be less interpretable. To confirm this exponential damping, we computed the proportion of eigenvalues above a threshold of 1e-14 as β varied from 0.1-100. Generally, this proportion started decreasing from 100% with β ≅ 1, to below 20% when β > 10. We set β = 1 in all of our analyses, but the main conclusions do not change as along as β was below the upper bound of 10.

Simulated spiking populations

We used synthetic spike trains to examine how well the proposed metrics reveal invariant properties between neural populations. For “target” spike trains, we simulated a population of 20 Poisson spiking neurons for 1000 time-steps with a rate parameter selected from the set (0.01, 0.02, 0.03, 0.04). We generated “test” spike trains by choosing 50 rate parameters evenly in the range of [0.001, 0.1] and simulating 25 population spike trains for each rate. We computed the simplicial KL divergence between each test spike train and the target spike train, then computed the average divergence between the target and the model at the given rate parameter.

We applied the same procedure to simulate heterogeneous Poisson populations, except that half of the neurons were simulated with rate parameter 0.02 and the other with 0.05. The assignment of neurons to rate-parameter subgroups was random. The above-described parameter sweep procedure was repeated over the two population rate parameters.

Simulated environments

We simulated physical environments similarly to the environments used in Curto and Itskov's original work (Curto and Itskov, 2008). We randomly placed 0-4 holes of radius 0.3 m inside a 2 × 2 m arena, ensuring hole centers did not overlap. We then randomly distributed 100 place fields of radius 0.2 m in the environment using an algorithm that ensured that the environment was fully tiled. To model movement, we constructed a random walk trajectory simulating a 10 min traverse of the environment sampled at 10 samples per second to ensure the whole arena was explored. We constructed spike trains by simulating Poisson spike trains at a constant rate for each neuron while the random walk trajectory was within the cell's place field, and silent outside of it.

We constructed five different environments with each specified number of holes. For each environment, 10 random walks were simulated. The pairwise 1-JS divergence between pairs of individual trials was computed by taking the simulated spike trains' 1-Laplacians and computing the spectral JS divergence as defined above. We averaged the 1-JS divergences over all walks for each environment to yield the plot in Figure 5a. Then, the JS divergences for each of the 25 environment pairs with given numbers of holes were averaged to yield the plot in Figure 5b.

In vivo JS divergences

Data from trained, anesthetized birds was recorded and preprocessed to yield population spike trains as described earlier. The JS divergence analyses were restricted to the four learned stimuli and four other unfamiliar stimuli. For an arbitrary pair of responses from two separate stimulus presentations, the 1-Laplacian of each response was computed by first converting the neural data into a simplicial complex as described above, using a bin size of 10 ms with 5 ms overlaps. Then, the JS divergences between the 1-Laplacians were computed as described above. It is important to note that the JS divergences were computed using the final simplicial complex from each trial. As a result, it may “miss” homological features that do not persist through the whole trial. We do not consider this a limitation in the present study, since the final simplicial complex reflects the aggregate effects of temporal coactivations through the entire trial, and thus tells us something about the population response to the entire stimulus considered as a unit. Nevertheless, it will be valuable for future work to explore the properties of these features that may be missed by the current analysis.

Correlation-based population response similarity

As an alternative to the simplicial Laplacian spectral entropy (SLSE) measure of population similarity, we also compute a correlation-based measure to quantify the trial-to-trial similarity in the response of a population, C(A,B), where A and B are the responses of the same population on different trials responding to either the same stimulus or to different stimuli. We first smooth the population spike trains with a Gaussian kernel with a width of 10 ms, then compute the correlation coefficients between the spike train of a neuron in response A and the spike train from the same neuron in response B, then average across all the correlation coefficients to a yield a raw scalar estimate of similarity, C, between the two populations. We use 1 – C in place of the raw correlation measure so that smaller values of both correlation and JS divergence indicated more similar spike trains, such that 1 – C = 0 when the response of every single neuron in the population for population response A is perfectly correlated with its corresponding response for population response B. Comparing this similarity measure with the SLSE measures allows us to test whether physically dissimilar responses (in terms of the actual spike trains) may evince similar topologies.

We also compute the distance between pairs of population responses using covariance matrices. For each trial, we calculated the pairwise covariance matrix between Gaussian-smoothed spike trains, again with a 10 ms window. This yields an N_cell × N_cell matrix for each trial in which each i,j entry is the temporal covariance between the smoothed spike trains from neurons i and j on that trial. We obtained a distance between pairs of trials by computing the Frobenius norm of the difference between each trial's covariance matrix. The distance matrix was normalized by its maximum value over all pairs of trials to yield a normalized distance measure. This measure places two trials close to each other if the second-order (pairwise) correlation structure in the population for response A is similar to the second-order correlation structure for response B.

The correlation between population responses with relabeled neurons was obtained by shuffling the rows of the Gaussian smoothed population response matrix and repeating the population correlation computation described above on the shuffled matrices.

JS divergence and correlation measures were compared using the relative accuracy of logistic regressions. We fit logistic regressions to predict whether a pair of responses to different stimuli came from stimuli that belonged to either the same or different behaviorally trained class based on the value of either the JS divergence or the specific correlation measure of interest between the responses. Regressions were fit to a random subset of 80% of the data and tested on the remaining 20%. A total of 240 separate regressions were performed for each condition.

We measured the similarity of topological structure across disjoint subpopulations by randomly splitting the full dataset from B1083 Population 1 into two disjoint subgroups each with 50 neurons. Then we conducted an identical analysis to the logistic regression described in the preceding paragraph but using the JS divergence between responses from the two disjoint populations. We repeated this 40 times with a different (randomly chosen) split of the original population. We assessed the difference between the means of the accuracies of the regressions fit to the JS divergence or correlations with a Gaussian GLM from the statsmodels python library with the type of measure (JS divergence/correlation) and the type of shuffle (shuffle/no shuffle) as regressors, and included an interaction term between the shuffle and measure.

Linear mixed-effects model

We used a linear mixed-effects model to examine the stimulus specificity of the empirically derived Betti curves associated with each stimulus, and those derived from the various shuffled versions of the spiking responses. The regression models predicted the correlation between Betti curves in a given dimension between pairs of trials, as a function of whether the pair came from two trials of a single stimulus or two trials of different stimuli. The model included population identity as a random effect. The model was implemented using MixedLM from the statsmodels python library. We included firing rate (and its interaction with stimuli) as a separate fixed effect in the model and exclude this from all reported effects on Betti curve correlations tied to stimulus specificity. Overall, the mean firing rate on a pair of trials tends to increase the strength of the correlation between Betti curves on those trials, in part because it tends to smooth over small moment-to-moment variations in a manner equivalent to computing simplicial complexes from spike trains quantized into larger time bins.