Data structure<\/p>\n
We used the International Brain Laboratory (IBL) public data release4<\/a>. For each experimental session, we collected time-series data on task, behaviour and electrophysiological recordings. These were segmented into trials based on key task events. The task recordings collected for each trial included information on the block prior as well as stimulus contrast and location. The behaviour recordings for each trial comprised the choice made, the outcome\/reward received and the time-varying movement such as wheel movement velocity, whisker motion energy and licks. Other behavioural variables, such as paw movement, body motion energy and pupil diameter traces, could be potentially included. However, we did not include them because many sessions had missing values. The electrophysiological recordings for each trial contained time-varying spike trains of recorded neurons. All of these recordings can be accessed directly through the IBL\u2019s open API. The following section describes the steps for preprocessing these raw data into data matrices for the encoding model.<\/p>\n Criteria for session inclusion<\/p>\n We iterated over all cortical regions and downloaded the related sessions. Sessions were included only if all of the behaviour recordings (wheel velocity, whisker motion energy and licks) and electrophysiological data were in place. A maximum of 30 sessions was included per cortical area to encourage a more balanced coverage. Some analyses required additional inclusion criteria, such as a minimum number of trials per condition. These analysis-specific criteria are discussed in the relevant sections below.<\/p>\n Criteria for trial inclusion<\/p>\n All trials from the left or right unbalanced blocks were included except when the animals did not respond to the stimulus in time (the first movement time was longer than 0.8\u2009s). Trials from the 50\u201350 balanced block were excluded from the analysis to avoid possible time artifacts arising from the fact that all of these trials were exclusively recorded in the first 90 trials of the session.<\/p>\n Criteria for neuron inclusion<\/p>\n All neurons were included in the downloaded data provided that their mean firing rate was higher than 0.5\u2009Hz and lower than 50\u2009Hz. For the selectivity and geometry analyses, we included only neurons whose activity was predicted accurately enough by the RRR model described below (above a minimal threshold of \\(\\min \\Delta {R}^{2}\\) with respect to a simple model that assumes that the activity is equal to the average firing rate for all conditions). Unless specified differently, we used \\(\\min \\Delta {R}^{2}=0.015\\). This threshold was necessary to avoid confounding effects from neurons that do not encode any relevant variable, including those recorded with a low signal-to-noise ratio. Although the main results of our article remain qualitatively the same for a broad range of values of \u0394R2 (Extended Data Fig. 6<\/a>), it is important to avoid the extreme cases (that is, when no neuron is discarded, or when too few neurons are selected) when studying whether the representations are categorical or not. Indeed, as already noticed previously13<\/a>, when all neurons are considered, there is the risk that \u2018junk\u2019 neurons with very low selectivity to all variables are over-represented, leading to a peak in the distribution around zero selectivity. This distribution would be significantly different from our null distribution (multivariate Gaussian), but that does not mean the representation is actually categorical or that there is any interesting structure in the selectivity distribution. Indeed, in that previous study13<\/a>, they used as a null distribution the superposition of two Gaussians: one representing the distribution of the selective neurons, and the other, peaked around zero selectivity, to describe the junk neurons. In our case, we decided to discard the worst neurons by selecting only cells that have a large enough \u0394R2. Again, the exact value is not important, but keeping all neurons would lead to an extra peak around zero selectivity, and a misleading, inflated number of brain areas that would pass the criterion for being considered categorical. Again, this is not a real, interesting structure and, in general, when performing this kind of analysis, we recommend checking that the structure revealed by a statistical test is not just due to the overrepresentation of junk neurons. Similarly, if only very few neurons are selected, the centre of the selectivity distribution might be depleted, leading again to the misleading conclusion that the representation is categorical.<\/p>\n RRR encoding model<\/p>\n In this section, we describe the RRR model used to analyse the selectivity profiles of single neurons. We start by describing the input and target variables of the model, followed by a description of the model itself and its fitting procedure. Finally, we introduce a few quantities resulting from the fitted model that are key to the follow-up analysis. The notation that will be used is summarized in the \u2018Notations\u2019 section. The code for implementing and fitting the encoding model is available at GitHub (https:\/\/github.com\/realwsq\/brainwide-RRR-encoding-model<\/a>).<\/p>\n Input and target variablesTarget variables<\/p>\n The target variables (y\u00a0in equation (1<\/a>)) were the preprocessed neuronal responses. The preprocessing steps were applied as follows:<\/p>\n 1.<\/p>\n For each trial, we used the spike trains of the time window \u22120.2 to 0.8\u2009s relative to stimulus onset, as the first movement time is typically less than 0.8\u2009s (ref. 4<\/a>). The activity of each neuron was first binned at 0.01\u2009s, divided by the size of the time bin and then smoothed with a Gaussian filter with a s.d. of 0.02\u2009s. We tried to apply the linear time warping technique51<\/a> so that the stimulus onset time and first movement or response time aligned across trials, but the results did not differ substantially.<\/p>\n 2.<\/p>\n The resulting activity of neuron n, denoted as frn, was organized into a matrix of shape Kn\u00a0\u00d7\u00a0T, where Kn is the number of trials and T\u00a0=\u00a0100 is the number of time steps per trial. Kn depends on n as neurons may have different numbers of trials if they were recorded in different sessions.<\/p>\n 3.<\/p>\n Finally, for each neuron and each time step, we z-scored the activity fr across trials to obtain the target variable y as follows (Extended Data Fig. 1b,c,d<\/a>): <\/p>\n $${y}_{n}(k,t)\\,=\\frac{{\\mathrm{fr}}_{n}(k,t)-{\\mu }_{n}(t)}{{\\sigma }_{n}(t)},\\,\\mathrm{where}\\,{\\mu }_{n}(t)=\\,\\frac{{\\sum }_{k}{\\mathrm{fr}}_{n}(k,t)}{{K}_{n}}\\,\\mathrm{and}\\,{\\sigma }_{n}(t)=\\,\\sqrt{\\frac{{\\sum }_{k}{({\\mathrm{fr}}_{n}(k,t)-{\\mu }_{n}(t))}^{2}}{{K}_{n}}}.$$<\/p>\n \n (1)\n <\/p>\n As the squared error between the preprocessed data and model predictions was used in the loss function for optimizing the model, the applied normalization prevented biases arising from inherent differences in activity scales and ensured that the predictions for all neurons and time steps were optimized equally. Notably, the z-score transformation is easily invertible, allowing the model\u2019s predictions to be mapped back to the original units of firing rate, therefore preserving interpretability (Extended Data Fig. 1b\u2013d<\/a>).<\/p>\n Examples of the processed neural activity are shown in Extended Data Fig. 1b<\/a>.<\/p>\n Input variables<\/p>\n The input variables (x; Extended Data Fig. 1a<\/a>) that we considered can be divided into two types: discrete task-based variables and continuous movement variables. Discrete task-based variables include task-related features, such as the block prior, stimulus contrast, stimulus side, choice and outcome. These are listed below:<\/p>\n Block: the prior probability for the stimulus to appear on the left side is either p(left)\u00a0=\u00a00.2 (right block) or p(left)\u00a0=\u00a00.8 (left block). We used one input variable to encode the block prior: \u22121. representing p(left)\u00a0=\u00a00.2, and +1. representing p(left)\u00a0=\u00a00.8. As noted above, we excluded trials from the p(left)\u00a0=\u00a00.5 unbiased block.<\/p>\n Contrast: the stimulus contrast is 0%, 6.25%, 12.5%, 25% or 100%. One variable was used to encode the stimulus contrast: 0, representing 0% contrast; 1, representing the low contrast (\u226412.5%); and 4, representing the high contrast (>12.5%).<\/p>\n Stimulus: the stimulus location is either on the left side (+1) or the right side (\u22121).<\/p>\n Choice: the choice is indicated by the turning of the wheel: clockwise (+1) or counterclockwise (\u22121).<\/p>\n Outcome: the outcome is either a water reward (+1) or negative feedback (\u22121).<\/p>\n As these values are static, all of the timepoints share the same values (Extended Data Fig. 1a<\/a>).<\/p>\n Continuous movement variables included both instructed (for example, licking and wheel velocity) and uninstructed (for example, whisker motion energy) movement. These were as follows:<\/p>\n Wheel: the velocity of the wheel movement (radian per second) per time bin.<\/p>\n Whisking: the whisker motion energy per time bin is calculated as the motion energy for a square of the left\/right camera roughly covering the whisker pad. The maximum value between the left and right whisker motion energy was used.<\/p>\n Lick: the number of licks per time bin.<\/p>\n A few preprocessing steps were applied separately to each movement variable:<\/p>\n (1)<\/p>\n For each trial, we first read out the continuous behaviour of the time window \u22120.2 to 0.8\u2009s relative to stimulus onset and interpolated into 0.01\u2009s time bins.<\/p>\n (2)<\/p>\n Then, to account for the activity that was shifted in time, for each session, we computed the mean time-lagged correlation between the neuronal activity and the movement traces averaged across neurons and trials and shifted the movement traces so that the zero-lagged correlation was maximized.<\/p>\n (3)<\/p>\n Last, significant differences were observed in the variance of input values across trials, both for different input variables and different time steps. To ensure optimal performance and clarity in interpretation, we z-scored the values for each input variable and time step across trials in the same way as in equation (1<\/a>). The reasoning behind this normalization step is twofold. From a performance perspective, as an extra regularization term was incorporated to penalize the high-value coefficients, the scales of the coefficients and, therefore, the scales of the input values should be comparable for the regularization to work properly. From an interpretability standpoint, the inherent differences in scales must be carefully addressed to enable the comparison of coefficients across input variables and time steps.<\/p>\n Examples of the resulting input variables are shown in Extended Data Fig. 1a<\/a>.<\/p>\n The model formulationLinear encoding model<\/p>\n For each neuron n, we describe its temporal responses as a linear, time-dependent combination of input variables (a visual illustration is shown in Fig. 2a<\/a> and Extended Data Fig. 1<\/a>):<\/p>\n $${y}_{n}(k,t)\\approx {\\widehat{y}}_{n}(k,t)=\\sum _{v}{\\beta }_{n}^{v}(t){x}^{v}(k,t),\\,\\,\\forall k,t,n,$$<\/p>\n \n (2)\n <\/p>\n where<\/p>\n yn(k,t) is the preprocessed neuronal activity of the trial k\u00a0\u2208\u00a0{1,\u2009…,\u2009Kn} and time step t\u00a0\u2208\u00a0{1,\u2009…,\u2009T}.<\/p>\n \\({\\widehat{y}}_{n}(k,t)\\) is the corresponding model prediction, given by the value of the equation on the right side.<\/p>\n v represents the relevant input variables included in the model. xv(k,t) is the preprocessed value of the input variable v for the trial k and time step t.<\/p>\n \\({\\beta }_{n}^{v}(t)\\) is the effect size of the input variable v at time step t. It is further referred to as the regression coefficient.<\/p>\n Low-rank coefficient matrix<\/p>\n The time-varying coefficients \\({{\\boldsymbol{\\beta }}}_{n}^{v}\\in {{\\rm{{\\mathbb{R}}}}}^{T}\\) are the weighted sum of a set of temporal basis vectors shared across all of the neurons and input variables, that is <\/p>\n $${{\\boldsymbol{\\beta }}}_{n}^{v}={{\\bf{U}}}_{n}^{v}{\\bf{V}},\\,\\,\\forall n,v.$$<\/p>\n \n (3)\n <\/p>\n Here, \\({{\\bf{U}}}_{n}^{v}\\in {{\\rm{{\\mathbb{R}}}}}^{d}\\) is the neuron n and input variable v-dependent loading of temporal basis vectors \\({\\bf{V}}\\in {{\\rm{{\\mathbb{R}}}}}^{d\\times T}\\). Specifically, we considered sharing a single set of temporal basis vectors across all of the neurons from all of the brain regions and across all of the input variables. We verified that this restriction did not compromise the goodness-of-fit. The rank d is generally a value much smaller than the number of time steps T. See Extended Data Fig. 1e<\/a> for an example decomposition. Sharing the temporal bases across neurons and input variables significantly reduces the number of parameters (Extended Data Fig. 1f<\/a>). Let N be the number of neurons, T be the number of time steps and \u2223v\u2223 be the number of input variables. An unconstrained full-rank coefficient matrix uses N\u00a0\u00d7\u00a0\u2223v\u2223\u00a0\u00d7\u00a0T parameters, while a reduced-rank coefficient matrix of the same shape only needs N\u00a0\u00d7\u00a0\u2223v\u2223\u00a0\u00d7\u00a0d\u00a0+\u00a0d\u00a0\u00d7\u00a0T parameters. As N and T are typically much larger than d, the reduction in parameters is on the order of T, that is, around 100-fold.<\/p>\n Comparison to previous regression models<\/p>\n Well-established linear encoding models include generalized linear model (GLM)4<\/a>,19<\/a> and kernel regression model20<\/a>. Below, we first describe these two models and then discuss how they compare to our RRR model. GLM is expressed as <\/p>\n $${y}_{n}(k,t)\\approx {\\hat{y}}_{n}(k,t)=\\sum _{v}\\sum _{t{\\prime} }{\\beta }_{n}^{v}(t{\\prime} ){x}^{v}(k,t-t{\\prime} ),\\forall k,t,n,$$<\/p>\n \n (4)\n <\/p>\n where the input filter vector, \\({{\\boldsymbol{\\beta }}}_{n}^{v}\\), composed of \\({\\beta }_{n}^{v}(t{\\prime} )\\) over a neighbouring time window, is factorized as <\/p>\n $${{\\boldsymbol{\\beta }}}_{n}^{v}={{\\bf{U}}}_{n}^{v}{{\\boldsymbol{\\kappa }}}^{v},\\forall n,v.$$<\/p>\n \n (5)\n <\/p>\n Here, \u03bav represents pre-specified temporal basis vectors, which are not trained. Adapted from a previous study4<\/a> (here we considered neural responses over a different time window from4<\/a> and an\u00a0enriched set of behaviour movements), we used the same set of raised cosine \u2018bump\u2019 functions in log space for Extended Data Fig. 1h<\/a>. The kernel regression model uses the same linear formulation as the GLM (equation (4<\/a>)). However, instead of relying on pre-defined temporal basis vectors, it allows these vectors to be trainable: <\/p>\n $${{\\boldsymbol{\\beta }}}_{n}^{v}={{\\bf{U}}}_{n}^{v}{{\\bf{V}}}^{v},\\forall n,v.$$<\/p>\n \n (6)\n <\/p>\n Both the kernel regression model and our RRR model fall into the category of RRR models. The primary distinction lies in the set of input features used by each approach.<\/p>\n In our RRR model, the influence of input variables on neural responses, \\({\\beta }_{n}^{v}(t)\\) can vary over time. By contrast, both the GLM and the kernel regression model assume that the influence is time-independent. Therefore, for example, whether the mouse response time is early or late makes no difference and will modulate the neural responses the same way. Instead, these models provide a more descriptive account of input effects by allowing neural responses to depend on inputs from neighbouring time steps. While this feature is absent in our current RRR model (equation (2<\/a>)), it could be incorporated in future extensions. A performance comparison is shown in Extended Data Fig. 1h<\/a>.<\/p>\n Estimation of the parameters<\/p>\n The parameters of the RRR model include a shared temporal bases matrix V of size d\u00a0\u00d7\u00a0T and loading vectors \\({{\\bf{U}}}_{n}^{v}\\) of length d for each input variable v and neuron n. The approach we adopted to fit the parameters was to minimize the ridge-penalized mean square loss: <\/p>\n $${\\mathcal{L}}({\\bf{V}},{\\{{{\\bf{U}}}_{n}^{v}\\}}_{n,v})=\\sum _{n}(\\sum _{k}\\sum _{t}{({y}_{n}(k,t)-{\\hat{y}}_{n}(k,t))}^{2}+\\lambda \\sum _{v}\\sum _{t}{\\beta }_{n}^{v}{(t)}^{2}).$$<\/p>\n \n (7)\n <\/p>\n Minimizing this particular loss function is straightforward as a closed-form solution exists52<\/a>. In practice, we chose to use the L-BFGS optimization algorithm to compute the optimum.<\/p>\n Moreover, to optimize the model hyperparameters, namely the rank d and the regularization penalty \u03bb, we implemented a threefold cross-validation technique across trials. First, the dataset was stratified based on a composite target label that included the block prior and stimulus contrast to ensure that each fold was representative of the entire dataset. Then, for each combination of d and \u03bb, the dataset was partitioned into three subsets by trials, using each subset in turn for testing the model while the remaining data served as the training set. Finally, the combination of d and \u03bb that yielded the lowest average test error across all folds was selected. d\u00a0=\u00a05 turned out to be the optimal number of temporal bases.<\/p>\n Estimating the goodness of fit<\/p>\n We used the threefold cross-validated R2 to measure the goodness-of-fit of single-trial predictions. For each session\u2019s data, we randomly sampled one-third of the trials as the test set held out during training. Once the model was trained\u2014using the remaining two-thirds of trials\u2014we computed the R2 between the model predictions \\({\\hat{\\mathrm{fr}}}_{n}(k,t)\\) ((\\({\\hat{\\mathrm{fr}}}_{n}(k,t)\\) is calculated by inversing the z-score transformation applied in the preprocessing step \\({\\hat{\\mathrm{fr}}}_{n}(k,t)={\\sigma }_{n}(t){\\hat{y}}_{n}(k,t)+{\\mu }_{n}(t)\\); Extended Data Fig. 1b\u2013d<\/a>) and the actual neuronal activity frn(k,t) in the test set. We repeated the whole split\u2013train\u2013test process three times and computed the mean of the three cross-validated R2 as the measure of goodness-of-fit.<\/p>\n Null model and selectively modulated neurons<\/p>\n Conceptually, we distinguish two types of task modulation: the selective modulation and the non-selective modulation. Selective modulation, captured by \\({\\hat{y}}_{n}(k,t)={\\sum }_{v}{\\beta }_{n}^{v}(t){x}^{v}(k,t)\\) is induced by the input variables and varied trial by trial. Non-selective modulation, captured by the mean time-varying response \\({\\mu }_{n}(t)=\\frac{{\\sum }_{k}{\\mathrm{fr}}_{n}(k,t)}{{K}_{n}}\\) is locked to the key events of the trial (stimulus onset in this case) and does not vary trial by trial. Both types have an important role in modulating neuronal responses. See Extended Data Fig. 2b<\/a> for example selectively modulated (left) and non-selectively modulated (right) neurons. In this work, we focus mostly on the neuronal responses selectively modulated by the task. To distinguish the variation explained by the selective modulation from the non-selective modulation, we use the trial-average estimate as the null model that does not consider any effect of the input variables: <\/p>\n $${\\hat{y}}_{n}^{{\\rm{null}}}(k,t)=0,\\,\\forall k,t,n.$$<\/p>\n \n (8)\n <\/p>\n The outperformance, \u0394R2, defined as: <\/p>\n $$\\Delta {R}^{2}({\\rm{model}})={R}^{2}({\\rm{model}})-{R}^{2}({\\rm{null}}),$$<\/p>\n \n (9)\n <\/p>\n captures the overall selective modulation of all of the input variables combined.<\/p>\n Only selectively modulated neurons, identified as \u0394R2(RRR)\u2009\u2265\u20090.015 (Extended Data Fig. 1g<\/a>), are included in the further analysis.<\/p>\n Computing the selectivity profiles of single neurons to the input variable<\/p>\n To compute the selectivity profiles of single neurons to the individual variables, we used the estimated coefficient \\({\\beta }_{n}^{v}(t)\\). For the clustering analysis of Fig. 3<\/a> and Extended Data Fig. 6<\/a>, we took the sum of the coefficients across time as a measure of the total selectivity of neuron n to input variable v.<\/p>\n $${\\alpha }_{n}^{v}=\\sum _{t}{\\beta }_{n}^{v}(t)$$<\/p>\n \n (10)\n <\/p>\n Note that by normalizing the neuronal responses and input variables in the preprocessing steps, we ensured that the unit-free coefficients \\({\\beta }_{n}^{v}(t)\\) are not affected by the neuron\u2019s mean firing rate or the inherently different scales in different input variables and can be compared directly across neurons, input variables and time steps. Thus, \\({\\beta }_{n}^{v}(t)\\) can be interpreted as the expected change in normalized neuronal activity yn per one s.d. change in the input variable xv at time t, and \\({\\alpha }_{n}^{v}\\) can therefore be thought as the expected total change across the whole trial.<\/p>\n The selectivity \\({\\alpha }_{n}^{v}\\) captures whether individual neurons are selectively modulated by the given variable v or not (examples of strongly selective neurons are shown in Extended Data Fig. 2c<\/a>). In the selectivity analysis (Fig. 2e,f<\/a>), when the goal is to estimate the absolute modulation of an input variable, \\({\\alpha }_{n}^{v}\\) is calculated as \\({\\alpha }_{n}^{v}={\\sum }_{t}| {\\beta }_{n}^{v}(t)| \\).<\/p>\n Computing the autocorrelation timescale of neural responses to task and behaviour variables<\/p>\n Given a matrix of single-neuron responses to task and behaviour variables \\({\\hat{y}}_{n}\\in {{\\rm{{\\mathbb{R}}}}}^{K\\times T}\\) with K being the number of trials and T the number of time steps per trial, we can compute the corresponding autocorrelation timescale.<\/p>\n To compute the timescale, we first calculate the time-lagged unnormalized autocorrelation sequence <\/p>\n $${c}_{n}(i)=\\frac{{\\sum }_{k}{\\sum }_{t}{\\hat{y}}_{n}(k,t){\\hat{y}}_{n}(k,t+i)}{K},\\,i\\ge 0.$$<\/p>\n We then linearly interpolate the autocorrelation sequence so that cn(i) is spaced at 1\u2009ms resolution (original 10\u2009ms). The timescale \u03c4n is approximated by the time the sequence first reaches half its peak value (that is, cn(0)). The timescale of brain area a is further determined by averaging the values over all of the selectively modulated neurons within this area.<\/p>\n We observed a significant correlation between a region\u2019s hierarchical position and its estimated autocorrelation timescale (Extended Data Fig. 2d<\/a>).<\/p>\n Notations<\/p>\n Indices:<\/p>\n n, N: index and number of neurons (n\u00a0\u2208\u00a0{1,\u00a0\u2026,\u00a0N}).<\/p>\n k, K: index and number of trials (k\u00a0\u2208\u00a0{1,\u00a0\u2026,\u00a0K}).<\/p>\n t, T: index and number of time steps (t\u00a0\u2208\u00a0{1,\u00a0\u2026,\u00a0T}).<\/p>\n v, \\(| v| \\): index and number of input variables (v\u00a0\u2208\u2009{block prior, stimulus contrast, stimulus side, choice, outcome, wheel velocity, whisker motion energy, lick}).<\/p>\n i, d : index and rank of the RRR model (i\u00a0\u2208\u00a0{1,\u00a0\u2026,\u00a0d}).<\/p>\n RRR model<\/p>\n \\({\\beta }_{n}^{v}(t)\\): regression coefficient (arbitrary units (a.u.)) of neuron n, input variable v at time step t.<\/p>\n \\({{\\bf{U}}}_{n}^{v}\\in {{\\mathbb{R}}}^{d}\\): neuron n and input variable v-dependent loading of temporal basis vectors (a.u.).<\/p>\n \\({\\bf{V}}\\in {{\\rm{{\\mathbb{R}}}}}^{d\\times T}\\): temporal basis vectors (a.u.).<\/p>\n Random variables:<\/p>\n xv(k,t): preprocessed value (a.u.) of input variable v, trial k at time step t.<\/p>\n yn(k,t): preprocessed neuronal response (a.u.) of neuron n, trial k at time step t.<\/p>\n \\({\\widehat{y}}_{n}(k,t)\\): model prediction of preprocessed neuronal response (a.u.) of neuron n, trial k at time step t.<\/p>\n frn(k,\u00a0t): smoothed, binned firing rate (Hz) of neuron n, trial k at time step t.<\/p>\n \\({\\hat{\\mathrm{fr}}}_{n}(k,t)\\): model prediction of smoothed, binned firing rate (Hz) of neuron n, trial k at time step t.<\/p>\n \u03bcn(t): mean of smoothed, binned firing rate (Hz) of neuron n at time step t.<\/p>\n \u03c3n(t): s.d. of smoothed, binned firing rate (Hz) of neuron n at time step t.<\/p>\n Clustering analysis<\/p>\n To test for the presence of functional clusters, we followed the steps explained below. The required inputs include each relevant neuron\u2019s response profile and original session ID. Two types of response profile can be considered: the estimated selectivity to individual input variables (equation (10<\/a>), referred to as clustering analysis in the variable selectivity space, used in Fig. 3<\/a> and Extended Data Fig. 6<\/a>) or the average response in each experimental condition (referred to as clustering analysis in the conditions space, used in Extended Data Fig. 6<\/a>). Performing clustering analysis in the selectivity space arguably has a few advantages: (1) It reduces the dimensionality in an interpretable and informed way. If we have \u2223v\u2223 variables, then there are at least 2\u2223v\u2223 conditions, assuming all of the variables are discrete and have more than one different value. (2) It mitigates the issue of unbalanced or even missing conditions. (3) It reduces the noise in the estimation of the response profile. As shown in Fig. 2b<\/a> and Extended Data Fig. 2<\/a>, neural responses are very noisy, and simple averaging may be non-satisfactory. By contrast, the selectivity estimated from the encoding model provides a more reliable account of the task-driven variance in the neural responses.<\/p>\n The code for the clustering analysis is available at GitHub (https:\/\/github.com\/realwsq\/clustering-analysis<\/a>).<\/p>\n Clustering analysis in the variable selectivity space<\/p>\n We summarize our clustering pipeline as follows (a schematic is shown in Fig. 3a<\/a>). (1) Check whether there are more than 50 neurons and only continue if so. (2) Given the selectivity profile of each neuron, run the k-means clustering algorithm (with 100 random initializations) with the number of clusters k varied from 3 to 20. (3) Then, select the optimal clustering result by maximizing the silhouette score. The silhouette score is defined as \\({ < \\frac{{b}_{i}-{a}_{i}}{\\max ({b}_{i},{a}_{i})} > }_{i}\\) where i is the index of the neuron, \\({a}_{i}=\\frac{1}{| {C}_{I}| -1}{\\sum }_{j\\in {C}_{I},i\\ne j}d(i,j)\\) is the mean Euclidean distance intracluster and \\({b}_{i}={\\min }_{J\\ne I}\\frac{1}{| {C}_{J}|}{\\sum }_{j\\in {C}_{J}}d(i,j)\\) is the minimum Euclidean distance outside cluster (Fig. 3a<\/a>). (4) Iterate over the resulting clusters and check whether there is a cluster whose total silhouette score summed over all neurons is mainly contributed by neurons from one single session (>90%). If so, remove neurons from that cluster and session, and repeat steps 1\u20134. (5) Sample the same number of datapoints from the Gaussian distribution with the mean and covariance matrix matched to the data values and compute the sampled data\u2019s null silhouette score according to steps 2 and 3. (6) Repeat step 5 100 times and pool the null silhouette scores to form the null distribution. Finally, compute the z score of the data silhouette score with respect to the null distribution.<\/p>\n Clustering analysis in the conditions space<\/p>\n When clustering in the space of mean firing rate, two additional preprocessing steps are required. First, we normalized each neuron\u2019s mean firing rates separately across conditions to prevent clustering driven solely by overall firing rate differences between neurons. Second, as the number of conditions and the dimensionality of the activity profile is high, we reduced the dimensionality using principal component analysis. Moreover, we modified the null model by replacing the multivariate Gaussian distribution with a multivariate log-normal distribution to better capture the lower-bounded and heavy-tailed nature of mean firing rates. Neurons with zero activity were excluded from this analysis.<\/p>\n The clustering analysis in the conditions space follows these steps (a schematic is shown in Extended Data Fig. 5<\/a>): (1) check whether there are more than 40 neurons and only continue if so. (2) z-Score the mean firing rates for each neuron. (3) Reduce the dimensionality using principal component analysis, retaining components that capture 90% of the total variance. (4) Run the k-means clustering algorithm (with 100 random initializations), varying the number of clusters k from 3 to 20. (5) Select the optimal clustering result by maximizing the silhouette score. (6) Iterate over the resulting clusters and check whether there is a cluster whose total silhouette score summed over all neurons is mainly contributed by neurons from one single session (>90%). If so, remove neurons from that cluster and session, and repeat steps 1\u20135. (7) Sample the same number of datapoints from the multivariate log-normal distribution with the mean and covariance matrix matched to the log data values, then compute the sampled data\u2019s null silhouette score following steps 2\u20135. (8) Repeat step 7 100 times and pool the null silhouette scores to form the null distribution. Finally, compute the z score of the data silhouette score with respect to the null distribution.<\/p>\n Measuring the similarity between cluster and area labels<\/p>\n The similarity between functional clusters and anatomical area labels was quantified using the RI, which measures the agreement between two labellings of the same dataset. For each pair of neurons, agreement occurs if both are assigned to the same cluster in both labellings, or to different clusters in both. The RI is defined as the fraction of agreeing pairs among all possible pairs. To assess statistical significance, we computed a z-scored RI by comparing the observed RI to a null distribution obtained from 10,000 random shufflings of the area labels. A significantly elevated z-scored RI indicates that functional clustering aligns closely with anatomical organization. To avoid biases towards areas or modules with larger neuron numbers, we considered the 100 best-encoded neurons from each module or area in the analysis of Fig. 3e,f<\/a> and Extended Data Fig. 7<\/a>.<\/p>\n Modified ePAIRS test<\/p>\n Given the average responses of single neurons in each experimental condition, we performed the ePAIRS test as follows (a schematic is shown in Extended Data Fig. 5<\/a>): (1) z-score the mean firing rates for each neuron. (2) Reduce the dimensionality using principal component analysis, retaining components that capture 90% of the total variance. (3) Calculate the cosine distance to its nearest neighbour for each neuron. (4) Calculate the empirical median of these nearest-neighbour distances as the aggregated nearest-neighbour angle. (5) Sample the same number of datapoints from the multivariate log-normal distribution with the mean and covariance matrix matched to the log data values, then compute the null nearest-neighbour angle of the sampled data according to steps 1\u20134. (6) Repeat step 5 5,000 times and pool the null nearest-neighbour angles to form the null distribution. Finally, the z-score of the data aggregated nearest-neighbour angle with respect to the null distribution is computed.<\/p>\n \u03b1-Diversity<\/p>\n To measure \u03b1-diversity, we took the participation ratio of the N\u00a0\u00d7\u00a0V matrix of \u03b1 coefficients resulting from the RRR analysis described above. The PR quantifies the effective dimensionality of a set of datapoints by measuring how evenly the variance is distributed across the eigenvalues of its principal component analysis decomposition28<\/a>. The PR is defined as: <\/p>\n $${\\rm{PR}}=\\frac{{\\left({\\sum }_{j}{\\lambda }_{j}\\right)}^{2}}{{\\sum }_{j}{\\lambda }_{j}^{2}},$$<\/p>\n \n (11)\n <\/p>\n where \u03bbj is the jth eigenvalue of the N\u00a0\u00d7\u00a0N covariance matrix. A higher PR indicates that the variance is more evenly spread across multiple dimensions, suggesting a higher effective dimensionality of the cloud of points in the high-dimensional space. Conversely, a lower PR implies that the variance is concentrated in fewer dimensions, indicating a lower effective dimensionality.<\/p>\n The number of neurons in individual areas was soft-equalized by taking a random subsample of N0\u00a0=\u00a0120 neurons when N\u00a0>\u00a0N0. In these cases, the participation ratio was computed over 100 random subsets, and the average was taken as the value of \u03b1-diversity.<\/p>\n Analysis of population representationsData preparation<\/p>\n For the analysis of population neural representations, we used the same sessions as described above. For each trial within a session, we therefore have a collection of N-dimensional population activity vectors ft,k, where k\u00a0\u2208\u00a0{1,\u00a0P} is the trial index within the session, and t\u00a0\u2208\u00a0{1,\u00a0T} is the time-bin index within each trial. For the analysis below, we used data from 0 to 1,000\u2009ms after the stimulus onset to capture a variety of sensory and behavioural variables. We then labelled each time bin according to the value of four binarized cognitive, sensory and movement variables:<\/p>\n Block: left (20\u201380) versus right (80\u201320) prior block.<\/p>\n Contrast: we binarized the contrast into low (0\u20130.125) versus high (0.25\u20131.0) values.<\/p>\n Stimulus: left versus right side of the screen.<\/p>\n Whisking: we binarized the whisking power using the distribution of whisking power values within each session. Time bins in which the mouse was whisking with a power larger than the 50th percentile across the distribution were annotated as high, while those below the 50th percentile were annotated as low.<\/p>\n These variables were chosen so that they span movement, cognitive and sensory variables while ensuring that all of the M\u00a0=\u00a016 conditions (combinations of the four variables) were well represented in the data. For example, we could not add Choice as a variable as mice are overtrained in the task and, as a consequence, make very few mistakes when block and stimulus are aligned (for example, choose \u2018left\u2019 when the block and the stimulus are both \u2018right\u2019).<\/p>\n For each condition c\u00a0\u2208\u00a0{1,\u00a0M} (for example, whisking\u00a0=\u2009high, contrast\u00a0=\u2009low, stimulus\u00a0=\u2009left, block\u00a0=\u2009right), we first identified those trials in which that specific combination of variables was present. We then defined a collection of \u2018conditioned trial\u2019 population activity vectors {fk,c} as the mean firing rate of the population of neurons conditioned to the specific condition in each trial. Given a trial k and a neuron index i, the mean firing rate was computed as <\/p>\n $${f}_{i}^{k,c}=\\frac{{\\sum }_{t}{\\delta }^{t,k}(c)\\,{f}_{i}^{t,k}}{{\\sum }_{t}{\\delta }^{t,k}(c)}\\,,$$<\/p>\n \n (12)\n <\/p>\n where i indicates the neuron index and \u03b4t,k(c)\u00a0=\u00a01 if the time bin t in trial k corresponds to the condition c, and 0 otherwise. These conditioned trial population vectors are the data samples that will be used for the dimensionality and decoding analyses below. Across all analyses, we considered only those recording sessions in which each condition was present in at least Mmin\u00a0=\u00a05 trials.<\/p>\n Representation dimensionality<\/p>\n To estimate the representation dimensionality of a neural geometry, we computed the PR of the centroids fc of the M conditions, defined as the average activity pattern across all trials of the same condition: <\/p>\n $${{\\bf{f}}}^{c}={ < {{\\bf{f}}}^{k,c} > }_{k},$$<\/p>\n \n (13)\n <\/p>\n To compute the PR for the set of centroids, we first calculated their covariance matrix, normalizing each neuron\u2019s mean activity vector by subtracting from each \\({f}_{i}^{c}\\) the mean across conditions c and dividing them by their s.d. We then performed principal component analysis on this covariance matrix to obtain its eigenvalues, \u03bbj.<\/p>\n The number of neurons in individual areas was soft-equalized by taking a random subsample of N0\u00a0=\u00a0120 neurons when N\u00a0>\u00a0N0. In these cases, the participation ratio was computed over 100 random subsets, and the average was taken as the value of representation dimensionality. The PR was computed on the set of centroids to highlight signal dimensions and prevent noise from dominating the measure. When computed across all trial vectors (without condition averaging), the PR ranged from 10 to 350 and was highly correlated with that of the centroids (Spearman Correlation coefficient\u2009=\u20090.86).<\/p>\n Cross-validated decoding<\/p>\n We used Decodanda53<\/a> (www.github.com\/lposani\/decodanda<\/a>) to perform a cross-validated, class-balanced decoding analysis of different combination of condition labels from the neural activity within individual trials (condition trial vectors fk,c). See the individual sections below for additional details on the data input structure of our decoding analyses. As a decoder, we used a scikit-learn SVM classifier with linear kernel54<\/a>. To ensure that results were comparable across regions, which might have a different number of recorded neurons, we created a pseudopopulation by resampling all of the recorded neurons within each region to a fixed number N\u00a0=\u00a04,000. Similarly, we resampled the same number of pseudopopulation for each analysis (T\u00a0=\u00a0100 patterns per condition). Note that simultaneously recorded neurons were always kept together during resampling to keep the noise correlations intact within the pseudopopulation53<\/a>. All cross-validated decoding analyses were performed using the following Decodanda parameters: training_fraction\u00a0=\u00a00.8, cross_validations\u00a0=\u00a0100, ndata\u00a0=\u00a0100.<\/p>\n Finding the independent conditions<\/p>\n To find the number of independent conditions encoded in the activity of a population of neurons, we developed an iterative algorithm based on linear decoding. The algorithm followed the steps below, and is shown in action on one example region in Extended Data Fig. 3<\/a>.<\/p>\n (1)<\/p>\n First, we performed a decoding analysis of the condition label c from trial population vectors fk,c using a set of binary linear classifiers. For each pair of conditions (ci,\u00a0cj), we estimated a cross-validated decoding performance \u03c6(ci,\u00a0cj), resulting in an initial M\u00a0\u00d7\u00a0M condition\u2013condition decoding matrix (C0; Extended Data Fig. 3<\/a>) defined as C0(ij)\u00a0=\u00a0\u03c6(ci,\u00a0cj).<\/p>\n (2)<\/p>\n We then chose a decoding threshold \\({\\varphi }_{\\min }=0.666\\); the pairs of conditions whose one-versus-one decoding performance was smaller than \\({\\varphi }_{\\min }\\) were defined as dependent. Using this threshold, we defined a binary dependency matrix D defined as D0(i,\u00a0j)\u00a0=\u00a01 if \\(\\varphi ({c}_{i},{c}_{j}) < {\\varphi }_{\\min }\\), and C0(i,\u00a0j)\u00a0=\u00a00 otherwise.<\/p>\n (3)<\/p>\n We then used the Bron\u2013Kerbosch algorithm55<\/a> to find all of the cliques, that is, subgroups of fully connected nodes, in the undirected graph defined by the dependency matrix D0. This process enables us to identify whether there are groups of conditions that are all non-decodable from each other (dark squares in the sorted matrix in Extended Data Fig. 3<\/a>).<\/p>\n (4)<\/p>\n We next identified the largest clique and grouped together all of the trials of the conditions within that group into a new, merged condition (see the arrows and \u2018merge\u2019 conditions in Extended Data Fig. 3<\/a>).<\/p>\n (5)<\/p>\n We repeated steps 1 and 2 with the new reduced set of conditions, yielding a new Ct and a new Dt matrix of a different size Mt, where t denotes the iteration step.<\/p>\n (6)<\/p>\n We then repeated steps 3 and 4, and iterated the whole process (1\u20134) until all of the merged and remaining conditions were found to be independent, that is, the dependency matrix \\({D}_{\\widetilde{{t}}}\\) was diagonal. The number of independent conditions was then defined as the size of the final dependency matrix: \\({M}_{\\mathrm{IC}}:={M}_{\\widetilde{t}}\\).<\/p>\n Separability and AD<\/p>\n Separability quantifies how many random dichotomies (equally sized groups) of experimental conditions can be decoded from neural activity using cross-validated linear classifiers. To estimate the separability of a neural population, we performed the following steps:<\/p>\n 1.<\/p>\n First, we randomly divided the set of M or MIC independent conditions into two equally sized groups (dichotomy).<\/p>\n 2.<\/p>\n Given the dichotomy d, we then measured the cross-validated decoding performance \u03c6d of a linear classifier trained to report whether individual condition trial vectors fk,c belonged to conditions within one or the other dichotomy groups. This decoding analysis was performed as described in the \u2018Cross-validated\u00a0decoding\u2019 section above, resampling a fixed large number of neurons (N\u00a0=\u00a04,000) and a fixed number of trials (T\u00a0=\u00a0100) per condition for all regions to ensure that performances could be compared across regions.<\/p>\n 3.<\/p>\n The random dichotomy assignment and decoding (steps 1 and 2) was then repeated n\u00a0=\u00a0200 times to obtain a set of decoding performances {\u03c6d}.<\/p>\n 4.<\/p>\n The decoding analysis was then repeated n\u00a0=\u00a0200 times with shuffled condition labels across the population vectors to obtain a distribution of null decoding performance values {\u03c6null}.<\/p>\n 5.<\/p>\n Separability was then defined as the fraction of decodable dichotomies, that is, the fraction of \u03c6d larger than the 99 percentile of the null population {\u03c6null}. AD was defined as the mean decoding performance across the n\u00a0=\u00a0200 random dichotomies.<\/p>\n The distributions of decoding performances for all of the analysed regions are shown in Extended Data Fig. 11<\/a>.<\/p>\n Synthetic dataModelling uneven and categorical selectivity<\/p>\n We generated synthetic population responses for N neurons to all 2V configurations of V binary variables xv\u00a0\u2208\u00a0{\u00a0\u2212\u00a01,\u00a01}. Each neuron\u2019s response combined linear and quadratic selectivity, <\/p>\n $${f}_{i}(\\vec{x})=\\mathop{\\sum }\\limits_{v=1}^{V}{\\alpha }_{iv}\\,{x}_{v}+\\gamma \\sum _{u < v}{\\beta }_{i,uv}\\,{x}_{u}{x}_{v}+{\\epsilon },$$<\/p>\n \n (14)\n <\/p>\n where \u03b3 controls the strength of non-linear interactions and \\({{\\epsilon }}_{i} \\sim {\\mathcal{N}}(0,\\sigma )\\) introduces trial-to-trial variability. For each condition \\(\\vec{x}\\), we generated T independent trials and used the data in the decoding analyses as performed for real spiking data.<\/p>\n Sampling the selectivity structure<\/p>\n The coefficients (\u03b1iv,\u00a0\u03b2i,uv) were sampled in a feature space of dimension \\(D=V+\\frac{V(V-1)}{2}\\), which groups all linear and quadratic features on an equal footing. To generate either categorical or uneven selectivity, we first drew k cluster centroids in the D-dimensional feature space and repeated each centroid kN\u00a0=\u00a0N\/k times, yielding N prototype vectors. These prototypes define the coarse structure of selectivity across neurons. To modulate within-cluster diversity, each prototype was perturbed by an additive diversity vector <\/p>\n $${\\delta }_{i}={\\sigma }_{d}\\,{\\xi }_{i},{\\xi }_{i} \\sim {\\mathcal{N}}(0,{I}_{D}),$$<\/p>\n where \\({\\sigma }_{d}=-\\log (c)\\) is set by a categorical specialization parameter c\u00a0\u2208\u00a0(0,\u00a01) (this is the x axis in Extended Data Fig. 8b<\/a>). Larger c produces more tightly clustered selectivity (strong categorical structure), whereas smaller c yields more dispersed coefficients. To model uneven selectivity, we introduced anisotropy across the D feature dimensions by using k\u00a0=\u00a01 (no clustering) and scaling the diversity terms by a geometric decay profile: <\/p>\n $${\\sigma }_{j}={\\sigma }_{d}\\,\\log {(1-r)}^{j-1},j=1,\\ldots ,D,$$<\/p>\n where r\u00a0\u2208\u00a0(0,\u00a01) is the global specialization parameter. When r\u00a0\u226a\u00a01, all dimensions contribute equally; when r\u00a0\u2248\u00a01, variance is concentrated in a low-dimensional subspace, producing a strongly uneven selectivity spectrum. The final coefficients, including categorical and\/or uneven structure, were <\/p>\n $$({\\alpha }_{i},{\\beta }_{i})={\\mathrm{centroid}}_{i}\\,+\\,\\sigma \\odot {\\xi }_{i},$$<\/p>\n with the quadratic coefficients additionally scaled by \u03b3 to control non-linearity.<\/p>\n Trial generation<\/p>\n For each of the 2V binary stimuli, we computed \\({f}_{i}(\\vec{x})\\) through the linear-quadratic form above, and generated T noisy samples per condition, <\/p>\n $${r}_{i,t}(\\vec{x})={f}_{i}(\\vec{x})+{{\\epsilon }}_{i,t},{{\\epsilon }}_{i,t} \\sim {\\mathcal{N}}(0,\\sigma ).$$<\/p>\n This procedure ensures that trial variability is independent across neurons and conditions, and that all structure in the population code arises exclusively from the geometry of the sampled coefficients.<\/p>\n Parameterizing categorical and uneven specialization<\/p>\n The two specialization parameters (c,\u00a0r) therefore independently control: the categorical structure of selectivity (the number and separation of clusters, through c and k); and the unevenness of selectivity across feature dimensions (anisotropy of coefficient variances, through r).<\/p>\n By sweeping either c (categorical specialization) or r (global specialization) while keeping the other fixed, we isolated the effects of clustered versus uneven selectivity on representational dimensionality, separability and independent condition structure. For the analyses in Extended Data Fig. 8<\/a>, we used N\u00a0=\u00a0100 neurons, V\u00a0=\u00a03 variables (for a total of P\u00a0=\u00a08 conditions), T\u00a0=\u00a020 samples per condition and \u03b3\u00a0=\u00a00.25.<\/p>\n Exploring the relationship between dimensionality and separability<\/p>\n To analyse how separability changes with the dimensionality of the geometry in the activity space, we performed a series of synthetic explorations shown in Extended Data Fig. 8<\/a>. In these simulations, P centroids are randomly sampled from a Gaussian distribution spanning an L-dimensional subspace of the N-dimensional activity space. Each centroid vector is normalized. Trial-to-trial variability is then added to the centroids with a s.d. \u03c3, scaled with \\(\\sqrt{L}\\) to keep the signal-to-noise level constant when pairwise distances between centroids increase with L. This obtains a T\u2009\u00d7\u2009N activity matrix for each of the P conditions. This synthetic activity is then analysed with the same pipeline used for the cortical data, yielding values of representation dimensionality, separability and AD shown in Extended Data Fig. 8g<\/a>, in which we used P\u2009=\u200916, N\u2009=\u2009100. For simulations in Extended Data Fig. 8i<\/a>, we fixed L\u2009=\u200914 and multiplied the first dimension of each centroid by a factor \u03b3 to stretch the geometry along a single axis.<\/p>\n Theoretical considerations on the relationship between dimensionality and clusteringThe conditions space and the neural space have the same dimensionality<\/p>\n As explained in Fig. 1a<\/a>, response profiles of single neurons can be thought of as rows of a matrix X for which the columns define the geometry of conditions in the neural space. The PR in the conditions space is computed from the eigenvalue spectrum of the covariance matrix of the rows of X, that is, XXT (assuming zero mean), while the PR in the activity space is computed from the eigenvalues of the covariance matrix of the columns of X, that is, XTX. Given a matrix, the eigenvalues of its covariance matrix are the squared singular values of X. Let X\u00a0=\u00a0USVT be the SVD decomposition of X, where S is the diagonal matrix with singular values on the diagonal, then XT\u00a0=\u00a0VSTUT. As S\u00a0=\u00a0ST, the eigenvalue spectrum of XTX and XXT is the same. Thus, the participation ratio of the conditions space is the same as that in the neural space.<\/p>\n Mathematical derivation of the PR of Gaussian clusters<\/p>\n We consider a data model with N features (neurons) and M observations (conditions), in which observations are\u00a0sampled from independent and identically distributed (i.i.d.)\u00a0random variables\u00a0as <\/p>\n $${{\\bf{x}}}^{\\mu }={{\\bf{z}}}^{\\mu }+{{\\boldsymbol{\\eta }}}^{\\mu },$$<\/p>\n \n (15)\n <\/p>\n where z\u03bc and \u03b7\u03bc are both vectors in \\({{\\mathbb{R}}}^{N}\\) and represent the clustered and heterogeneous part of the data, respectively. More precisely, z\u03bc is sampled from a normal distribution \\({\\mathcal{N}}(0,{\\bf{B}})\\) that has a clustered covariance matrix, that is, Bij\u00a0=\u00a01 if i and j belong to the same cluster and Bij\u00a0=\u00a00 otherwise. We call k the number of clusters and assume that all clusters have the same number of neurons Nc\u00a0=\u00a0N\/k. By contrast, the heterogenous part \u03b7\u03bc is sampled from \\({\\mathcal{N}}(0,{\\sigma }^{2}{\\bf{I}})\\), where I is the identity matrix. Our goal is to compute the participation ratio (PR) of this representation, which we define as <\/p>\n $$\\mathrm{PR}=\\frac{\\mathrm{Tr}{({\\bf{C}})}^{2}}{\\mathrm{Tr}{{\\bf{C}}}^{2}}=\\frac{N{\\langle {C}_{ii}\\rangle }^{2}}{\\langle {C}_{ii}^{2}\\rangle +(N-1)\\langle {C}_{ij}^{2}\\rangle },$$<\/p>\n \n (16)\n <\/p>\n where the averages are across neurons and the matrix C is the sample neuron-by-neuron covariance matrix, that is \\(C=\\frac{1}{M}{\\sum }_{\\mu =1}^{M}{{\\bf{x}}}^{\\mu }{({{\\bf{x}}}^{\\mu })}^{T}\\). We note that this definition assumes that the sample mean of both z and \u03b7 are negligible or have been subtracted.<\/p>\n We are interested in the regime in which N\u00a0\u2192\u00a0\u221e while M is allowed to be small, as it often happens in controlled experiments. Small M might cause the sample covariance matrix to differ substantially from the true covariance matrix. Defining Cc, Ch, and Cch as the sample covariance matrices of z, \u03b7, and the cross-covariance between z and \u03b7, respectively, we have that <\/p>\n $$\\mathrm{PR}=N\\frac{{\\langle {C}_{ii}^{{\\rm{h}}}+{C}_{ii}^{{\\rm{c}}}+2{C}_{ii}^{\\mathrm{ch}}\\rangle }^{2}}{\\langle {({C}_{ii}^{{\\rm{h}}}+{C}_{ii}^{{\\rm{c}}}+2{C}_{ii}^{\\mathrm{ch}})}^{2}\\rangle +(N-1)\\langle {({C}_{ij}^{{\\rm{h}}}+{C}_{ij}^{{\\rm{c}}}+2{C}_{ij}^{\\mathrm{ch}})}^{2}\\rangle }.$$<\/p>\n \n (17)\n <\/p>\n We therefore need to evaluate the first and second moments of both diagonal and off-diagonal elements of all covariance and cross-covariance matrices. The diagonal elements of these matrices have the following statistics: <\/p>\n $$\\begin{array}{c}{C}_{ii}^{{\\rm{c}}}=\\frac{1}{M}\\mathop{\\sum }\\limits_{\\mu =1}^{M}{({z}_{i}^{\\mu })}^{2}\\,\\Rightarrow \\,\\langle {C}_{ii}^{{\\rm{c}}}\\rangle =1\\,,\\quad \\langle {({C}_{ii}^{{\\rm{c}}})}^{2}\\rangle =\\frac{M+2}{M}\\\\ {C}_{ii}^{{\\rm{h}}}=\\frac{1}{M}\\mathop{\\sum }\\limits_{\\mu =1}^{M}{({\\eta }_{i}^{\\mu })}^{2}\\,\\Rightarrow \\,\\langle {C}_{ii}^{{\\rm{h}}}\\rangle ={\\sigma }^{2}\\,,\\quad \\langle {({C}_{ii}^{{\\rm{h}}})}^{2}\\rangle =\\frac{M+2}{M}{\\sigma }^{4}\\\\ {C}_{ii}^{\\mathrm{ch}}=\\frac{1}{M}\\mathop{\\sum }\\limits_{\\mu =1}^{M}{z}_{i}^{\\mu }{\\eta }_{i}^{\\mu }\\,\\Rightarrow \\,\\langle {C}_{ii}^{\\mathrm{ch}}\\rangle =0\\,,\\quad \\langle {({C}_{ii}^{ch})}^{2}\\rangle =\\frac{1}{M}{\\sigma }^{2},\\end{array}$$<\/p>\n \n (18)\n <\/p>\n and <\/p>\n $$\\langle {C}_{ii}^{{\\rm{c}}}{C}_{ii}^{{\\rm{h}}}\\rangle ={\\sigma }^{2},\\,\\langle {C}_{ii}^{{\\rm{c}}}{C}_{ii}^{\\mathrm{ch}}\\rangle =0,\\,\\langle {C}_{ii}^{{\\rm{h}}}{C}_{ii}^{\\mathrm{ch}}\\rangle =0.$$<\/p>\n \n (19)\n <\/p>\n For the off-diagonal elements, we have <\/p>\n $$\\begin{array}{c}{C}_{ij}^{{\\rm{c}}}=\\frac{1}{M}\\mathop{\\sum }\\limits_{\\mu =1}^{M}{z}_{i}^{\\mu }{z}_{j}^{\\mu }\\,\\Rightarrow \\,\\langle {C}_{ij}^{{\\rm{c}}}\\rangle =\\frac{1}{k}\\,,\\quad \\langle {({C}_{ij}^{{\\rm{c}}})}^{2}\\rangle =\\frac{1}{M}+\\frac{1}{kM}+\\frac{1}{k}\\\\ {C}_{ij}^{{\\rm{h}}}=\\frac{1}{M}\\mathop{\\sum }\\limits_{\\mu =1}^{M}{\\eta }_{i}^{\\mu }{\\eta }_{j}^{\\mu }\\,\\Rightarrow \\,\\langle {C}_{ij}^{{\\rm{h}}}\\rangle =0\\,,\\quad \\langle {({C}_{ij}^{{\\rm{h}}})}^{2}\\rangle =\\frac{1}{M}{\\sigma }^{4}\\\\ {C}_{ij}^{\\mathrm{ch}}=\\frac{1}{M}\\mathop{\\sum }\\limits_{\\mu =1}^{M}{z}_{i}^{\\mu }{\\eta }_{j}^{\\mu }\\,\\Rightarrow \\,\\langle {C}_{ij}^{\\mathrm{ch}}\\rangle =0\\,,\\quad \\langle {({C}_{ij}^{\\mathrm{ch}})}^{2}\\rangle =\\frac{1}{M}{\\sigma }^{2},\\end{array}$$<\/p>\n \n (20)\n <\/p>\n and <\/p>\n $$\\langle {C}_{ij}^{{\\rm{c}}}{C}_{ij}^{{\\rm{h}}}\\rangle =0,\\,\\langle {C}_{ij}^{{\\rm{c}}}{C}_{ij}^{\\mathrm{ch}}\\rangle =0,\\,\\langle {C}_{ij}^{{\\rm{h}}}{C}_{ij}^{\\mathrm{ch}}\\rangle =0.$$<\/p>\n \n (21)\n <\/p>\n Most of the expressions above can be straightforwardly derived by writing down the definition of the sample covariance matrix for a zero-mean variable and then performing the average over neurons. To illustrate this procedure, let us consider one of the most involved terms: <\/p>\n $$\\begin{array}{c}\\langle {({C}_{ij}^{{\\rm{c}}})}^{2}\\rangle =\\frac{1}{{M}^{2}}\\mathop{\\sum }\\limits_{\\mu ,\\nu =1}^{M}\\langle {z}_{i}^{\\mu }{z}_{j}^{\\mu }{z}_{i}^{\\nu }{z}_{j}^{\\nu }\\rangle \\\\ =\\frac{1}{{M}^{2}}\\mathop{\\sum }\\limits_{\\mu =\\,1}^{M}\\langle {({z}_{i}^{\\mu })}^{2}{({z}_{j}^{\\mu })}^{2}\\rangle +\\frac{1}{{M}^{2}}\\langle {z}_{i}^{\\mu }{z}_{j}^{\\mu }\\rangle \\langle {z}_{i}^{\\nu }{z}_{j}^{\\nu }\\rangle \\end{array}$$<\/p>\n \n (22)\n <\/p>\n The probability that zi and zj are part of the same cluster is given, for large N, by \\(\\frac{1}{k}\\). The expression for \\(\\langle {({C}_{ij}^{c})}^{2}\\rangle \\) then becomes <\/p>\n $$\\begin{array}{c}\\langle {({C}_{ij}^{{\\rm{c}}})}^{2}\\rangle =\\frac{1}{k{M}^{2}}\\mathop{\\sum }\\limits_{\\mu =1}^{M}\\langle {({z}_{i}^{\\mu })}^{4}\\rangle +\\frac{1}{{M}^{2}}\\left(1-\\frac{1}{k}\\right)\\mathop{\\sum }\\limits_{\\mu =1}^{M}{\\langle {({z}_{i}^{\\mu })}^{2}\\rangle }^{2}+\\frac{1}{k{M}^{2}}\\mathop{\\sum }\\limits_{\\mu \\ne \\nu }^{M}{\\langle {({z}_{i}^{\\mu })}^{2}\\rangle }^{2}\\\\ \\,=\\,\\frac{3}{kM}+\\frac{1}{M}-\\frac{1}{kM}+\\frac{1}{kM}(M-1)\\\\ \\,=\\,\\frac{1}{M}+\\frac{1}{kM}+\\frac{1}{k}.\\end{array}$$<\/p>\n \n (23)\n <\/p>\n The other terms can be computed following the same steps.<\/p>\n Given that k,\u00a0M are finite, we can approximate the PR for large N as: <\/p>\n $$\\mathrm{PR}\\simeq \\frac{{\\langle {C}_{ii}^{{\\rm{h}}}+{C}_{ii}^{{\\rm{c}}}+2{C}_{ii}^{\\mathrm{ch}}\\rangle }^{2}}{\\langle {({C}_{ij}^{{\\rm{h}}}+{C}_{ij}^{{\\rm{c}}}+2{C}_{ij}^{\\mathrm{ch}})}^{2}\\rangle }.$$<\/p>\n \n (24)\n <\/p>\n Expanding the square and using the results above for the first and second moments of the covariance matrices, we get to our final expression: <\/p>\n $${\\rm{PR}}\\,simeq\\,frac{(1+{\\sigma }^{2})}^{2}\\frac{1}{k}+\\frac{1}{kM}+\\frac{1}{M}{(1+{\\sigma }^{2})}^{2}=\\,M\\frac{k{(1+{\\sigma }^{2})}^{2}}{1+M+k{(1+{\\sigma }^{2})}^{2}}$$<\/p>\n \n (25)\n <\/p>\n From the mathematical expression in equation (25<\/a>), we can see that, in the limit of perfect clusters (\u03c3\u00a0\u2192\u00a00), the function is either limited by the number of rows-neurons (in this case, the k perfect clusters) or columns-conditions M, coherently with the intuition above: <\/p>\n $${\\rm{PR}}\\mathop{\\to }\\limits_{k\\to \\infty }\\,M{\\rm{PR}}\\mathop{\\to }\\limits_{M\\to \\infty }\\,k\\,,$$<\/p>\n \n (26)\n <\/p>\n However, things become more nuanced when k, M and \u03c3 are finite and non-zero. First, as shown in Fig. 5e<\/a>, when k and M are kept fixed, the representation dimensionality decreases with the clustering quality (expressed as the average silhouette score of a population of Gaussian clusters with given k, M and \u03c3). Second, if we fix the quality of clusters and the number of conditions (Fig. 5e<\/a> (left)), we see that dimensionality increases with the number of clusters, with a magnitude that is larger for high silhouette scores (categorical representations). Finally, when fixing the number and quality of clusters, the dimensionality is determined by the number of conditions, with a magnitude that is larger for low silhouette scores (non-categorical representations; Fig. 5e<\/a> (right)).<\/p>\n Different measures of dimensionality<\/p>\n Measuring dimensionality in the presence of noise and determining whether it is high or low can be quite challenging7<\/a>. This is why, in neuroscience, multiple methods are used to assess dimensionality. Each approach is different and, as dimensionality is expressed as a single number, it inevitably emphasizes only specific aspects of the representational geometry.<\/p>\n In our Article, we always refer to the embedding dimensionality of the set of points representing different experimental conditions in the activity space56<\/a>. These points represent patterns of activity recorded at the same time (not trajectories). To discount the dimensions due to the noise, we computed the representation dimensionality of the average positions of the conditions. An alternative approach is to consider a cross-validated measure of representation dimensionality57<\/a>. Notice that separability and other measures related to it6<\/a>, which consider the computational consequences of high dimensionality, are typically insensitive to noise, as they are also cross-validated measures.<\/p>\n Other recent works have introduced dimensionality measures that depend on the spatial and temporal scales considered58<\/a>. For large scales, they get an estimate of the embedding dimensionality and, for short scales, they get an estimate of the intrinsic dimensionality. While these quantities can reveal many other interesting aspects of the representational geometry, here we focused only on the embedding dimensionality because it is the one most relevant to the performance of a linear readout.<\/p>\n Reporting summary<\/p>\n