dPCA
An implementation of demixed Principal Component Analysis (a supervised linear dimensionality reduction technique)
view repo
Neurons in higher cortical areas, such as the prefrontal cortex, are known to be tuned to a variety of sensory and motor variables. The resulting diversity of neural tuning often obscures the represented information. Here we introduce a novel dimensionality reduction technique, demixed principal component analysis (dPCA), which automatically discovers and highlights the essential features in complex population activities. We reanalyze population data from the prefrontal areas of rats and monkeys performing a variety of working memory and decision-making tasks. In each case, dPCA summarizes the relevant features of the population response in a single figure. The population activity is decomposed into a few demixed components that capture most of the variance in the data and that highlight dynamic tuning of the population to various task parameters, such as stimuli, decisions, rewards, etc. Moreover, dPCA reveals strong, condition-independent components of the population activity that remain unnoticed with conventional approaches.
READ FULL TEXT VIEW PDF
All important decisions about the variety of rice grain end product are ...
read it
The allocation of a (treatment) condition-effect on the wrong principal
...
read it
In areas of application, including actuarial science and demography, it ...
read it
The human brain contains approximately 10^9 neurons, each with
approxima...
read it
Independent component analysis (ICA) has often been used as a tool to mo...
read it
We have developed an efficient information-maximization method for compu...
read it
Motivation
The genotype assignment problem consists of predicting, from ...
read it
An implementation of demixed Principal Component Analysis (a supervised linear dimensionality reduction technique)
Preprocessing and analysis scripts for the dPCA paper (eLife 2016)
In many state of the art experiments, a subject, such as a rat or a monkey, performs a behavioral task while the activity of tens to hundreds of neurons in the animal’s brain is monitored using electrophysiological or imaging techniques. The common goal of these studies is to relate the external task parameters, such as stimuli, rewards, or the animal’s actions, to the internal neural activity, and to then draw conclusions about brain function. This approach has typically relied on the analysis of single neuron recordings. However, as soon as hundreds of neurons are taken into account, the complexity of the recorded data poses a fundamental challenge in itself. This problem has been particularly severe in higher-order areas such as the prefrontal cortex, where neural responses display a baffling heterogeneity, even if animals are carrying out rather simple tasks (Brody et al., 2003; Machens, 2010; Mante et al., 2013; Cunningham and Yu, 2014).
Traditionally, this heterogeneity has often been ignored. In neurophysiological studies, it is common practice to pre-select cells based on particular criteria, such as responsiveness to the same stimulus, and to then average the firing rates of the pre-selected cells. This practice eliminates much of the richness of single-cell activities, similar to imaging techniques with low spatial resolution, such as MEG, EEG, or fMRI. While population averages can identify the information that higher-order areas process, they ignore how exactly that information is represented on the neuronal level (Wohrer et al., 2013). Indeed, most neurons in higher cortical areas will typically encode several task parameters simultaneously, and therefore display what has been termed “mixed selectivity” (Rigotti et al., 2013).
Instead of looking at single neurons and selecting from or averaging over a population of neurons, neural population recordings can be analyzed using dimensionality reduction methods (for a review, see Cunningham and Yu, 2014). In recent years, several such methods have been developed that are specifically targeted to electrophysiological data, taking into account the binary nature of spike trains (Yu et al., 2009; Pfau et al., 2013), or the dynamical properties of the population response (Buesing et al., 2012; Churchland et al., 2012). However, these approaches reduce the dimensionality of the data without taking task parameters, i.e., sensory and motor variables controlled or monitored by the experimenter, into account. Consequently, mixed selectivity remains in the data even after the dimensionality reduction step, impeding interpretation of the results.
A few recent studies have sought to solve these problems by developing methods that reduce the dimensionality of the data in a way that is informed by the task parameters (Machens, 2010; Machens et al., 2010; Brendel et al., 2011; Mante et al., 2013). One possibility is to adopt a parametric approach, i.e. to assume a specific (e.g. linear) dependency of the firing rates on the task parameters, and then use regression to construct variables that demix the task components (Mante et al., 2013). While this method can help to sort out the mixed selectivities, it still runs the risk of missing important structures in the data if the neural activities do not conform to the dependency assumptions (e.g. because of nonlinearities).
Here, we follow an alternative route by developing an unbiased dimensionality reduction technique that fulfills two constraints. It aims to find a decomposition of the data into latent components that (a) are easily interpretable with respect to the experimentally controlled and monitored task parameters; and (b) preserve the original data as much as possible, ensuring that no valuable information is thrown away. Our method, which we term demixed principal component analysis (dPCA), improves our earlier methodological work (Machens, 2010; Machens et al., 2010; Brendel et al., 2011) by eliminating unnecessary orthogonality constraints on the decomposition. In contrast to several recently suggested algorithms for decomposing firing rates of individual neurons into demixed parts (Pagan and Rust, 2014; Park et al., 2014), our focus is on dimensionality reduction.
There is no a priori guarantee that neural population activities can be linearly demixed into latent variables that reflect individual task parameters. Nevertheless, we applied dPCA to spike train recordings from monkey prefrontal cortex (PFC) (Romo et al., 1999; Qi et al., 2011) and from rat orbitofrontal cortex (OFC) (Feierstein et al., 2006; Kepecs et al., 2008), and obtained remarkably successful demixing. In each case, dPCA automatically summarizes all the important, previously described features of the population activity in a single figure. Importantly, our method provides an easy visual comparison of complex population data across different tasks and brain areas, which allows us to highlight both similarities and differences in the neural activities. Demixed PCA also reveals several hitherto largely ignored features of population data: (1) most of the activity in these datasets is not related to any of the controlled task parameters, but depends only on time (“condition-independent activity”); (2) all measured task parameters can be extracted with essentially orthogonal linear readouts; and (3) task-relevant information is shifted around in the neural population space, moving from one component to another during the course of the trial.
We illustrate our method with a toy example (Figure 1). Consider a standard experimental paradigm in which a subject is first presented with a stimulus and then reports a binary decision. The firing rates of two simulated neurons are shown in Figure 1A. The first neuron’s firing rate changes with time, with stimulus (at the beginning of the trial), and with the animal’s decision (at the end of the trial). The second neuron’s firing rate also changes with time, but only depends on the subject’s decision (not on the stimulus). As time progresses within a trial, the joint activity of the two neurons traces out a trajectory in the space of firing rates (Figure 1B; decision not shown for visual clarity). More generally, firing rates of a population of neurons at any moment in time is represented by a point in an -dimensional space.
One of the standard approaches to reduce the dimensionality of complex, high-dimensional datasets is principal component analysis (PCA). In PCA, the original data (here, firing rates of the neurons) are linearly transformed into several latent variables, called principal components. Each principal component is therefore a linear combination, or weighted average, of the actual neurons. We interpret these principal components as linear “read-outs” from the full population (Figure 1C, “decoder” ). This transformation can be viewed as a compression step, as the number of informative latent variables is usually much smaller than the original dimensionality. The resulting principal components can then be de-compressed with another linear mapping (“encoder” ), approximately reconstructing the original data. Geometrically, when applied to a cloud of data points (Figure 1D), PCA constructs a set of directions (principal axes) that consecutively maximize the variance of the projections; these axes form both the decoder and the encoder.
More precisely, given a data matrix
, in which each row contains the firing rates of one neuron for all task conditions, PCA finds decoder and encoder that minimize the loss function
under the constraint that the principal axes are normalized and orthogonal, and therefore (Hastie et al., 2009, section 14.5). However, applying PCA to a neural dataset like the one sketched in Figure 1A generally results in principal components exhibiting mixed selectivity: many components will vary with time, stimulus, and decision (see examples below). Indeed, information about the task parameters does not enter the loss function.
To circumvent this problem, we sought to develop a new method by introducing an additional constraint: the latent variables should not only compress the data optimally, but also demix the selectivities (Figure 1E). As in PCA, compression is achieved with a linear mapping and decompression with a linear mapping . Unlike PCA, the decoding axes are not constrained to be orthogonal to each other and may have to be non-orthogonal, in order to comply with demixing. The geometrical intuition is presented in Figure 1F. Here, the same cloud of data points as in Figures 1B and 1D has labels corresponding to time and stimulus. When we project the data onto the time-decoding axis, information about time is preserved while information about stimulus is lost, and vice versa for the stimulus-decoding axis. Hence, projections on the decoding axes demix the selectivities. Since the data have been projected onto a non-orthogonal (affine) coordinate system, the optimal reconstruction of the data requires a de-compression along a different set of axes, the encoding axes (which implies that , the pseudo-inverse of , see Methods).
Given the novel demixing constraint, we term our method demixed PCA (dPCA), and provide a detailed mathematical account in the Methods. Briefly, in the toy example of Figure 1, dPCA first splits the data matrix into a stimulus-varying part and a time-varying part , so that , and then finds separate decoder and encoder matrices for stimulus and time by minimizing the loss function
See Methods for a more general case.
We first applied dPCA to recordings from the PFC of monkeys performing a somatosensory working memory task (Romo et al., 1999; Brody et al., 2003). In this task, monkeys were required to discriminate two vibratory stimuli presented to the fingertip. The stimuli were separated by a 3 s delay, and the monkeys had to report which stimulus had a higher frequency by pressing one of two available buttons (Figure 2A). Here, we focus on the activity of 832 prefrontal neurons from two animals, see Methods. For each neuron, we obtained the average time-dependent firing rate (also known as peri-stimulus time histogram, PSTH) in 12 conditions: 6 values of stimulus frequency F1, each paired with 2 possible decisions of the monkey. We only analyzed correct trials, so that the decision of the monkey always coincided with the actual trial type (F2 F1 or F1 F2).
As is typical for PFC, each neuron has a distinct response pattern and many neurons show mixed selectivity (Figure S1). Several previous studies have sought to make sense of these heterogeneous response patterns by separately analyzing different task periods, such as the stimulation and delay periods (Romo et al., 1999; Brody et al., 2003; Machens et al., 2010; Barak et al., 2010), the decision period (Jun et al., 2010), or both (Hernández et al., 2010). With dPCA, however, we can summarize the main features of the neural activity across the whole trial in a single figure (Figure 2B). Just as in PCA, we can think of these demixed principal components as the “building blocks” of the observed neural activity, in that the activity of each single neuron is a linear combination (weighted average) of these components. These building blocks come in four distinct categories: some are condition-independent (Figure 2B, top row); some depend only on stimulus F1 (second row); some depend only on decision (third row); and some depend on stimulus and decision together (bottom row). The components in the first three categories demix the parameter dependencies, which is exactly what dPCA aimed for. We stress that applying standard PCA to the same dataset leads to components that strongly exhibit mixed selectivity (Figure S2; see Methods for quantification).
What can we learn about the population activity from the dPCA decomposition? First, we find that information about the stimulus and the decision can be fully demixed at the population level, even though it is mixed at the level of individual neurons. Importantly, the demixing procedure is linear, so that all components could in principle be retrieved by the nervous system through synaptic weighting and dendritic integration of the neurons’ firing rates. Note that our ability to linearly demix these different types of information is a non-trivial feature of the data.
Second, we see that most of the variance of the neural activity is captured by the condition-independent components (Figure 2B, top row) that together amount to 90% of the signal variance (see the pie chart in Figure 2D and Methods). These components capture the temporal modulations of the neural activity throughout the trial, irrespective of the task condition. Their striking dominance in the data may come as a surprise, as such condition-independent components are usually not analyzed or shown. Previously, single-cell activity in the delay period related to these components has often been dubbed “ramping activity” or “climbing activity” (Durstewitz, 2004; Rainer et al., 1999; Komura et al., 2001; Brody et al., 2003; Janssen and Shadlen, 2005; Machens et al., 2010); however, analysis of the full trial here shows that the temporal modulations of the neural activity are far more complex and varied. Indeed, they spread out over many components. The origin of these components can only be conjectured, see Discussion.
Third, our analysis also captures the major findings previously obtained with these data: the persistence of the F1 tuning during the delay period — component #6 (Romo et al., 1999; Machens et al., 2005), the temporal dynamics of short-term memory — components ##6, 10, 11 (Brody et al., 2003; Machens et al., 2010; Barak et al., 2010), the ramping activities in the delay period — components ##1–3, 6 (Singh and Eliasmith, 2006; Machens et al., 2010); and pronounced decision-related activities — component #5 (Jun et al., 2010). We note that the decision components resemble derivatives of each other; these higher-order derivatives likely arise due to slight variations in the timing of responses across neurons (Figure S3).
Fourth, we observe that the three stimulus components show the same monotonic (“rainbow”) tuning but occupy three distinct time periods: one is active during the S1 period (#10), one during the delay period (#6), and one during the S2 period (#11). Information about the stimulus is therefore shifted around in the high-dimensional firing rate space. Indeed, if the stimulus were always encoded in the same subspace, there would only be a single stimulus component. In other words, if we perform dPCA analysis in a short sliding time window, then the main stimulus axis will rotate with time instead of remaining constant.
We furthermore note several important technical points. First, the overall variance explained by the dPCA components (Figure 2C, red line) is very close to the overall variance captured by the PCA components (Figure 2C, black line). This means that by imposing the demixing constraint we did not lose much of the variance, and so the population activity is accurately represented by the obtained dPCA components (Figure 2B). Second, as noted above, the demixed principal axes are not constrained to be orthogonal. Nonetheless, most of them are quite close to orthogonal (Figure 2E, upper-right triangle). There are a few notable exceptions, which we revisit in the Discussion section. Third, pairwise correlations between components are all close to zero (Figure 2E, lower-left triangle), as should be expected as the components are considered to represent independent signals.
We next applied dPCA to recordings from the PFC of monkeys performing a visuospatial working memory task (Qi et al., 2011; Meyer et al., 2011; Qi et al., 2012). In this task, monkeys first fixated at a small white square at the centre of a screen, when a square S1 appeared in one of the eight locations around the centre (Figure 3A). After a 1.5 s delay, a second square S2 appeared in either the same (“match”) or the opposite (“non-match”) location. Following another 1.5 s delay, a green and a blue choice targets appeared in locations orthogonal to the earlier presented stimuli. Monkeys had to saccade to the green target to report a match condition, and to the blue one to report a non-match.
We analyzed the activity of 956 neurons recorded in the lateral PFC of two monkeys performing this task. Proceeding exactly as before, we obtained the average time-dependent firing rate of each neuron for each condition. Following the original studies, we eliminated the trivial rotational symmetry of the task by collapsing the eight possible stimulus locations into five locations that are defined with respect to the preferred direction of each neuron (0, 45, 90, 135, or 180 away from the preferred direction, see Methods). As a consequence, we obtained 10 conditions: 5 possible stimulus locations, each paired with 2 possible decisions of the monkey. We again only analyzed correct trials, so that the decision of the monkey always coincided with the actual trial type.
The dPCA results are shown in Figure 3. Just as before, the components fall into four categories: they can be condition-independent (Figure 3B, top row); dependent only on stimulus location (second row); only on decision (third row); or dependent on stimulus-decision interactions (bottom row).
We note several similarities to the somatosensory working memory task. First, stimulus and decision can be separated at the population level just as easily as before, despite being intermingled at the single-neuron level. Second, most of the variance in the neural activity is again captured by the condition-independent components (Figure 3B, top row, and 3D, pie chart). Third, all stimulus components are active in different time periods, and the same is true for the decision, indicating rotation of the stimulus and decision representations in the firing rate space. Fourth, and maybe most surprisingly, we find that the overall structure of population activity up to the second stimulus (S2) is almost identical to that found in the somatosensory task. In both cases, the leading condition-independent components are quite similar (compare Figure 2B with Figure 3B). Furthermore, in both tasks we find a stimulus component with strong activity during the F1/S1 period, a stimulus component with persistent activity during the delay period, and a decision component with activity during the F2/S2 period.
One notable difference between Figures 2 and 3 is the presence of strong interaction components in Figure 3B. However, interaction components in Figure 3B are in fact stimulus components in disguise. In match trials, S2 and S1 appear at the same location, and in non-match trials at opposite locations. Information about S2 is therefore given by a non-linear function of stimulus S1 and the trial type (i.e. decision), which is here captured by the interaction components.
The differences in the population activity between these two tasks seem to stem mostly from the specifics of the tasks, such as the existence of a second delay period in the visuospatial working memory task. There are also differences in the overall amount of power allocated to different components, with the stimulus (and interaction) components dominating in the visuospatial working memory task, whereas stimulus and decision components are more equally represented in the somatosensory task. However, the overall structure of the data is surprisingly similar.
Here again, our analysis summarizes previous findings obtained with this dataset. For instance, the first and the second decision components show tuning to the match/non-match decision in the delay period between S2 and the cue (first component) and during the S2 period (second component). Using these components as fixed linear decoders, we achieve cross-validated single-trial classification accuracy of match vs. non-match of 75% for (Figure S4), which is approximately equal to the state-of-the-art classification performance reported previously (Meyers et al., 2012).
Constantinidis et al. have also recorded population activity in PFC before starting the training (both S1 and S2 stimuli were presented exactly as above, but there were no cues displayed and no decision required). When analyzing this pre-training population activity with dPCA, the first stimulus and the first interaction components come out close to the ones shown on Figure 3, but there are no decision and no “memory” components present (Figure S5), in line with previous findings (Meyers et al., 2012). These task-specific components appear in the population activity only after extensive training.
Next, we applied dPCA to recordings from the OFC of rats performing an odour discrimination task (Feierstein et al., 2006). This behavioral task differs in two crucial aspects from the previously considered tasks: it requires no active storage of a stimulus, and it is self-paced. To start a trial, rats entered an odour port, which triggered delivery of an odour with a random delay of 0.2–0.5 s. Each odour was uniquely associated with one of the two available water ports, located to the left and to the right from the odour port (Figure 4A). Rats could sample the odour for as long as they wanted up to 1 s, and then had to move to one of the water ports. If they chose the correct water port, reward was delivered following an anticipation period of random length (0.2–0.5 s).
We analyzed the activity of 437 neurons recorded in five rats in 4 conditions: 2 stimuli (left and right) each paired with 2 decisions (left and right). Two of these conditions correspond to correct (rewarded) trials, and two correspond to error (unrewarded) trials. Since the task was self-paced, each trial had a different length; in order to align events across trials, we restretched the firing rates in each trial (see Methods). Alignment methods without restretching led to similar results (Figure S6).
Just as neurons from monkey PFC, neurons in rat OFC exhibit diverse firing patterns and mixed selectivity (Feierstein et al., 2006). Nonetheless, dPCA is able to demix the population activity, resulting in the condition-independent, stimulus, decision, and interaction components (Figure 4). In this dataset, interaction components separate rewarded and unrewarded conditions (thick and thin lines on Figure 4B, bottom row), i.e., correspond to neurons tuned either to reward, or to the absence of reward.
We note several similarities to the monkey PFC data. The largest part of the total variance is again due to the condition-independent components (over 60% of the total variance falls into the first three components). Also, for both stimulus and decision, the main components are localized in distinct time periods, meaning that information about the respective parameters is shifted around in firing rate space.
The overall pattern of neural tuning across task epochs that we present here agrees with the findings of the original study
(Feierstein et al., 2006). Interaction components are by far the most prominent among all the condition-dependent components, corresponding to the observation that many neurons are tuned to presence/absence of reward. Decision components come next, with a caveat that decision information may also reflect the rat’s movement direction and/or position, as was pointed out previously (Feierstein et al., 2006). Stimulus components are less prominent, but nevertheless show clear stimulus tuning, demonstrating that even in error trials there is reliable information about stimulus identity in the population activity.Finally, we note that the first interaction component (#4) shows significant tuning to reward already in the anticipation period. In other words, neurons tuned to presence/absence of reward start firing before the reward delivery (or, on error trials, before the reward could have been delivered). We return to this observation in the next section.
Kepecs et al. (2008) extended the experiment of Feierstein et al. (2006) by using odour mixtures instead of pure odours, thereby varying the difficulty of each trial (Uchida and Mainen, 2003; Kepecs et al., 2008). In each trial, rats experienced mixtures of two fixed odours with different proportions (Figure 5A). Left choices were rewarded if the proportion of the “left” odour was above 50%, and right choices otherwise. Furthermore, the waiting time until reward delivery (anticipation period) was increased to 0.3–2 s.
We analyzed the activity of 214 OFC neurons from three rats recorded in 8 conditions, corresponding to 4 odour mixtures, each paired with 2 decisions (left and right). During the presentation of pure odours (100% right and 100% left) rats made essentially no mistakes, and so we had to exclude these data from the dPCA computations (which require that all parameter combinations are present, see Methods). Nevertheless, we displayed these additional 2 conditions in Figure 5.
The dPCA components shown in Figure 5B are similar to those presented above in Figure 4B. The condition-independent components capture most of the total variance; the interaction components (corresponding to the reward) are most prominent among the condition-dependent ones; the decision components are somewhat weaker and show tuning to the rat’s decision, starting from the odour port exit and throughout the rest of the trial; and the stimulus components are even weaker, but again reliable.
Here again, some of the interaction components (especially the second one, #6) show strong tuning already during the anticipation period, i.e. before the actual reward delivery. The inset in Figure 5B shows the mean value of the component #6 during the anticipation period, separating correct (green) and incorrect (red) trials for each stimulus. The characteristic U-shape for the error trials and the inverted U-shape for the correct trials agrees well with the predicted value of the rat’s uncertainty in each condition (Kepecs et al., 2008). Accordingly, this component can be interpreted as corresponding to the rat’s uncertainty or confidence about its own choice, confirming the results of Kepecs et al. (2008). In summary, both the main features of this dataset, as well as some of the subtleties that have been pointed out before (Kepecs et al., 2008) are picked up and reproduced by dPCA.
All dPCA components in each of the data sets are distributed across the whole neural population. For each component and each neuron, the corresponding encoder weight shows how much this particular component is exhibited by this particular neuron. For each component, the distribution of weights is strongly unimodal and centred at zero (Figure 6). In other words, there are no distinct sub-populations of neurons predominantly expressing a particular component; rather, each individual neuron can be visualized as a random linear combination of these components. We confirmed this observation by applying a recently developed clustering algorithm (Rodriguez and Laio, 2014)
to the population of neurons in the 15-dimensional space of dPC weights. In all cases, the algorithm found only one cluster (Figure S7; see Figure S8 for an analysis with Gaussian mixture models).
The complexity of responses in higher-order areas such as the prefrontal or orbitofrontal cortices has plagued researchers for quite some time. Here we have introduced a new dimensionality reduction method (demixed PCA) that is specifically targeted to resolve these difficulties. Our method summarizes the data by extracting a set of latent components from the single-cell activities. The critical advantage of dPCA, compared to standard methods such as PCA or factor analysis (FA), is that individual components do not show mixed selectivity, but are instead “demixed”. This demixing can greatly simplify exploration and interpretation of neural data. Indeed, in all cases presented here, all the major aspects of the population activity that had previously been reported are directly visible on the dPCA summary figure.
In the following discussion we first compare dPCA with various alternative approaches to the analysis of high-dimensional neural datasets, and then recapitulate what dPCA teaches us about neural activity in higher-order areas.
The most conventional and wide-spread method of analysis is to count the number of neurons that show significant response to a particular parameter of interest in a particular time period, and to then report the resulting percentage. This method was used in all of the original publications whose data we re-analyzed here: Romo et al. (1999) showed that 65% of the neurons (from those showing any task-related responses) exhibited monotonic tuning to the stimulation frequency during the delay period. Meyer et al. (2011) reported that the percentage of neurons with increased firing rate in the delay period increased after training from 15–21% to 27%, and the percentage of neurons tuned to match/non-match difference increased from 11% to 21%. Feierstein et al. (2006) found that 29% of the neurons were significantly tuned to reward, and 41% to the rat’s decision. Similarly, Kepecs et al. (2008) found that 42% of the neurons were significantly tuned to reward and observed that in the anticipation period these neurons demonstrated tuning to the expected uncertainty level. Note that in all of these cases our conclusions based on the dPCA analysis qualitatively agree with these previous findings.
Even though the conventional approach is perfectly sound, it does have several important limitations that dPCA is free of:
The conventional analysis focuses on a particular time bin (or on the average firing rate over a particular time period) and therefore provides no information about the time course of neural tuning. In contrast, dPCA results in time-dependent components, which highlight the time course of neural tuning.
If a parameter has more than two possible values, such as the vibratory frequencies in the somatosensory working memory task, then reporting a percentage of neurons with significant tuning disregards the shape of the tuning curve. On the other hand, a vertical slice through any stimulus-dependent dPCA component results in a tuning curve. In the case of “rainbow-like” stimulus components (Figure 2B, 3B, or 5B), these tuning curves are easy to imagine.
To count significantly tuned neurons, one chooses an arbitrary p-value cutoff (e.g. ). This can potentially create a false impression that neurons are separated into two distinct sub-populations: those tuned to a given parameter (e.g. stimulus) and those that are not tuned. This, however, is not the case in any of the datasets considered here: each demixed component is expressed in the whole population, albeit with varying strength (Figure 6).
If neural tuning to several parameters is analyzed separately (e.g. with several t-tests or several one-way ANOVAs instead of one multi-way ANOVA), then the results of the tests can get confounded due to different number of trials in different conditions (“unbalanced” experimental design). In case of dPCA, all parameters are analyzed together, and the issue of confounding does not arise.
An extended version of this approach uses multi-way ANOVA to test for firing rate modulation due to several parameters of interest, and repeats the test in a sliding time window to obtain the time course of neural tuning (see e.g. Sul et al., 2010, 2011). This avoids the first and the fourth limitations listed above, but the other two limitations still remain.
In many studies, the time-varying firing rates, or PSTHs, of individual neurons are simply averaged over a subset of the neurons that were found to be significantly tuned to a particular parameter of interest (e.g. Rainer et al., 1999; Kepecs et al., 2008). While this approach can highlight some of the dynamics of the population response, it ignores the full complexity and heterogeneity of the data (Figure S1), which is largely averaged out. The approach thereby fails to faithfully represent the neural activities, and may lead to the false impression pointed out in the third issue above, namely that there are distinct subpopulations of neurons with distinct response patterns.
Mante et al. (2013)
, have recently introduced a demixing approach based on multiple regression. The authors first performed a linear regression of neural firing rates to several parameters of interest, and then took the vectors of regression coefficients as demixing axes. While this method proved sufficient for the purposes of that study, it does not aim to faithfully capture all of the data, which leads to several disadvantages. Specifically, the approach (a) ignores the condition-independent components, (b) assumes that all neural tuning is linear, (c) finds only a single demixed component for each parameter, and (d) cannot achieve demixing when the axes for different parameters are far from orthogonal. To illustrate these disadvantages, we applied this method to all our datasets and show the results in Figure S9.
A significant advantage of the approach by Mante et al. (2013)
, is that it can deal with missing data or continuous task parameters, which dPCA currently can not. Future extensions of dPCA may therefore consider replacing the non-parametric dependencies of firing rates on task parameters with a parametric model, which would combine the advantages of both methods.
Another multivariate approach relies on linear classifiers that predict the parameter of interest from the population firing rates. The strength of neural tuning can then be reported as the cross-validated classification accuracy. For example,
Meyers et al. (2012), built linear classifiers for stimulus and match/non-match condition for the visuospatial working memory task analyzed above. Separate classifiers were used in each time bin, resulting in a time-dependent classification accuracy. The shape of the accuracy curve for the match/non-match condition follows the combined shapes of our decision components (and the same is true for stimulus). While the time-dependent classification accuracy provides an important and easily understandable summary of the population activity, it is far removed from the firing rates of individual neurons and is not directly representative of the neural tuning. The dPCA approach is more direct.Linear classifiers can also be used to inform the dimensionality reduction. For instance, linear discriminant analysis (LDA) reduces the dimensionality of a dataset while taking class labels into account. Whereas PCA looks for linear projections that maximize total variance (and ignores class labels), LDA looks for linear projections that maximize class separation, i.e. with maximal between-class and minimal within-class variance. Consequently, LDA is related to dPCA. However, LDA is a one-way technique (i.e. it considers only one parameter) and is not concerned with reconstruction of the original data, which makes it ill-suited for our purposes. See Supplementary Materials for an extended discussion on the differences between dPCA and LDA.
One interesting outcome of the analysis concerns the strength of the condition-independent components. These components capture 70–90% of the total, task-locked variance of neural firing rates. They are likely explained in part by an overall firing rate increase during certain task periods (e.g. during stimulus presentation). More speculatively, they could also be influenced by residual sensory or motor variables that vary rhythmically with the task, but are not controlled or monitored (Renart and Machens, 2014). The attentional or motivational state of animals, for instance, often correlates with breathing (Huijbers et al., 2014), pupil dilation (Eldar et al., 2013), body movements (Gouvêa et al., 2014), etc.
The second important observation is that parameter tuning “moves” during the trial from one component to another. To visualize these movements, we set up the dPCA algorithm such that it preferentially returns projections occupying only localized periods in time (see Methods). As a result, there are e.g. three separate stimulus components on Figure 2: one is active during the S1 period, one during the delay period, and one during the S2 period. The same can be observed in all other data sets as well: consider stimulus components in Figure 3 or decision components in Figures 4–5. The possibility to separate such components means that a neural subspace carrying information about a particular task parameter, changes (rotates) with time during the trial.
Our third finding is that most encoding axes turn out to be almost orthogonal to each other (only 22 out of 420 pairs are marked with stars on Figures 2–5E), even though dPCA, unlike PCA, does not enforce orthogonality. Non-orthogonality between two axes means that neurons expressing one of the components tend also to express the other one. Most examples of non-orthogonal pairs fall into the following three categories:
Non-orthogonality between a condition-dependent component and a condition-independent one (9 pairs), e.g. the first interaction (#4) and the second condition-independent (#2) components in Figure 4. This means that the neurons that are tuned to the presence/absence of reward will tend to have a specific time course of firing rate modulation, corresponding to component #2. Note that the same effect is observed in Figure 5.
Non-orthogonality between components describing one parameter (10 pairs), e.g. stimulus components in Figure 2. As many neurons express all three components, their axes end up being strongly non-orthogonal.
A special example is given by Figure 3, where the first stimulus component (#7) and the first interaction component (#10) are strongly non-orthogonal. As the interaction components in that particular data set are actually tuned to stimulus S2, this example can be seen a special case of the previous category.
Only two pairs do not fall into any of the categories listed above (first stimulus and first decision components in both olfactory datasets). In all other cases, the condition-dependent components turn out to be almost orthogonal to each other. In other words, task parameters are represented independently of each other. Indeed, if latent components are independently, i.e. randomly mapped to a space of neurons, then the encoding axes will be nearly orthogonal to each other, because in a high-dimensional space any two random vectors are close to being orthogonal (unlike 2D or 3D cases). Substantial deviations from orthogonality can only occur if different parameters are encoded in a population in a non-independent way. Our results indicate that this is mostly not the case. In addition to that, orthogonal readouts are arguably the most effective, as they maximize signal-to-noise ratio, ensuring successful demixing in downstream brain areas.
One limitation of dPCA is that it works only with discrete parameters (and all possible parameter combinations must be present in the data) and would need to be extended to be able to treat continuous parameters as well. Another is that the number of neurons needs to be high enough: we found that at least 100 neurons are usually needed to achieve satisfactory demixing in the data considered here. Furthermore, here we worked with trial-averaged PSTHs and ignored trial-to-trial variability. The datasets presented here could not have been analyzed differently, because the neurons were recorded across multiple sessions, and noise correlations between neurons were therefore not known. Demixed PCA could in principle also be used for simultaneous recordings of sufficient dimensionality. However, dPCA does not specifically address the question of how to properly treat trial-to-trial variability, and this problem may require further extensions in the future.
Here we argued that dPCA is an exploratory data analysis technique that is well suited for neural data, as it provides a succinct and immediately accessible overview of the population activity. It is important to stress its exploratory nature: dPCA enables a researcher to have a look at the full data, ask further questions and perform follow-up analyses and statistical tests. It should therefore be a beginning, not the end of the analysis.
The dPCA code for Matlab and Python is available at https://github.com/wielandbrendel/dPCA.
Brief descriptions of experimental paradigms are provided in the Results section and readers are referred to the original publications for all further details. Here we describe the selection of animals, sessions, and trials for the present manuscript. In all experiments neural recordings were obtained in multiple sessions, so most of the neurons were not recorded simultaneously.
Somatosensory working memory task in monkeys (Romo et al., 1999; Brody et al., 2003). We used the data from two monkeys (code names RR14 and RR15) that were trained with the same frequency set, and selected only the sessions where all six frequencies {10, 14, 18, 26, 30, 34} Hz were used for the first stimulation (other sessions were excluded). Monkeys made few mistakes (overall error rate was 6%), and here we analyzed only correct trials. Monkey RR15 had additional 3 s delay after the end of the second stimulation before it was cued to provide the response. Using the data from monkey RR13 (that experienced different frequency set) led to very similar dPCA components (data not shown).
Visuospatial working memory task in monkeys (Qi et al., 2011; Meyer et al., 2011; Qi et al., 2012). We used the data from two monkeys (code names AD and EL) that were trained with the same spatial task. Monkeys made few mistakes (overall error rate was 8%), and here we analysed only correct trials. The first visual stimulus was presented at 9 possible spatial locations arranged in a 33 grid (Figure 3A); here we excluded all the trials where the first stimulus was presented in the centre position.
Olfactory discrimination task in rats (Feierstein et al., 2006). We used the data from all five rats (code names N1, P9, P5, T5, and W1). Some rats were trained with two distinct odours, some with four, some with six, and one rat experienced mixtures of two fixed odours in varying proportions. In all cases each odour was uniquely associated with one of the two available water ports (left/right). Following the original analysis (Feierstein et al., 2006), we grouped all odours associated with the left/right reward together as a “left/right odour”. For most rats, caproic acid and 1-hexanol (shown on Figures 4–5A) were used as the left/right odour. We excluded from the analysis all trials that were aborted by rats before reward delivery (or before waiting 0.2 s at the reward port for the error trials).
Olfactory categorization task in rats (Kepecs et al., 2008). We used the data from all three rats (code names N1, N48, and N49). Note that recordings from one of the rats (N1) were included in both this and previous datasets; when we excluded if from either of the datasets, the results stayed qualitatively the same (data not shown). We excluded from the analysis all trials that were aborted by rats before reward delivery (or before waiting 0.3 s at the reward port for the error trials).
In each of the datasets each trial can be labeled with two parameters: “stimulus” and “decision”. Note that a “reward” label is not needed, because its value can be deduced from the other two due to the deterministic reward protocols in all tasks. For our analysis it is important to have recordings of each neuron in each possible condition (combination of parameters). Additionally, we required that in each condition there were at least
trials, to reduce the standard error of the mean when averaging over trials, and also for cross-validation purposes. The cutoff was set to
for both working memory datasets, and to for both olfactory datasets (due to less neurons available).We have further excluded very few neurons with mean firing rate over 50 Hz, as neurons with very high firing rate can bias the variance-based analysis. Firing rates above 50 Hz were atypical in all datasets (number of excluded neurons for each dataset: 5 / 2 / 1 / 0). This exclusion had a minor positive effect on the components. We did not apply any variance-stabilizing transformations, but if the square-root transformation was applied, the results stayed qualitatively the same (data not shown).
No other pre-selection of neurons was used. This procedure left 832 neurons (230 / 602 for individual animals, order as above) in the somatosensory working memory dataset, 956 neurons (182 / 774) in the visuospatial working memory dataset, 437 neurons in the olfactory discrimination dataset (166 / 30 / 9 / 106 / 126), and 214 neurons in the olfactory categorization dataset (67 / 38 / 109).
Spike trains were filtered with a Gaussian kernel ( ms) and averaged over all trials in each condition to obtain smoothed average peri-stimulus time histograms (PSTHs) for each neuron in each condition.
In the visuospatial working memory dataset we identified the preferred direction of each neuron as the location that evoked maximum mean firing rate in the 500 ms time period while the first stimulus was shown. The directional tuning was shown before to have a symmetric bell shape (Qi et al., 2011; Meyer et al., 2011), with each neuron having its own preferred direction. We then re-sorted the trials (separately for each neuron) such that only 5 distinct stimuli were left: preferred direction, 45, 90, 135, and 180 away from the preferred direction.
In both olfactory datasets trials were self-paced. This means that trials could have very different length, making averaging of firing rates over trials impossible. We used the following re-stretching procedure (separately in each dataset) to equalize the length of all trials and to align several events of interest (see Figure S10 for the illustration). We defined five alignment events: odour poke in, odour poke out, water poke in, reward delivery, and water poke out. First, we aligned all trials on odour poke in () and computed median times of four other events (for the time of reward delivery we took the median over all correct trials). Second, we set to be the minimal waiting time between water port entry and reward delivery across the whole experiment ( s for the olfactory discrimination task and s for the olfactory categorization task). Finally, for each trial we computed the PSTH as described above, set , to be the times of alignment events on this particular trial (for error trials we took ), and stretched along the time axis in a piecewise-linear manner to align each with the corresponding .
We made sure that this procedure did not introduce any artifacts by considering an alternative procedure, where short (450ms) time intervals around each were cut out of each trial and concatenated together; this procedure is similar to the pooling of neural data performed in the original studies (Feierstein et al., 2006; Kepecs et al., 2008). The dPCA analysis revealed qualitatively similar components (Figure S6).
For each neuron (out of ), stimulus (out of ), and decision (out of ) we have a collection of trials. For each trial we have a filtered spike train sampled at time points. The number of trials can vary with , , and , and so we denote it by . For each neuron, stimulus, and decision we average over these trials to compute the mean firing rate . These data can be thought of as time-dependent neural trajectories (one for each condition) in a -dimensional space (Figure 1A). The number of distinct data points in this -dimensional space is . We collect them in one matrix with dimensions (i.e. rows, columns).
Let be centred, i.e. the mean response of any single neuron over all times and conditions is zero. As we previously showed in (Brendel et al., 2011), can be decomposed into independent parts that we call “marginalizations” (Figure S11):
Here denotes the time-varying, but stimulus- and decision-invariant part of , which can be obtained by averaging over all stimuli and decisions (i.e. over all conditions), i.e. (where the angle brackets denote the average over the subscripted parameters). The stimulus-dependent part is an average over all decisions of the part that is not explained by . Similarly, the decision-dependent part is an average of the same remaining part over all stimuli. Finally, the higher-order, or “interaction”, term is calculated by subtracting all simpler terms, i.e. . As a matrix, each marginalization has exactly the same dimensions as (e.g. is not a -dimensional matrix, but a matrix with identical values for every time point).
In this manuscript we are not interested in separating neural activity that varies only with stimulus (but not with time) from neural activity that varies due to interaction of stimulus with time. More generally, however, we can treat all parameters on an equal footing. In this case, data labeled by parameters can be marginalized into marginalization terms. For the most general decomposition is given by
See Supplementary Materials for more details.
We have additionally separated stimulus, decision, and interaction components into those having variance in different time periods of the trial. Consider the somatosensory working memory task (Figure 2). Each trial can be reasonably split into three distinct periods: up until the end of S1 stimulation, delay period, and after beginning of S2 stimulation. We aimed at separating the components into those having variance in only one of these periods. For this, we further split e.g. into , , such that and each of the parts equals zero outside of the corresponding time interval. This was done with stimulus, decision, and interaction marginalizations (but not with the condition-independent one).
As a result, three stimulus components shown on Figure 2B have very distinct time course. Note, however, that the angles between corresponding projection axes are far from orthogonal (Figure 2E). This means that most of the neurons expressing one of these components express other two as well, and so if the splitting had not been enforced, these components would largely have joined together (Figure S11).
For the visuospatial working memory dataset we split trials into four parts: before the end of S1, delay period, S2 period, second delay period. For both olfactory datasets we split trials into two parts on water poke in.
Importantly, we made sure that this procedure did not lead to any noticeable loss of explained variance, as compared to dPCA without time period separation (Figure S12). On the other hand, it often makes individual components easier to interpret (e.g. second stimulus components on Figures 2 and 3 are only active during the delay periods and so are clear “memory” components), and highlights the fact that parameter (e.g. stimulus) representations tend to rotate with time in the firing rate space.
Given a decomposition , dPCA aims at finding directions in that capture as much variance as possible, with an additional restriction that the variance in each direction should come from only one of . Let us assume that we want the algorithm to find directions for each marginalization (i.e. directions in total). The cost function for dPCA is
where is an encoder matrix with columns and is a decoder matrix with rows. Here and below matrix norm signifies Frobenius norm, i.e. . Without loss of generality we assume that has orthonormal columns and that components are ordered such that their variance (row variance of ) is decreasing. The term regularizes the solution to avoid overfitting (Figure S13).
For the optimization problem can be understood as a reduced rank regression problem (Izenman, 1975; Reinsel and Velu, 1998): minimize with and . The minimum can be found in three steps:
Compute the standard regression solution to the unconstrained problem of minimizing ;
Project on the -dimensional subspace with highest variance to incorporate the rank constraint, i.e. , where is a matrix of leading singular vectors of ;
Factorize with and to recover the encoder and decoder.
Finally, notice that the regularized problem can be reduced to the unregularized case by replacing and where and are zero and unit matrices. See Supplementary Materials for a more detailed mathematical treatment.
We found that a very good approximation to the optimal solution can be achieved in a simpler way that can provide some further intuition. Let be equal to a matrix of the first principal axes (singular vectors) of , join all together horizontally to form one matrix , and take , pseudo-inverse of . This works well, provided that is chosen to capture most of the signal variance of , but not larger. Choosing too small results in poor demixing, and choosing too large results in overfitting. We found that provides a good trade-off in all datasets considered here. However, the general method described above is a preferred approach, as it does not depend on the choice of , provides a more accurate regularization and is derived from an explicit objective (hence avoids hidden assumptions).
This approximate solution highlights the conditions under which dPCA will work, i.e. will result in well-demixed components (components having variance in only one marginalization) that together capture most of the variance of the data. For this to work, the main principal axes of different marginalizations should be non-collinear, or in other words, principal subspaces of different marginalizations should not overlap.
We used cross-validation to find the optimal regularization parameter for each dataset. To separate the data into training and testing sets, we held out one random trial for each neuron in each condition as a set of test “pseudo-trials” (as the neurons were not recorded simultaneously, we do not have recordings of all neurons in any actual trial). Remaining trials were averaged to form a training set . Note that and have the same dimensions. We then performed dPCA on for various values of , selected 10 components in each marginalization (i.e. 40 components in total) to obtain and , and computed
as a fraction of variance not explained by the held-out data. We repeated this procedure 10 times for different train-test splittings and averaged the resulting functions . In all cases the average function had a clear minimum (Figure S14) that we selected as the optimal . An alternative formula for using only yielded the same results, see Supplementary Information.
The values used for each dataset were / / / times the total variance of the corresponding dataset.
It can be argued that only the decoding axes are of interest (and not the encoding axes). In the toy example shown on Figure 1F, decoding axes roughly correspond to discriminant axes that linear discriminant analysis (LDA) would find when trying to decode stimulus and time. This remains true in the case of real data as well (see Supplementary Materials). However, without encoding axes there is no way to reconstruct the original dataset and therefore no way in assigning “explained variance” to principal components.
On the other hand, it can be argued that only the encoding axes are of interest. Indeed, in the same toy example shown on Figure 1F, encoding axes roughly correspond to principal axes of the stimulus and time marginalizations. In other words, they show directions along which stimulus and time are varying the most. Correspondingly, instead of performing dPCA, one can perform standard PCA in each marginalization and analyze resulting components inside each marginalization (they will, by definition, be perfectly demixed). However, this makes a transformation data dPC complex and involving a series of multi-trial averages. The brain has arguably to rely on standard linear projections to do any kind of trial-by-trial inference, and so to gain insight into the utility of the population code for a biological system, we prefer the method involving only linear projections of the full data. This is achieved with a decoder.
We believe therefore that both encoder and decoder are essential for our method and treat them on equal footing, in line with the schematic shown on Figure 1E.
The marginalization procedure ensures that the total covariance matrix is equal to a sum of covariance matrices from each marginalization:
This means that the the variance of in each direction can be decomposed into a sum of variances due to to time, stimulus, decision, and stimulus-decision interaction. This fact was used to produce bar plots shown on Figures 2–5D. Namely, consider a dPC with a decoding vector (that does not necessarily have a unit length). Total variance of this dPC is and it is equal to the sum of marginalized variances .
This allows us to define a “demixing index” of each component as . This index can range from 0.25 to 1, and the closer it is to 1, the better demixed the component is. As an example, for the somatosensory working memory dataset, the average demixing index over the first 15 PCA components is 0.760.16 (meanSD), and over the first 15 dPCA components is 0.970.02, which means that dPCA achieves much better demixing (, Mann-Whitney-Wilcoxon ranksum test). For comparison, the average demixing index of individual neurons in this dataset is 0.550.18. In other datasets these numbers are similar.
To compute the cumulative variance explained by the first components, we cannot simply add up individual variances, as the demixing axes are not orthogonal. Instead, on Figures 2–5C we show “explained variance”, computed as follows: Let the decoding vectors for these axes be stacked as rows in a matrix and encoding vectors as columns in a matrix . Then the proportion of total explained variance is given by
Note that for the standard PCA when , the standard explained variance (sum of the first eigenvalues of the covariance matrix over the sum of all eigenvalues) can be given by an analogous formula:
Following (Machens et al., 2010), in Figures 2–5C we applied a correction to show the amount of explained “signal variance”. Assuming that each PSTH consists of some trial-independent signal and some random noise , , the average PSTH will be equal to . If the number of trials
is not very large, some of the resulting variance will be due to the noise term. To estimate this variance for each neuron in each condition, we took two random trials
and and considered . These data form a data matrix of the same dimensions as , which has no signal but approximately the same amount of noise as . The following text assumes that is centred.Singular values of provide an approximate upper bound of the amount of noise variance in each successive PCA or dPCA component. Therefore, the cumulative amount of signal variance for PCA is given by
where and are singular values of and respectively. For dPCA, the formula becomes
Pie charts in Figures 2–5D show the amount of signal variance in each marginalization. To compute it, we marginalize , obtain a set of , and then compute signal variance in marginalization as . The sum over all marginalizations is equal to the total amount of signal variance .
On Figures 2–5E stars mark the pairs of components whose encoding axes and are significantly and robustly non-orthogonal. These were identified as follows: In Euclidean space of
dimensions, two random unit vectors (from a uniform distribution on the unit sphere
) have dot product (cosine of the angle between them) distributed with mean zero and standard deviation
. For large the distribution is approximately Gaussian. Therefore, if , we say that the axes are significantly non-orthogonal ().Coordinates of quantify how much this component contributes to the activity of each neuron. Hence, if cells exhibiting one component also tend to exhibit another, the dot product between the axes is positive (note that is approximately equal to the correlation between the coordinates of and ). Sometimes, however, the dot product has large absolute value only due to several outlying cells. To ease interpretation, we marked with stars only those pairs of axes for which the absolute value of Spearman (robust) correlation was above 0.2 with (in addition to the above criterion on ).
We used decoding axis of each dPC in stimulus, decision, and interaction marginalizations as a linear classifier to decode stimulus, decision, or condition respectively. Black lines on Figures 2–5B show time periods of significant classification. A more detailed description follows below.
We used 100 iterations of stratified Monte Carlo leave-group-out cross-validation, where on each iteration we held out one trial for each neuron in each condition as a set of test “pseudo-trials” and averaged over remaining trials to form a training set (see above). After running dPCA on , we used decoding axes of the first three stimulus/decision/interaction dPCs as a linear classifier to decode stimulus/decision/condition respectively. Consider e.g. the first stimulus dPC: first, for each stimulus, we computed the mean value of this dPC separately for every time-point. Then we projected each test trial on the corresponding decoding axis and classified it at each time-point according to the closest class mean. Proportion of test trials (out of ) classified correctly resulted in a time-dependent classification accuracy, which we averaged over 100 cross-validation iterations. Note that this is a stratified procedure: even though in reality some conditions have much fewer trials than others, here we classify exactly the same number of “pseudo-trials” per condition. At the same time, as the coordinates of individual data points in each pseudo-trial are pooled from different sessions, the influence of noise correlations on the classification accuracies is neglected, similar to Meyers et al. (2012).
We then used 100 shuffles to compute Monte Carlo distribution of classification accuracies expected by chance. On each iteration for each neuron we shuffled all available trials between conditions, respecting the number of trials per condition (i.e. all trials were shuffled and then randomly assigned to the conditions such that all values stayed the same). Then exactly the same classification procedure as above (with 100 cross-validation iterations) was applied to the shuffled dataset to find mean classification accuracy for the first stimulus, decision, and interaction components. All 100 shuffling iterations resulted in a set of 100 time-dependent accuracies expected by chance.
The time periods when actual classification accuracy exceeded all 100 shuffled decoding accuracies in at least 10 consecutive time bins are marked by black lines on Figures 2–5. Components without any periods of significant classification are not shown. See Figure S4 for classification accuracies in each dataset. Monte Carlo computations took 24 hours for each of the larger datasets on a 6 core 3.2 Ghz Intel i7-3930K processor.
Journal of Multivariate Analysis
, 5(2):248–264.
Comments
There are no comments yet.