Recordings from higher-order areas in awake behaving animals often yield a large variety of neural responses (see e.g., Miller, 1999; Churchland and Shenoy, 2007; Jun et al., 2010; Machens et al., 2010). These observations at the level of individual cells could imply a complicated and intricate response at the population level for which a simplified description does not exist. Alternatively, the large heterogeneity of responses may be the result of a simple mixing procedure. For instance, response variety can come about if the responses of individual neurons are random, linear mixtures of a few generic response components (see e.g., Eliasmith and Anderson, 2003).

To illustrate this insight, we will construct a simple toy model. Imagine an animal which performs a two-alternative-forced choice task (Newsome et al., 1989; Uchida and Mainen, 2003). In each trial of such a task, the animal receives a sensory stimulus, s, and then makes a binary decision, d, based on whether s falls into one of two response categories. If the animal decides correctly, it receives a reward. We will assume that the activity of the neurons in our toy model depends only on the stimulus s and the decision d.

To obtain response heterogeneity, we construct the response of each neuron as a random, linear mixture of two underlying response components, one that represents the stimulus, z₁(t,s), and one that represents the decision, z₂(t,d), see Figure 1A. The time-varying firing rate of neuron i is then given by

Here, the parameters a_i1 and a_i2 are the mixing coefficients of the neuron, the bias parameter c_i describes a constant offset, and the term η_i(t) denotes additive, white noise. We assume that the noise of different neurons can be correlated so that

where the angular brackets denote averaging over time, and H_ij is the noise covariance between neuron i and j. We will assume that there are N neurons and, for notational compactness, we will assemble their activities into one large vector, r(t,s,d) = (r₁(t,s,d),…,r_N(t,s,d))^T. After doing the same for the mixing coefficients, the constant offset, and the noise, we can write equivalently,

Without loss of generality, we can furthermore assume that the mixing coefficients are normalized so that aiTai=1 for i∈{1,2}. Since we assume that the mixing coefficients are drawn at random, and independently of each other, the first and second coefficient will be uncorrelated, so that on average, a1Ta2=0, implying that a₁ and a₂ are approximately orthogonal.

With this formulation, individual neural responses mix information about the stimulus s and the decision d, leading to a variety of responses, as shown in Figure 1B. While with only two underlying components, the overall heterogeneity of responses remains limited, the response heterogeneity increases strongly when more components are allowed (see Figures 3A,B for an example with three components).

The standard approach to deal with such data sets is to sort cells into categories. In our example, this approach may yield two overlapping categories of cells, one for cells that respond to the stimulus and one for cells that respond to the decision. While this approach tracks down which variables are represented in the population, it will fail to quantify the exact nature of the population activity, such as the precise co-evolution of the neural population activity over time.

A common approach to address these types of problems are dimensionality reduction methods such as PCA (Nicolelis et al., 1995; Friedrich and Laurent, 2001; Hastie et al., 2001; Zacksenhouse and Nemets, 2008; Machens et al., 2010). The main aim of PCA is to find a new coordinate system in which the data can be represented in a more succinct and compact fashion. In our toy example, even though we may have many neurons with different responses (N = 50 in Figure 1, with five examples shown in Figure 1B), the activity of each neuron can be represented by a linear combination of only two components. In the N-dimensional space of neural activities, the two components, z₁(t,s) and z₂(t,d), can be viewed as two coordinates of a coordinate system whose axes are given by the vectors of mixing coefficients, a₁ and a₂. Since the first two coordinates capture all the relevant information, the components live in a two-dimensional subspace. Using PCA, we can retrieve the two-dimensional subspace from the data. While the method allows us to reduce the dimensionality and complexity of the data dramatically, PCA will in general only retrieve the two-dimensional subspace, but not the original coordinates, z₁(t,s) and z₂(t,d).

To see this, we will briefly review PCA and show what it does to the data from our toy model. PCA commences by computing the covariances of the firing rates between all pairwise combination of neurons. Let us define the mean firing rate of neuron i as the average number of spikes that this neuron emits, so that

We will use the angular brackets in the second line as a short-hand for averaging. The variables to be averaged over are indicated as subscript on the right bracket. Here, the average runs over all time points t, all stimuli s, and all decisions d. For the vector of mean firing rates we write r = (r₁,…,r_N)^T.

The covariance matrix of the data summarizes the second-order statistics of the data set,

and has size N × N where N is the number of neurons in the data set. Given the covariance matrix, we can compute the firing rate variance that falls along arbitrary directions in state space. For instance, the variance captured by a coordinate axis given by a normalized vector u is simply L = u^TCu. We can then look for the axis that captures most of the variance of the data by maximizing the function L with respect to u subject to the normalization constraint u^Tu = 1. The solution corresponds to the first axis of the coordinate system that PCA constructs. If we are looking for several mutually orthogonal axes, these can be conveniently summarized into an N × n orthogonal matrix, U = [u₁,…,u_n]. To find the maximum amount of variance that falls into the subspace spanned by these axes, we need to maximize

where the trace-operation, tr(·), sums over all the diagonal entries of a matrix, and I_n denotes the n × n identity matrix.

Mathematically, the principal axes u_i correspond to the eigenvectors of the covariance matrix, C, which can nowadays be computed quite easily using numerical methods. Subsequently, the data can be plotted in the new coordinate system. The new coordinates of the data are given by

These new coordinates are called the principal components. Note that the new coordinate system has a different origin from the old one, since we subtracted the vector of mean firing rates, r. Consequently, the principal components can take both negative and positive values. Note also that the principal components are only defined up to a minus sign since every coordinate axis can be reflected along the origin. For our artificial data set, only two eigenvalues are non-zero, so that two principal components suffice to capture the complete variance of the data. The data in these two new coordinates, y₁(t,s,d) and y₂(t,s,d), are shown in Figure 1C.

Our toy model shows how PCA can succeed in summarizing the population response, yet it also illustrates the key problem of PCA: just as the individual neurons, the components mix information about the different task parameters (Figure 1C), even though the original components do not (Figure 1A). The underlying problem is that PCA ignores the causes of firing rate variability. Whether firing rates have changed due to the external stimulus s, due to the internally generated decision d, or due to some other cause, they will enter equally into the computation of the covariance matrix and therefore not influence the choice of the coordinate system constructed by PCA.

To make these notions more precise, we compute the covariance matrix of the simulated data. Inserting Eq. 3 into Eq. 6, we obtain

where M₁₁ and M₂₂ denote firing rate variance due to the first and second component, respectively, M₁₂ denotes firing rate variance due to a mix of the two components, and H is the covariance matrix of the noise. Using the short-hand notations z₁(t) = 〈z₁(t,s)〉_s, z₂(t) = 〈z₂(t,d)〉_d, and z_i = 〈z_i(t)〉_t for i∈[1,2], the different variances are given by

Principal component analysis will only be able to segregate the stimulus- and decision-dependent variance if the mixture term M₁₂ vanishes and if the variances of the individual components, M₁₁ and M₂₂, are sufficiently different from each other. However, if the two underlying components z₁(t,s) and z₂(t,d) are temporally correlated, then the mixture term M₁₂ will be non-zero. Its presence will then force the eigenvectors of C away from a₁ and a₂. Moreover, even if the mixture term vanishes, PCA may still not be able to retrieve the original mixture coefficients, if the variances of the individual components, M₁₁ and M₂₂ are too close to each other when compared to the magnitude of the noise: in this case the eigenvalue problem becomes degenerate. In general, the covariance matrix therefore mixes different origins of firing rate variance rather than separating them. While PCA allows us to reduce the dimensionality of the data, the coordinate system found may therefore provide only limited insight into how the different task parameters are represented in the neural activities.

To solve these problems, we need to separate the different causes of firing rate variability. In the context of our example, we can attribute changes in the firing rates to two separate sources, both of which contribute to the covariance in Eq. 6. First, firing rates may change due to the externally applied stimulus s. Second, firing rates may change due to the internally generated decision d.

To account for these separate sources of variance in the population response, we suggest to estimate one covariance matrix for every source of interest. Such a covariance matrix needs to be specifically targeted toward extracting the relevant source of firing rate variance without contamination by other sources. Naturally, this step is somewhat problem-specific. For our example, we will first focus on the problem of estimating firing rate variance caused by the stimulus separately from firing rate variance caused by the decision. When averaging over all stimuli, we obtain the marginalized firing rates r(t,d) = 〈r(t,s,d)〉_s. The covariance caused by the stimulus is then given by the N × N matrix

We will refer to C_s as the marginalized covariance matrix for the stimulus. We can repeat the procedure for the decision-part of the task. Marginalizing over decisions, we obtain r(t,s) = 〈r(t,s,d)〉_d and

Having two different covariance matrices, one may now perform two separate PCAs, one for each covariance matrix. In turn, one obtains two separate coordinate systems, one in which the principal axes point into the directions of state space along which firing rates vary if the stimulus is changed, the other in which they point into the directions along which firing rates vary if the decision changes.

For the toy model, it is readily seen that the marginalized covariance matrices are given by Cs=a1a1TMs,11+H and Cd=a2a2TMd,22+H with M_s,11 = 〈(z₁(t,s) − z₁(t))²〉 and M_d,22 = 〈(z₂(t,d) − z₂(t))²〉. Consequently, the principal eigenvectors of C_s and C_d will be equivalent to the mixing coefficients a₁ and a₂, at least as long as the variances M_s,11 and M_d,22 are much larger than the size of the noise, which is given by tr(H).

If the noise term is not negligible, it will force the eigenvectors away from the actual mixing coefficients. This problem can be alleviated by using the orthogonality condition, a1Ta2=0, which implies that there are separate sources of variance for the stimulus- and decision-components. To this end, we can seek to divide the full space into two subspaces, one that captures as much as possible about the stimulus-dependent covariance C_s, and another, that captures as much as possible about the decision-dependent covariance C_d. Our goal will then be to maximize the function

with respect to the two orthogonal matrices U₁ and U₂ whose columns contain the basis vectors of the respective subspaces. The first term in Eq. 15 captures the total variance falling into the subspace spanned by the columns of U₁, and the second term the total variance falling into the subspace given by U₂. Writing U = [U₁,U₂], we obtain an orthogonal matrix for the full space, and the orthogonality conditions are neatly summarized by UU^T = I. As shown in the Appendix, the maximization of Eq. 15 under these orthogonality constraints can be solved by computing the eigenvectors and eigenvalues of the difference of covariance matrices,

In this case, the eigenvectors belonging to the positive eigenvalues of D form the columns of U₁ and the eigenvectors belonging to the negative eigenvalues of D form the columns of U₂. As with PCA, the positive or negative eigenvalues can be sorted according to the amount of variance they capture about C_s and C_d.

For the simulated example, we obtain

where the noise term H has now dropped out. Diagonalization of D results in two clearly separated eigenvalues, M_s,11 and −M_d,11, and in two eigenvectors, a₁ and a₂, that correspond to the original mixing coefficients.

As a result of the above method, we obtain a new coordinate system, whose basis vectors are given by the columns of the matrix U. This coordinate system provides simply a different, and hopefully useful, way of representing the population response. One major advantage of orthogonality is that one can easily move back and forth between the single cell and population level description of the neural activities. Just as in PCA, we can project the original firing rates of the neurons onto the new coordinates,

and the two leading coordinates for the toy model are shown in Figure 1D. These components correspond approximately to the original components, z₁(t,s) and z₂(t,d). In turn, we can reconstruct the activity of each neuron by inverting the coordinate transform,

For every neuron this yields a set of N reconstruction coefficients which correspond to the rows of U.

Since two coordinates were sufficient to capture most of the variance in the toy example, the firing rate of every neuron can be reconstructed by a linear combination of these two components, y₁(t,s,d) and y₂(t,s,d). For each neuron, we thereby obtain two reconstruction coefficients, u_i1 and u_i2. The set of all reconstruction coefficients constitutes a cloud of points in a two-dimensional space. The distribution of this cloud, together with the activities of several example neurons are shown in Figure 2. This plot allows us to link the single cell with the population level by visualizing how the activity of each neuron is composed out of the two underlying components.

In our toy example, we have assumed that each task parameter is represented by a single component. We note that this is a feature of our specific example. In more realistic scenarios, a single task parameter could potentially be represented by more than one component. For instance, if one set of neurons fires transiently with respect to a stimulus s, but another set of neurons fires tonically, then the firing rate dynamics of the stimulus representation are already two-dimensional, even without taking the decision into account. In such a case, we can still use the method described above to retrieve the two subspaces in which the respective components lie.

However, the number of task parameters will often be larger than two. In the two-alternative-forced choice task, there are at least four parameters that could lead to changes in firing rates: the timing of the task, t, potentially related to anticipation or rhythmic aspects of a task, the stimulus, s, the decision, d, and the reward, r. Even more task parameters could be of interest, such as those extracted from previous trials etc.

These observations raise the question of how the method can be generalized if there are more than two task parameters to account for. To do so, we write the relevant parameters into one long vector θ = (θ₁,θ₂,…,θ_M), and assume that the firing rates of the neurons are linear mixtures of the form

where each task parameter is now represented by more than one component. For each parameter, θ_i, we can compute the marginalized covariance matrix,

which measures the covariance in the firing rates due to changes in the parameter θ_i. Diagonalizing each of these covariance matrices will retrieve the various subspaces corresponding to the different mixture coefficients. For instance, when diagonalizing C₁, we obtain the subspace for the components that depend on the parameter θ₁. The relevant eigenvectors of C₁ will therefore span the same subspace as the mixture coefficients a₁₁, a₁₂, etc., in Eq. 22.

As before, the method's performance under additive noise can be enhanced by maximizing a single function (see Appendix)

subject to the orthogonality constraint U^TU = I for U = [U₁,U₂,…,U_M]. Maximization of this function will force the firing rate variance due to different parameters θ_i into orthogonal subspaces (as required by the model). If M = 1, then maximization results in a standard PCA. In the case M = 2, maximization requires the diagonalization of the difference of covariance matrices C₁ − C₂, as in Eq. 16. In the case M > 2, various algorithms can be constructed to find local maxima of L (see e.g., Bolla et al., 1998). To our knowledge, a full understanding of the global solution structure of the maximization problem does not exist for M > 2. In the Appendix, we show how to maximize Eq. 24 with standard gradient ascent methods. In any case, it may often be a good idea to use PCA on the full covariance matrix of the data, Eq. 6, to reduce the dimensionality of the data set prior to the demixing procedure. Indeed, this preprocessing step was applied in Machens et al. (2010).

The above formulation of the problem may be further generalized by allowing individual components to mix parameters in non-trivial ways. To study this scenario in a simple example, imagine that in the above two-alternative-forced choice task, in addition to the stimulus- and decision-dependent component, there were a purely time-dependent component, z₃(t), locked to the time structure of the task, so that

This scenario is illustrated in Figures 3A,B. As before, we can compute marginalized covariance matrices, that capture the covariance due to the stimuli s, the decisions d, or the time points t. While the marginalized covariance matrices for the stimuli and decisions, C_s and C_d, have one significant eigenvalue each, and thereby capture the relevant component (Figure 3C), the marginalized covariance matrix for time, C_t, now has three significant eigenvalues, and therefore does not allow us to retrieve the purely time-dependent component z₃(t). The reason for this failure is that all three components in Eq. 25 have a time-dependence that cannot be averaged out. By design, the stimulus-averaged first component, z₁(t) = 〈z₁(t,s)〉_s, and the decision-averaged second component, z₂(t) = 〈z₂(t,d)〉_d do not vanish. In other words, the stimulus- and decision-components have intrinsic time-dependent variance that cannot be separated from the stimulus- or decision-induced variance.

Consequently, the subspace spanned by the first three eigenvectors of C_t overlaps with the respective subspaces spanned by the first eigenvectors of C_s and C_d. One way to visualize this overlap is to take the five relevant eigenvectors (three for C_t, one for C_s, and one for C_d) and compute how much of the variance of each marginalized covariance matrix they capture. To do so, we compute the “confusion matrix”

This confusion matrix measures what percentage of the variance attributed to the j-th cause is captured by the i-th coordinate. For the above example, it is illustrated in Figure 3D. If in one row of this matrix, more than one entry is significantly above 0, then more than one covariance matrix has significant variance along that direction of state space. Whereas the eigenvectors of the C_s and C_d matrix do not interfere with each other, i.e., they are approximately orthogonal, the eigenvectors of the C_t matrix interfere with both the C_s and C_d eigenvectors, i.e., the respective subspaces overlap. The method introduced above will still yield a result in this case, however, the new coordinate system will generally not retrieve the original components.

An ad hoc solution to this problem may be to section the three-dimensional eigenvector subspace of C_t, and identify a direction that is orthogonal to the first eigenvectors of C_s and C_d, which will then correspond to the purely time-dependent component z₃(t). Alternatively, we could restrict the estimation of C_t to the time before stimulus onset, so that the covariance matrix is no longer contaminated by time-dependent variance from the stimulus- or decision-components. The rank of C_t then reduces to one, and the different components separate nicely (Figures 3E,F,G). While feasible in our toy scenario, these ad hoc procedures are not guaranteed to work for real data, when more dimensions are involved, and more complex confusion matrices may result. However, the latter solution demonstrates that by a judicious choice of marginalized covariance matrices, one may sometimes be able to avoid such problems of non-separability.

In all of these scenarios, we assumed that the firing rates r are linear mixtures of a set of underlying sources z, each with mean 0, so that

The problem that we have been describing then consists in estimating the unknown sources, z, the unknown mixture coefficients, A, and the unknown bias parameters c from the observed data, r. Without loss of generality, we can assume that the sources are centered so that 〈z〉 = 0. Ours is therefore a specific version of the much-studied blind source separation problem (see e.g., Molgedey and Schuster, 1994; Bell and Sejnowski, 1995). In many standard formulations of this problem, one assumes that the sources are uncorrelated, or even statistically independent, which implies that the covariance matrix of the sources, M = 〈zz^T〉_t, is diagonal.

In our case, we do not want to make this assumption, which rules out the use of many blind source separation methods, such as independent component analysis (Hyvärinen et al., 2001). On the upside, we do have additional information, in the form of n task parameters, that provide indirect clues toward the underlying sources. More specifically, we assume that the sources are of the form z_k(t,θ_k) where θ_k denotes a single task parameter, or a specific combination of task parameters. For each task parameter, we can estimate the marginalized covariance matrix C_i, which in turn is given by C_i = AM_iA^T with

As long as different task parameters are distributed over different components, the matrix M_i will be block-diagonal. In the most general case, however, as discussed above, this will not be true. If one parameter is shared among several components, then the respective marginalized covariance matrix will capture variance from all of these components, and maximization of Eq. 24 will not necessarily retrieve the original components. Future work may show how this general, semi-blind source separation problem can be solved by using knowledge about the structure of the marginalized M-matrices. For now, we suggest that in many practical scenarios, a judicious choice of covariance measurements, for instance, by focusing on particular time intervals of a task etc., may help to partly reduce the problem to those that are completely separable, as in Eq. 22.