Introduction

Controlling movement is a critical function of the brain, and a sizable fraction of mammalian neocortex is involved in this task (Currie et al. 2022; Li et al. 2024; Muñoz-Castañeda et al. 2021). Moving in dynamic environments requires processes ranging from trajectory planning and pattern generation to error monitoring and feedback integration. Because of the importance of movement control, a wide variety of approaches have been used to study motor and somatosensory areas: from examining neural biophysics and connectivity, to determining neurons’ direct influence on movement, to relating activity of neurons in relevant subregions to behavior. However, to date the field has not systematically surveyed activity across sensorimotor cortex at single-cell resolution during behavior.

Anatomical mapping has formed the foundation of our understanding of the biological pieces of the brain’s puzzle, and has classically leveraged cytoarchitecture and projection tracing to describe how structure varies across the cortical sheet. These approaches support the view that cortex is divided into many areas that serve different functions, while remaining strongly interconnected (Brodmann 2005; Majka et al. 2021; Yeterian et al. 2012) especially in the mouse (Gămănuţ et al. 2018; Harris et al. 2019; Oh et al. 2014; Winnubst et al. 2019; Zingg et al. 2014). More recent genetic tools have complemented these methods with finer-grained insight into neuron-type distributions and connectivity, providing additional distinctions between parts of the brain (DeNardo et al. 2015; Muñoz-Castañeda et al. 2021; Wang et al. 2020) and revealing multiple parallel circuits within brain areas (Carmona et al. 2024; Park et al. 2022; Ueno et al. 2018). Subpopulations with distinct projection patterns and layer-specific connectivity suggest the existence of multiple functions coexisting within both somatosensory (Kerr et al. 2007) and motor areas (Carmona et al. 2024; Hausmann et al. 2022; Hira et al. 2013b; Muñoz-Castañeda et al. 2021; Zingg et al. 2014). Moreover, strong projections to subcortical, thalamic, and spinal targets originate from subregions of these cortical areas, some of which closely abide by anatomical borders determined by other methods (Carmona et al. 2024; Jeong et al. 2016) while others cross them (Hausmann et al. 2022; Ueno et al. 2018).

To relate activity to behavior more directly, other approaches have instead correlated neural activity with self-generated behavior. Large-scale recordings with simple behaviors have suggested that motor signals are widespread across cortex (Musall et al. 2019; Steinmetz et al. 2019; Stringer and Pachitariu 2019; Wang et al. 2023b), while studies with slightly more challenging behaviors have shown a greater proportion of motor signals in areas closer to the motor or sensory periphery (Li et al. 2024; Wang et al. 2023b). However, the details of what motor signals reside where have not been distinguished, and even in a more challenging motor task neural populations in primary motor cortex (M1) and the primary somatosensory forelimb area yield similarly strong kinematic decoding of many tracked joint angles (Grier et al. 2025). These findings describe a distributed system, where signals are densely shared and where areas’ functional distinctions remain unresolved despite their sharp differences in connectivity with the periphery (Carmona et al. 2024; Ueno et al. 2018). A key question then remains: how are the activity patterns of individual cortical cells organized across sensorimotor cortex during behavior, and how does this relate to previous anatomical mapping?

Here, at single-cell resolution, we densely mapped single neuron activity across mouse sensorimotor cortex during a reach-to-grasp behavior to systematically capture the spatial distribution of tuning properties. We used two-photon calcium imaging to record from over 39,000 neurons as mice performed a 15-target reach-to-grasp task. We focused on five areas involved in this behavior: secondary motor cortex (M2), M1, and the primary sensory areas corresponding to the forelimb (S1-fl), hindlimb (S1-hl), and trunk (S1-tr). We focused on the format of how neurons encoded information, more than correlating activity with the behavior directly. Some tuning properties related to response onsets and peak timing closely followed anatomical borders. Other properties were more closely related to somatotopy, such as tuning sharpness, response duration, or the linearity of responses with respect to target location. Using data-driven approaches, we were also able to identify subpopulations that shared response features, which were salt-and-pepper intermingled within anatomically-defined areas but had characteristic spatial footprints across cortex that spanned areas. Together, these results clarify the organization of the response properties of neural populations at multiple scales, including identifying previously-unknown, overlapping and intermingled area-spanning subpopulations.

Results

Dense sampling of sensorimotor cortex reveals heterogeneous tuning

We implanted cranial windows covering sensorimotor cortex (Fig. 1A-B) in six male mice expressing GCaMP8s in most cortical pyramidal cells (Zhang et al. 2023). Dense sampling of layer 2/3 neurons was achieved using two-photon calcium imaging across 98 fields of view (FOVs) spanning the forelimb, hindlimb, and trunk representations (Fig. 1D). FOVs were aligned to the Allen Mouse Brain Common Coordinate Framework version 3 (Allen CCF) (Wang et al. 2020) using a combination of anatomical landmarks and widefield imaging with paw vibration mapping (Methods; Fig. 1B). A total of 39,398 neurons were extracted across the five relevant regions in the sensorimotor cortex (Fig. 1E). A larger number of neurons were imaged in M2 and M1 (31,453 neurons) than in somatosensory areas (7,945; Fig. 1C).

Figure 1.

Behavior

While neural activity was acquired the mice performed a trained, motorically-challenging, delayed reach-to-grasp-to-drink task (Methods). The task was a modified version of the Galiñanes and Huber (2023) multi-target reaching task. Mice were head-fixed and held two paw rests to initiate a trial (Fig. 2A,B). A trial began with a spout moving into one of 15 positions. After a variable delay period an auditory Go cue (“cue”) was given and a water droplet released, and the mouse could use its right paw to reach toward the target and retrieve its reward. Stereo high speed camera data enabled markerless tracking (DeepLabCut; Mathis et al. (2018)) offline. Animals were successful on 97.1% of trials, and an average of 373 successful trials were collected per session (264-450; 39,229 successful trials total in the dataset).

Figure 2.

As a simple characterization of the behavioral kinematics, single-trial paw centroid trajectories were plotted and colored by their corresponding target position (Fig. 2C). Two relevant features were readily apparent. First, the paw centroid location was strongly target-specific near the spout, demonstrating that reaches were in fact targeted. Second, a looping trajectory near the spout was typically present, which corresponds with repeated reaches to collect all the water from the spout. This was due to a combination of our custom-shaped spouts, which were designed to cause the water to cling, and our use of touch sensors which encouraged strong contact with the spout itself (Methods). Although variability was present in both reaction time (Fig. 2D) and movement time (Fig. 2E), as is typical in reaching data, both intervals were generally fast and indicated that animals were highly engaged in the task. A fuller characterization of how neural activity relates to the kinematic output will be detailed in forthcoming work.

Neural data visualization

To begin investigating the neural activity, we identified which neurons modulated their firing rates during the task. To do so, we computed a modulation score for all cells using a variant of the ZETA method (Methods; Montijn et al. (2021)), which asks whether the activity for a neuron changes over time consistently in relation to any target. Modulated neurons were common and scattered across all five regions of interest in the sensorimotor cortex (Fig. 3A). To quantify their spatial distribution, we computed the fraction of neurons that were modulated across the cortical sheet. The five areas of interest varied in the fraction of neurons that were modulated: M2 had 14%, M1 had 23%, S1-fl had 30%, S1-hl had 25%, and S1-tr had 27% (cue-aligned; Methods). Given that these areas have known spatial organization within them and structure was apparent by eye in the spatial scatterplot of modulated neurons (Fig. 3A), we then produced a high-resolution map of the fraction of neurons that were modulated by location (Fig. 3B). Neurons in two zones were more often modulated than elsewhere: a wide strip along the border between S1 and M1, and an anteromedial patch of M2. The probability of a neuron being task-modulated was highest around the border between M1 and S1-fl, corresponding to the center of the forelimb portion of M1 and the proprioceptive region of S1-fl (Alonso et al. 2023). Probability of modulation was particularly low in a wide strip along the border between M1 and M2.

Figure 3.

Modulation data are shown here computed with neural activity aligned to the Go cue. When aligning data to the paw lift (“Lift”) or spout contact (“Contact”), a mostly-overlapping set of cells were identified as modulated (Fig. 3C, Fig. 3S1). Although more cells were identified as modulated when using the Lift or Contact event alignments, all results in this work were similar across all three possible alignments. Results are shown through the rest of this work aligned to Cue because that alignment best preserved some unique activity structure in M2 (noted below) while preserving all other findings.

We next visualized the tuning profiles of individual neurons through peri-event time histograms (PETHs). Five representative example neurons for each of the five areas of interest are shown (Fig. 3D-H). One striking feature of the neural responses was that neurons were often strongly tuned, with large differences in the activity associated with different targets, despite the fact that the targets were relatively close together in physical space. Several systematic differences between the neurons in each area were also prominent, both in these examples and in our qualitative assessment of the recorded population overall. Neurons in both M2 (Fig. 3D) and M1 (Fig. 3E) were often deeply tuned, with firing rates ranging from no measurable activity for some targets to strong and reliable responses to other targets. M1 neurons typically exhibited briefer responses than M2 neurons, while M2 neurons were often more sharply tuned – that is, strongly active for a smaller number of targets. S1-fl (Fig. 3F) typically exhibited the briefest responses, with a consistent peak of activity for most or all targets. S1-hl (Fig. 3G) and S1-tr (Fig. 3H) neurons had a large variety of response profiles, frequently showing sustained and temporally complex responses and having activity peaks at substantially different times for reaches to different targets. While examples of almost any response profile could be found in each of the recorded areas, these observations reflect systematic differences in the distributions of responses across these areas, quantified below.

Onset of neural activity varied with somatotopy and subregion

To quantify the overall temporal profile of single neurons, we measured when activity for each neuron ramped up. To compute these onset times, we first averaged the activity for each modulated neuron across all trials, creating a single grand-mean trace per neuron, then found the time to half-max. We next asked whether onset times differed by anatomical area (Fig. 4A). These distributions over neurons revealed clear differences in the overall profile of activation: the earliest onsets were present in S1 trunk and hindlimb regions, perhaps due to the animal stabilizing itself; then M2, and finally S1-fl and M1.

Figure 4.

Next, we leveraged our dense spatial sampling to build a map of the time of activation across cortex. The temporal onset metric was computed for each neuron, neurons were binned into pixels, then this map was lightly smoothed (Methods). The resulting temporal onset map revealed a number of sub-area features (Fig. 4B). The earliest activations, occurring during the delay (after the target had moved in but well before the Go cue), were in S1 at the juncture of the trunk, hindlimb, and forelimb regions; and in the portions of M1 and M2 that correspond somatotopically to the hindlimb. As suggested by the area-level result above, these zones presumably reflect the animal preparing its body to support the reach. Shortly after, there was activation in a small zone of S1 that corresponded closely with the tactile subregion of the forelimb region, and in the vibrissal subregion of M2 (see Fig. 1A; Esmaeili et al. (2021); Mayrhofer et al. (2019)). This latter activation might relate to the animal whisking to find the target. Next, activation spread through much of M2, part of the forelimb portion of M1, and into other parts of S1; and then finally to more lateral parts of M1 and M2 that relate to the face and oromanual feeding (Fig. 1A; An et al. (2025)).

Activity profiles of single neurons varied systematically across sensorimotor cortex

To quantify the heterogeneity of tuning profiles across single neurons, we developed several summary metrics that captured the most salient variations in properties that we observed in our inspection of thousands of PETHs. This provides greater interpretability than using unsupervised methods such as Principal Component Analysis (PCA), and is consistent with the practice of feature engineering to extract known-meaningful structure from data where more generic methods fail to capture this structure (e.g., in spikesorting; Caro-Martin et al. (2018), Lee et al. (2021), Lewicki (1994)). Brief intuition is given for each metric below, with complete explanations of each in the Methods.

Response duration

To capture whether a neuron was active only briefly or for a more extended time, we measured the width of the autocorrelation for each neuron. Autocorrelations were averaged across all targets for which the neuron was modulated. Larger values indicate that a neuron remained active for longer. This metric was highest in the posterior regions S1-tr, S1-hl, and the portions of M1 and M2 corresponding to the hindlimb (Fig. 5A).

Figure 5.

Tuning sharpness

To ask whether a neuron was active for only a few targets or many, we computed a metric we refer to as PETH tuning sharpness (Fig. 5B). Neurons attained values near zero if they were similarly active for most targets (i.e., active but broadly tuned), while neurons attained values near unity if they were strongly active for few targets (i.e., sharply tuned). This metric showed that neurons in forelimb-related subregions had sharper tuning than neurons in hindlimb-related regions, with M2 neurons being particularly sharply tuned.

Target tuning linearity

To quantify how linearly a neuron’s activity related to target location in physical space, we correlated the mean activity of the neuron with the targets’ locations (physical distance from an optimized anchor target). The resultant metric map shows that the motor areas were much more linearly tuned with respect to target location than were the somatosensory areas, with M2 even more linear than M1 (Fig. 5C).

Tuning persistence

To capture whether a neuron’s tuning profile was similar across time points or varied at different time points, we computed a tuning persistence metric. Neurons with high persistence maintained the same tuning across longer stretches of time than neurons with low values. Tuning persistence was highest in the anterior of M2, was moderate in M1, and was low in the sensory regions (Fig. 5D).

Peak time variation

To quantify whether a neuron’s firing peaked at the same time for every target or varied by target, we found the peak firing rate of the response to each target, then computed the standard deviation of these peak times across targets. This value is therefore higher if the peak time varied and nearly zero if the timing was consistent. This metric was quite low in M1, indicating highly consistent timing of the activity peak, and was highest in the posteromedial part of M2 (presumably corresponding to the hindlimb representation) and the posterior tip of S1-hl (Fig. 5E).

Local dimensionality

To ask how heterogeneous small neighborhoods of neurons were, we computed a “local dimensionality” metric: for each neuron, we computed the participation ratio (an imperfect but simple eigenspectrum-based measure of dimensionality) of its 30 nearest neighbors. Higher local dimensionality reflects locations where nearby neurons exhibited more heterogeneity in their PETHs with respect to one another. This local dimensionality was highest in a patchy portion of anterior M2, and in anterior S1-fl (possibly corresponding to the tactile region (Alonso et al. (2023), see Fig. 1A) (Fig. 5F). Note that because this metric involved neighborhoods instead of single neurons, it was not included in subsequent analyses.

We also compared our metric maps simply using the top 20 PCs, rotated using VARIMAX to identify a sparser basis (Musall et al. 2019). These maps also exhibited spatial structure, but as expected, when compared with our bespoke PETH features created using knowledge of the problem domain, the spatial structure was weaker (Fig. 5S1).

PETH features changed more strongly at anatomical and somatotopic boundaries

Each of our PETH feature maps varied across the sensorimotor areas, and appeared to change more abruptly in some locations. We wished to ask whether there were in fact relatively discrete transitions in the PETH feature values across cortical space, or whether these values changed only smoothly. To do so, for each PETH feature we computed the gradient: the vector-valued derivative of the PETH feature with respect to two-dimensional cortical space. We then took the magnitude of this gradient computed at every point in the map (Methods; Fig. 6A). Intuitively, this simply asks how quickly the PETH feature is changing at each point on the map. Each map revealed a ridge at the edge of the high-value region (Fig. 6A-E), suggesting that the PETH features maps had edges, albeit soft ones, as opposed to smoothly falling off in value. In the target tuning linearity’s gradient map (Fig. 6A), this boundary was strongest along the border between M1 and S1-fl, cut through the anterior horn of S1-hl (consistent with projection neuron patterns; Ueno et al. (2018) and Yang et al. (2023), and separated what is presumably the forelimb representation from both the hindlimb representation posteromedial and the orofacial representation anterolateral. Other PETH features strongly separated M1 and M2 and somatotopic zones (Fig. 6B,D,E). Each of the candidate functional boundaries was present in at least two of the gradient maps. We therefore pooled the identified boundaries by computing an overall gradient magnitude (the quadratic mean of the individual PETH feature maps; Fig. 6F). This map makes clear that the PETH features of single cells changed most sharply along both area and somatotopic boundaries. Notably, these boundaries did not correspond with ‘edges’ in the map of where more or fewer cells were modulated (Fig. 3B).

Figure 6.

This gradient method has the advantage of producing a smooth map with simple operations. However, to validate the result, we also tested a decoding approach. For each of 10,000 randomly-placed candidate boundaries, we attempted to classify which side of the boundary nearby neurons were on using their PETH features (Fig. 6G; Methods). We expected that classification would be better when the random boundary happened to match places where PETH features were changing more rapidly. The top-performing 2.5% of random boundaries are plotted in Figure 6H. They aligned well with the boundaries revealed by the gradient method.

Areas exhibited distinct and multimodal activity profiles

From the PETH feature maps and distributions (Fig. 5), it was clear that areas exhibited distinct response profiles on average. However, within each area, neurons’ responses were heterogeneous. To better understand the distribution of response profiles within and across areas, we used nonlinear dimensionality reduction (t-distributed Stochastic Neighbor Embedding, t-SNE; Hinton and Roweis (2002)). Specifically, we summarized each neuron’s response with a five-dimensional feature vector comprising the five PETH features described above, then applied t-SNE to these vectors for the full neural population spanning all five brain areas (Fig. 7A). Finally, we found where the neurons from each area were embedded in the t-SNE space (Fig. 7B). The resulting two-dimensional t-SNE projection led to three findings. First, single neurons from each region were found across most of the t-SNE space (Fig. 7B). This makes clear that much of the diversity of neuron activity patterns was observed in each area. Second, when looking at where neurons in each area most commonly appeared in the t-SNE space, most areas had peaks centered at different locations – though with overlapping distributions (Fig. 7C). The largest differences were between motor and somatosensory areas, with S1-fl bridging the other S1 subregions and M1. This implies that areas had distinct distributions of their PETH features, but also that there was an interpretable spectrum of typical properties from M2 to M1 to forelimb S1 to hindlimb and trunk S1. Third, the distributions of PETH features within-area were strongly multimodal for several areas. This suggests that the heterogeneity of neurons within each area may form loose response ‘types,’ especially in M2, M1, and S1-fl. This multimodality motivated further analysis targeted at understanding the range of neural response profiles in each area, and whether neurons with specific response profiles were preferentially found in localized regions of each cortical area.

Figure 7.

Neighboring sensorimotor regions shared PETH features in complex spatial patterns

The multimodality of feature vector distributions within-area opens a range of possibilities for how PETH response features could be organized across the cortex, both within and across areas. To work with these multimodal distributions, we fit a generative model to the PETH feature vectors defined above for each area’s population of neurons. This approach modeled the feature vectors of an area’s neurons as arising from a Gaussian Mixture Model (GMM). Doing so explicitly models different modes as arising from different Gaussian distributions, and retains probabilistic interpretability. In agreement with the t-SNE results, M1, M2, and S1-fl were best modeled by the GMM as having multiple components (Methods): four for M2, three for M1, and two for S1-fl, while S1-hl and S1-tr were best modeled as being unimodal. The Gaussians fitting different modes were well-separated, as expected when fitting a multimodal distribution (Fig. 8S1).

Figure 8.

We hypothesized that even though the distributions of response profiles were multimodal within-area, different areas might still be well separated. To test this, we computed a data likelihood for each neuron – p(neuron i’s PETH features | GMM_area) – as arising from each area’s distribution (Fig. 8A). Intuitively, this produces a measure of how similar each neuron is to the population of neurons in an anatomically-defined area. We then binned and smoothed over neurons by cortical location to produce the maps in Figure 8B-F.

When computing this likelihood map using neurons in M2 to define the distribution (Fig. 8B), the highest likelihood zone was, as expected, in M2. However, we found that a patch of M1 was strongly similar. If we instead defined the distribution using M1 neurons (Fig. 8C), again the highest-likelihood zone was within the defining area. However, two additional features were salient. First, the anterior portion of S1-fl was also strongly similar to M1, in line with previous results (Grier et al. 2025). Second, the posteromedial portion of M1 differed strongly from the rest of M1, consistent with a difference between forelimb and hindlimb somatotopic regions (Maurer et al. 2023). Considering other areas, defining the distribution with S1-fl led to high likelihoods in adjacent M1 and S1-hl (Fig. 8D), while defining the distribution with S1-hl or S1-tr produced similar spatial layouts and bled over mostly into the other sensory areas (Fig. 8E,F). These results make clear that the working hypothesis – of areas with well-separated feature distributions – is incorrect. Instead, neurons in portions of some areas were similar to neurons in portions of other areas, and two dissimilar areas could both be similar to the same third location in cortex. These neural populations therefore warranted further dissection.

Subpopulations spatially overlapped and spanned anatomical areas

The likelihood maps above suggested that populations of neurons with similar response profiles were not always confined to the anatomically-defined areas. For example, portions of M1 and S1-fl were strongly similar (Fig. 8BC,D). This hinted at a second hypothesis: that some of the ‘modes’ (groups of neurons) discovered separately in each area might correspond. If this were due to a subpopulation that was shared across the areas, we would expect that the M1 GMM would contain a component (a one-area estimate of the subpopulation) with similar characteristics to a component in the S1-fl GMM (Fig. 9A). We therefore computed a distance metric (Methods) between each pair of GMM components (Fig. 9B). This analysis revealed that many components identified independently in different areas were indeed highly similar, and formed four clusters. We therefore combined components within each of the four clusters (Methods) and took the neurons belonging to each cluster as the subpopulations.

Figure 9.

These four subpopulations differed in two major ways. First, by construction, different subpopulations were characterized by distinct PETH features (Fig. 9G). Second, their spatial likelihood maps differed strongly as well (Fig. 9C-F). One subpopulation, which we refer to as the “Anterior” subpopulation, was almost unique to M2 (Fig. 9C). It had uniquely high tuning persistence and tuning sharpness along with highly linear target tuning, arguing for a relatively abstract and external target representation (Galiñanes and Huber 2023). A second subpopulation, the “Forelimb motor” subpopulation, spanned M2 and the forelimb portion of M1 (Fig. 9D). These neurons were also strongly linear in their target tuning, but less sharply tuned and with tuning that was somewhat less persistent (i.e., changed more over the course of the movement). The “Forelimb somatomotor” subpopulation spanned M1 forelimb and orofacial zones and the anterior portion of S1-fl (Fig. 9E). This subpopulation was more broadly tuned, was less linearly tuned to target location, changed tuning over the course of the movement (low tuning persistence), and responded very briefly (low response duration). Nevertheless, each neuron in this subpopulation still responded at its own characteristic time in the movement (low peak time variation). Finally, a “Hindlimb somatomotor” subpopulation (Fig. 9F) was centered on S1-hl and spilled into S1-tr, the posterior portion of S1-fl, and the posterior tips of M1 and M2. This subpopulation was similar to the Forelimb somatomotor subpopulation, but each neuron no longer fired at a specific time point in the movement; instead, neurons fired at different times for different targets (high peak time variability).

To confirm that this identification of area-spanning subpopulations was robust, we validated these findings with several other methods. To determine whether our finding of multimodality (i.e., subpopulations) was due to our use of hand-designed PETH features, we repeated some of our analyses using the top 20 PCs of the PETHs instead of the PETH feature vectors. Using PCA as a preprocessing step in this way is more common in the field. Performing t-SNE on the PCs, analogously to Figure 7, again yielded multimodality (Fig. 9S1A). Similarly, repeating the analysis of Figure 8 on the PCs produced results that were consistent with our feature-based likelihood maps, just marginally less distinct (Fig. 9S1C-G). This indicates that although our choice of PETH features enabled a sharper view of the distinctions between subpopulations, these features were not necessary to detecting them.

To validate the clustering, we also fit a GMM to all neurons without consideration of their area membership. This process produced three components with maps that were nearly identical to three of the four subpopulations identified with the area-aware method (Fig. 9S2). The only noticeable difference was that the Anterior subpopulation was lost when all areas were fitted together, likely due to having a lower membership than the other subpopulations.

To determine whether our distinct subpopulations required recordings from a large area, we directly examined a smaller patch of M1 (480 µm diameter) in an overlap zone (Fig. 9S1B, left). Applying t-SNE to either our PETH features (Fig. 9S1B, top right) or to the top 20 PCs of the PETHs (Fig. 9S1B, bottom right) again yielded multimodality. This demonstrates that even with a single ordinary-size field of view, distinct subpopulations in overlapping zones could be readily detected in the context of our behavior.

Notably, the high-likelihood regions for different subpopulations spatially overlapped. The Anterior and Forelimb motor subpopulations overlapped throughout M2; the Forelimb motor and the Forelimb somatomotor subpopulations overlapped in the forelimb portion of M1; and the Forelimb somatomotor and Hindlimb somatomotor subpopulations overlapped in a band along the S1-fl and S1-hl border. This indicates that although several of the PETH feature boundaries found in Figure 6 aligned with area borders, several of these boundaries corresponded to the edges of subpopulations that spanned other borders. To confirm this interpretation, we computed the gradients of the likelihood maps directly (Fig. 9S3). These maps closely matched those computed from the PETH feature gradients directly (Fig. 6), as expected.

Spatial overlaps between subpopulations indicate intermingled members

The spatial overlap of subpopulations strongly suggests that members of different subpopulations were salt-and-pepper intermingled. However, the use of likelihoods above mixes two properties: the fraction of neurons in a given neighborhood that are members of each subpopulation, and how similar each neuron in a subpopulation is to a prototypical member of that subpopulation. If these properties varied systematically within the zone covered by a subpopulation, additional spatial structure could be present.

To search for such additional spatial structure, we first classified each neuron into its maximum-likelihood subpopulation, then mapped the subpopulation identity of neurons (Fig. 10A). As expected, members of the different subpopulations were locally salt-and-pepper intermingled, including in overlap zones. Also consistent with the likelihood maps, membership prevalence varied by area (Fig. 10B). Next, we asked whether the prevalence of membership in each subpopulation followed the likelihood maps. Indeed, the map of subpopulation membership probability matched the likelihood maps closely (Fig. 10D). Finally, if different parts of a subpopulation varied substantially in their PETH features, a model built on neurons from one end of the spatial zone (e.g., the M2 end) would assign lower likelihoods to neurons at the other end (e.g., the M1 end), and vice versa (Fig. 10C). To determine if this was the case, we produced likelihood maps using the GMM components computed from each area separately, thresholded the likelihoods to produce a contour, and compared their shapes. Subpopulations had nearly identical spatial footprints no matter which area was used to identify them (Fig. 10E). This was true even in cases where only a relative handful of members of that subpopulation were present in the area: for example, the Forelimb somatomotor and Hindlimb somatomotor subpopulations could be identified using M2 data, even though there was no high-likelihood zone for these subpopulations within M2. Note that this consistency also validates the use of clustering to combine components and identify the subpopulations in the first place. Together, these results argue that the overlap zones for different subpopulations are due to locally intermingled members, and the PETH properties within the subpopulation were consistent across the territory covered by that subpopulation.

Figure 10.

Discussion

By densely sampling layer 2/3 of mouse sensorimotor cortex during a 15-target reach-to-grasp behavior, we mapped how single cell response characteristics were distributed across five areas. Neurons responded heterogeneously within all imaged areas, and within M2 and M1 in particular, but the average response properties of areas were distinct. In parallel, we identified four subpopulations of neurons with unique activity profiles, which spanned multiple areas and were spatially intermingled in the zones where different subpopulations overlapped. One subpopulation was unique to M2, one spanned M2 and the forelimb portion of M1, one was centered in S1-fl but extended into S1-hl and the forelimb and orofacial portions of M1, and one spanned the likely hindlimb representations of S1, M1, and M2. Members of these subpopulations did not vary substantially with their spatial location, suggesting that the average differences observed between areas were largely, though probably not exclusively, a function of their containing different proportions of the subpopulations. These activity-based maps thus reveal new structure in sensorimotor cortex and complement a variety of existing mapping methods.

The investigation of coding here was somewhat different than the traditional approach. We did not focus on traditional questions of “what” variables are encoded, such as external target location vs. kinematics or reach direction vs. reach duration. Instead, we focused on “how” task-related activity was structured, such as response duration, timing of peaks across targets and tuning linearity. Perhaps surprisingly, this generally yielded the same answer as an unsupervised approach to coding, implemented here by using PCA on the PETHs. It also makes clear, though, that there are potentially important differences not just in the content of the activity, but also its format. In particular, we found that anterior areas tended to have responses that were sharply tuned for target but broader in time, while posterior areas tended to be broadly tuned for target but have very brief responses that consistently peaked at the same time point in the movement.

What, then, might be the functional roles of these subpopulations? While future work relating neural responses to behavior will be essential to answering this question, the present analyses suggest some possibilities. M2 contained neurons that were sharply and linearly tuned for target location, as well as having earlier onsets, longer responses, and consistent tuning over the movement. This might suggest a more abstract (Bernardi et al. 2020) or external encoding (Galiñanes and Huber 2023). S1-hl and S1-tr also had early responses, with activity spanning the movement but differing in complex ways across targets. This may suggest activity related to bracing or stabilizing the animal (Disse et al. 2023; Nandakumar et al. 2021). The responses of S1-fl were briefest and typically weakly tuned for target, but peaked at a tightly-timed moment in the movement for most or all targets, with the peak timing being neuron-specific. These responses may correspond to some specific sensory signal like afferent feedback during paw lift (Grier et al. 2025; Mathis et al. 2017), but understanding what that signal would be will require future study. The shape of the Forelimb motor subpopulation is similar to the shape of a recently proposed organization based on corticofugal projection patterns, the Medial Anterior Cortex (MAC) (Yang et al. 2023). This would be consistent with MAC’s proposed role in reaching. Finally, M1 was at the center of the overlap between subpopulations. This is reminiscent of its suggested role in taking top-down commands from M2 and elsewhere and bottom-up sensory input from S1 to actively control the limb (Galiñanes and Huber 2023; Heindorf et al. 2018; Terada et al. 2022).

The area-spanning nature of the subpopulations could arise from a number of possible organizations of neural circuits. One hypothesis is that each subpopulation represents a single recurrent circuit that is distributed across areas. In this hypothesis multiple circuits exist in parallel within each area, each strongly connected within the subpopulation and less strongly connected with the other subpopulations. An alternative hypothesis is that a single area contains the circuit that generates the signals characteristic of a single subpopulation, then this activity is broadcast to other areas through the strong connectivity of sensorimotor cortex (Muñoz-Castañeda et al. 2021; Zingg et al. 2014). This broadcast would then cause that characteristic activity pattern to appear in the recipient areas. Recording neurons that project from one area to another, and recording from other layers corresponding to later stages of processing and outputs, will be critical to distinguishing these hypotheses and more fully understanding these circuits.

The maps produced here, using response properties of neurons during task performance, aligned with anatomical and stimulation-based maps in some respects but differed in others. We identified response-property boundaries that aligned closely with the anatomical borders between areas M2, M1, and S1. We also identified response-property boundaries consistent with somatotopy, separating forelimb from hindlimb and forelimb from orofacial representations. Interestingly, these somatotopic boundaries were generally stronger than those corresponding to “true” area borders, despite the known coarseness of somatotopic maps in primate motor areas (Mitz and Wise 1987; Neafsey et al. 1983; Qi et al. 2010; Van Acker et al. 2013); whether the maps are similarly coarse in mouse remains to be tested. Some of our maps also distinguished the tactile portion of S1-fl from the remainder of the area, and grouped the proprioception-receiving portions of S1-fl and M1 (Alonso et al. 2023). We could not, however, readily identify boundaries between other known subareas of M2, such as the Rostral Forelimb Area (RFA; Li et al. (2024); Morandell and Huber (2017); Tennant et al. (2011)), Anterolateral Motor area (ALM; Chen et al. (2017); Inagaki et al. (2018); Li et al. (2016)), or the LAC/MAC distinction (Yang et al. 2023). It is possible that these distinctions would appear with other tasks, which might require different structure in the coordination between areas, or that they would be apparent in deeper layers.

Crafting additional PETH features, or using end-to-end neural network approaches to discover other features, might enable the discovery of additional structure (Minderer et al. 2019; Wang et al. 2023b). For example, our PETH features were chosen to be invariant to the onset time of activity, but these onset times were markedly later in lateral M1 than in adjacent M2 or S1-fl. Including onset times or other features would presumably result in finer and finer subdivisions of sensorimotor cortex. This suggests that although our current methods find a large amount of structure in sensorimotor cortex, more structure in the responses surely exists.

Our findings of larger-scale structure, in the form of area-spanning subpopulations, may intersect with three prior findings. First, several widefield calcium imaging studies have also found that S1-fl activity groups with activity in portions of M1 and S1-hl (Costa et al. 2025; Musall et al. 2019; Saxena et al. 2020; Zatka-Haas et al. 2021). Relatedly, orofacial representations share similar relationships with behavior across M2, M1, and S1 (An et al. 2025; Musall et al. 2019); and both these spatial patterns of coactivation reflect anatomical connectivity (Zingg et al. 2014). These previous findings of larger-scale grouping are consistent with our results. Second, in several cortical areas, populations with different projection targets have known functional differences (Chen et al. 2013; Currie et al. 2022; Hwang et al. 2019; Kim et al. 2018). However, these specific populations reside mainly in layer 5, while imaging here was performed more superficially. Third, the existence of subpopulations of neurons in sensorimotor cortices have been found independently by using genetic identification of excitatory neuron subtypes, and these subpopulations exhibit distinct activity patterns (Hira et al. 2013a; Li et al. 2024; Mohan et al. 2023; Musall et al. 2023). Whether those subpopulations correspond to the ones found here cannot yet be reconciled. However, many of the genetically-identified neuron subtypes in those studies are present only in deeper cortical layers, while the imaging performed here was exclusively superficial. This suggests that the subpopulations we have identified are unlikely to correspond closely with the particular subtypes found thus far via genetic tools.

The distributions of PETH properties we observed were multimodal, but these distributions were continuous; we did not observe tight clusters in feature space separated by gaps. In interpreting the observed multimodal distribution, we note that a uniform distribution can be transformed into a multimodal one, or vice versa, by a continuous nonlinear function. This is a fundamental problem whenever searching for structure in the distribution of neural activity: even examining the raw spike rate of a neuron involves the spike threshold nonlinearity, and what property counts as the ‘fundamental, linear’ one is almost a philosophical question. Nonetheless, the potential for nonlinear measurement makes it challenging to determine whether there are true “categories” of neurons (Kaufman et al. 2022; Posani et al. 2025). If the interpretation of separate modes as subpopulations were incorrect under some relevant nonlinear transformation, our results would nevertheless imply that M2 has the highest diversity of responses, with shrinking variety toward the posterior as well as changes in the average response profile. If this alternate interpretation were the strongest conclusion possible, it would still be essential to understanding the layout of function across sensorimotor cortex. Confirming that these modes are more deeply meaningful will require future external validation.

Several limitations of the present work should be noted. First, the data here include only superficial layer 2/3. Activity patterns in other layers, or even in the deeper portion of layer 2/3, may differ (Currie et al. 2022; Heindorf et al. 2018; Masamizu et al. 2014). Second, while we took great care in aligning our FOVs to the Allen CCF, we did not analyze the anatomy of fixed tissue for each mouse. Limitations in our alignment, differences between the anatomy of the strain of mice used here from the strain used for building the atlas, and possible mouse-to-mouse variation in the borders of areas all introduce jitter. The sharpness of the boundaries we identified therefore serves as a lower bound. Third, although we achieved coverage of a large swath of sensorimotor cortex, we sampled motor areas more densely than sensory and not every area was densely sampled in every mouse. We therefore may have underestimated heterogeneity in S1, which was sampled less. Fourth, although we obtained high-speed video of the behavior, this work does not consider the kinematics. The relationship of these neural data to the detailed kinematics will be treated at length in future work. Fifth, this work considered only trial-averaged responses, not single trials. This choice enabled us to compare much larger areas of cortex by pooling across sessions, but precluded many types of analysis including of noise correlations (Kiani et al. 2015; Kohn et al. 2016; Zohary et al. 1994). Future work could leverage simultaneous recording of multiple areas to investigate single-trial co-fluctuations in neural activity within subpopulations and across areas.

In summary, these findings provide a response-based atlas of sensorimotor cortex and reveal overlapping subpopulations with intermingled members. The nature of the circuits that underlie these subpopulations, and the content of communication between areas more broadly, remain questions of fundamental importance for understanding how movement is generated and feedback incorporated.

Methods

Subjects and surgical procedures

All procedures were approved by the University of Chicago Institutional Animal Care and Use Committee. For two-photon imaging, six male mice aged 8-12 weeks were used. Mice were crosses of Slc17a7-IRES2-Cre-D (JAX strain 037512; Harris et al. (2014)) and TIGRE2-jGCaMP8s-IRES-tTA2 (JAX strain 037719; (Sweeney et al. 2025)) which express GCaMP8s in most glutamatergic neurons (Wang et al. 2023a). Each animal underwent a single surgery, then were individually housed on a reverse 12-h light/dark cycle, with an ambient temperature of 71.5 °F and a humidity of 58%. Experiments were conducted in the afternoon, during the animal’s dark cycle.

Anesthesia and basic surgical preparation were as described in Grier et al. (2025). A circular craniotomy was outlined with a 5-mm biopsy punch and completed with a dental drill, centered over the left forelimb M1 at approximately 0.4 mm anterior and 1.5 mm lateral of bregma. The craniotomy was cleaned with SurgiFoam (Ethicon) soaked in phosphate-buffered solution (PBS). If the dura remained intact, a durotomy was performed with a needle.

A retrograde virus carrying a tdTomato construct (AAVretro-tdTomato, diluted to 1 × 10¹³ particles per mL in PBS, Addgene stock 59462-AAVrg) was injected in S1-fl (mice 1-3; 0.0 mm A, 2.65 mm L relative to Bregma) or the center of the forelimb representation of M1 (mice 4-6; 0.4 mm A, 1.5 mm L). 100 nL of virus was pressure injected (NanoJect III, Drummond Scientific) at a depth of 300 µm from the cortical surface over 5 min. The injection pipette was kept in place for another 3-6 minutes to ensure viral dispersion before being removed slowly over 2-3 minutes. This injection retrogradely labeled cell bodies in other areas that sent projections to the injected area. In this paper the labeling was used solely for stabilizing the imaging plane (see below) throughout the cortex.

The craniotomy was sealed with a 5 mm #1 round cover glass (Thomas Scientific) glued in place first with tissue glue (VetBond, 3M) and then with cyanoacrylate glue (Krazy Glue) mixed with black dental acrylic powder (Ortho Jet; Lang Dental). In all mice two layers of MetaBond (Parkell) were applied, and then a custom water jet-cut stainless steel head bar was affixed to the skull with black dental acrylic. Animals were awoken by administering atipamezole via intraperitoneal injection and allowed to recover at least 3 days before water restriction.

Behavior and training

The behavioral task (Fig. 2A) was a variant of the multi-target water-reaching task of Galiñanes and Huber (2023). Water-restricted mice were head-fixed with their forepaws on paw rests (bent partially threaded screw) and the hindpaws and body supported by a custom 3D-printed translucent seat. Trials were triggered after the animal held both paw rests for 100 ms. After an additional 100 ms hold period a water spout (22 gauge, modified 90-degree bent, 1-in blunt dispensing needle, McMaster) was moved into one of 15 positions in front of the mouse’s face via a set of three motors (LSM050B-T4A, Zaber). The spout was customized, with ∼2 mm of the side facing away from the animal ground down; this caused the water to cling to the last few mm of the spout and not just hang from the tip. The animals were then required to hold the paw rests continuously for an additional 400-600 ms delay period to trigger the delivery of a water droplet on the spout. If the animals released the paw rests during this time the delay was reset; this occurred at variable rates across animals (20-85% of trials). Animals were allowed a maximum of 2000 ms to achieve paw rest contact for the entire delay; if this did not occur the spout was removed to the neutral location and a 5000 ms time-out was triggered (<3% of all trials across animals). Once the paw rests were held continuously for the variable delay time, a 4000 Hz tone was played by stereo speakers as the Go cue and a 2- to 3-μl droplet of water appeared on the spout. The cue tone lasted 500 ms or until the mouse made contact with the water spout, whichever came first. The mouse could grab the water droplet and bring it to its mouth to drink anytime after the tone began. The spout was removed after 4000 ms, triggering a variable ITI of 3000-7000 ms, and automatically wiped against a small sponge to remove any remaining water. Both the paw rests and spouts were wired with capacitive touch sensors (Teensy 3.2, PJRC). Control software was custom-written in MATLAB R2018a using PsychToolbox 3.0.14. Touch event monitoring and task control were performed at 60 Hz using the Teensy.

The mice were trained to make all reaches with the right paw and to keep the left paw on the paw rest during reaching. For the first 1-3 days, mice were head-fixed and provided water directly to their mouths to acclimate them to head-fixation and reward collection. During this period, the animals were encouraged to make consistent contact with both paw rests to trigger water delivery to their mouths. After acclimation and consistent trial initiation, the spout was withdrawn to a position on the right side of their face intended to elicit forelimb reaching with the right arm. If animals were reluctant to reach, small droplets of water were placed on their right whisker pad with a blunt syringe to induce grooming and water collection. After animals reliably performed reaches toward the rightward spout location (80-90% success), its position was moved medially and anteriorly towards the midline. Once mice could reliably reach towards this new central spout position, two new target locations were introduced on either side. Once in this three-target training stage, each target was cued in blocks. The size of the blocks was reduced as the mice progressed in their training. After 80-90% success at the three target configuration, two additional targets were added on the left and right sides, creating a five target semicircle. The targets were then progressively spaced apart along the semicircle until they reached their final azimuthal spacing. Variation in the Z position was then introduced and gradually increased until the final Z locations were achieved and the animal could reliably reach towards all 15 target positions. This training procedure took 2-3 weeks, although the behavior continued to solidify for at least two more weeks. Data presented here were collected after 4-6 weeks’ experience with the task.

Two-photon calcium imaging

Calcium imaging procedures were similar to those described previously (Grier et al. 2025). Briefly, imaging was performed with a Neurolabware 2p microscope running Scanbox 4.1 and a pulsed Ti:sapphire laser (Vision II, Coherent). Depth stability of the imaging plane was maintained using a custom plugin that made automatic movements of the objective up to every ∼10s based on comparing red-channel data with an initially-acquired image stack. Offline, images were run through Suite2p to perform motion correction, region-of-interest (ROI) detection, and fluorescence extraction. ROIs were manually curated using the Suite2p GUI. Fluorescence was neuropil-corrected using a coefficient optimized by minimizing the mutual information between the corrected trace and the neuropil trace (Allen 2025). As previously described (Grier et al. 2025; Zhu et al. 2022), fluorescence was then detrended and passed through OASIS using the ‘thresholded’ method, AR1 event model, and limiting the tau parameter to be between 300 and 800 ms. To put events on a more useful scaling, for each ROI, we found the distribution of event sizes, smoothed the distribution (ksdensity in MATLAB, with an Epanechnikov kernel and log transform), found the peak of the smoothed distribution, and divided all event sizes by this value. This rescales the peak of the event size distribution to have a value of unity.

As previously described, when time-locking the data for visualization, modulation testing, or further analyses, the deconvolved events were resampled into 10-ms bins. This resampling was performed by assigning a fraction of each event into the new 10-ms bins proportionally to how much the 10-ms bins overlapped the original ∼32.2 ms frames, and taking into account the exact time each ROI centroid was sampled based on its position within the field of view (FOV).

All cells were imaged between 140 and 210 µm below the cortical surface, in superficial layer 2/3. We chose this depth range both because it yielded large numbers of modulated neurons, and because GCaMP expression in these animals was lower in deeper layer 2/3 of S1. Imaging in this limited depth range therefore produced more comparable data across areas.

Measuring modulation

To determine whether a neuron was modulated in the task, we used an adapted ZETA method (Montijn et al. 2021) as described in Grier et al. (2025). Briefly, trial-averaged responses were compared to shuffled trial-averages computed by circularly permuting each trial by a random offset. Modulation was computed as the maximum difference between the real and shuffled data. This procedure was performed separately for each target and for locking data to the cue (−300 ms to +300 ms), lift (−200 ms to +400 ms), and first contact (−300 ms to 300 ms). For each locking event (cue, lift, and first contact), neurons with p-values ≤ 0.01 were considered modulated. For all analyses, only neurons modulated by the relevant locking event were included.

Mapping cell positions onto a common coordinate frame (Allen atlas)

Each 2P FOV was aligned based on the vasculature to an image of the entire window captured via widefield imaging. S1-fl was identified by paw vibration (Alonso et al. 2023), with imaging and preprocessing performed exactly as described in Grier et al. (2025) except that the camera was upgraded to a pco.edge 5.5 CLHS sCMOS camera (Excelitas). The vibration map was aligned to the widefield vasculature image via shared vasculature landmarks. This composite image was then scaled and translated to best align the S1-fl activation shape to its CCF boundary, while aligning the central sinus vasculature to the midline overlapping bregma. See Figure 1S1 for vibration map to Allen CCF alignment across animals. 2P images from the cortical surface from each session were then scaled to 720 µm width by 980 µm height and aligned to the registered vasculature via rigid transformations. This series of transformations produced global coordinates for each ROI from its location within the 2P FOV.

Summarizing subregions from prior literature

To produce the map in Figure 1A, outlines of subregions described by prior literature were extracted. Our review of the literature is non-exhaustive and serves simply to illustrate the heterogeneity of function that has been described in previous studies. Corresponding figures were overlaid with the Allen CCF and rigidly transformed to qualitatively match anatomical boundaries as closely as possible. Outlines were then extracted by tracing the subregion boundaries by hand.

From anterior to posterior we describe the acronyms and the corresponding citations as follows: Secondary motor cortex (M2), or MOs in the Allen CCF. Anterior lateral motor area (ALM), as originally described in Komiyama et al. (2010) and further delineated in Guo et al. (Guo et al. 2015) and Chen et al. (2017). Rostral Forelimb Area (RFA) as described in Tennant et al. (2011) and further delineated in Carmona et al. (2024). Rostral Forelimb Orofacial area (RFO) as described in An et al. (2024), and referred to as tjM1 in Mayrhofer et al. (2019). Secondary vibrissal motor area (M2-vb) as described in Mayrhofer et al 2019 as the frontal region activated by C2 vibrissal stimulation, and called medial motor area or MM in Chen et al. (2017). Primary motor cortex (M1) or MOp in Allen CCF. Caudal Forelimb Area (CFA) as described in Tennant et al. (2011) via stimulation and Munoz-Casteneda et al. (2021) via tracing. Primary hindlimb motor area (M1-hl) as traced from hindlimb musculature by Maurer et al. (2023), and sensory cortical areas by Zingg et al. (2014). Primary forelimb somatosensory area (S1-fl), or SSp-ul in Allen CCF. Primary proprioceptive and tactile forelimb somatosensory areas (S1-pro and S1-tac) as originally identified in Alonso et al. (2023). Primary hindlimb somatosensory area (S1-hl), or SSp-ll in Allen CCF. Primary trunk somatosensory area (S1-tr), or SSp-tr in Allen CCF. Primary mouth somatosensory area (S1-m), or SSp-m in Allen CCF. Primary nose somatosensory area (S1-n), or SSp-n in Allen CCF. Primary vibrissal somatosensory area (S1-bfd) or SSp-bfd in Allen CCF. Anterior visual cortex (A) or VISa in Allen CCF. Retrosplenial cortex (RSC), or RSPagl, RSPd and RSPv in Allen CCF. Cortical cervical spinal cord (CSC) and lateral rostral medulla (latRM) projections as traced in Yang et al. (2023).

Computing map images from single-cell properties

To quantify and visualize the spatial distribution of properties of imaged neurons, we used a binning procedure to prevent biases due to uneven sampling. We binned the neurons into 70 µm-wide square pixels, then averaged the property of interest over all neurons within the pixel. Pixels that contained zero neurons were excluded from further analysis. The images were then Gaussian smoothed (s.d. 150 µm). To Gaussian-smooth the images while excluding empty pixels and accommodating irregular map edges, the final normalization for each pixel was the sum of the weights of contributing pixels, excluding empty pixels. For modulation images, each neuron was assigned a zero if not statistically modulated or a one if modulated. For PETH feature maps, each neuron was assigned its PETH feature value. The one-dimensional pixel values were then mapped onto a perceptually-linear colormap (Kovesi 2015).

PETH feature computations

Tuning sharpness

Tuning sharpness measures how selective each neuron is to some targets over others. For each neuron, we obtained the peak response v(c) for each condition (i.e., target) c. The tuning sharpness then computes the average response as a fraction of the max response, and subtracts this from one

This metric ranges from near one, meaning a neuron is highly selective for a single condition, to near zero, indicating a neuron achieves a similar peak firing rate for all targets.

Target tuning linearity

To assess how linearly related a neuron’s response profile was to the target locations in physical space, we performed an analysis related to Representational Similarity Analysis (RSA; Kriegeskorte et al. (2008)) but limited it to linearity with respect to an anchor target. To do so, we first computed a matrix containing the Euclidean distances in real-world coordinates between all pairs of targets, and a matrix of the ‘distances’ between the firing rates for each pair of target responses (averaged over the time window from −200 ms prior to the cue to 400 ms post-cue). We then computed the pairwise correlation (R²) between row m of the target distance matrix and the corresponding row m of the neural distance matrix, for all rows m. Finally, we took as the target tuning linearity the highest of the 15 R² values. Intuitively, this looks for a linear relationship between firing rate and target distance with respect to a ‘preferred’ target. Values near one indicate a neuron whose firing scaled perfectly linearly with target location assessed in this way.

Response duration

We quantified the response duration of single neurons by analyzing the autocorrelation structure of their responses over time. For each neuron, we first computed the autocorrelation function r(τ) of the trial-averaged response over time x(t) to each target separately:

where τ represents time lag, x̄ and is the mean of the target-specific response over time. Only targets for which the neuron was statistically significantly modulated were included. From this autocorrelation, we determined the full width at half maximum (FWHM). The mean FWHM was calculated across all the targets for which the neuron was modulated (up to 15), resulting in a scalar measure of response duration for each neuron. Higher values of mean FWHM indicate neuronal responses that are more sustained (with a broader peak), whereas lower values indicate a more transient response.

Peak time variation

Some neurons exhibited firing rate peaks that were very similarly timed for all targets, while other neurons did not. To quantify this, each neuron’s PETH was first computed as described above, then for each condition (target) c, the peak response time was identified as the time bin at which the trial-averaged trace reached its maximum amplitude. To prevent the inclusion of conditions with minimal modulation, we restricted this analysis to conditions whose peak amplitudes were at least 25% of the maximum peak amplitude across all 15 conditions for each neuron. The variability in peak timing was then quantified as the standard deviation of the selected peak times:

where n is the number of conditions meeting our inclusion criterion and t̄^peak is the mean peak time across conditions above the threshold. A higher σ_peak indicates greater variability in the timing of peak responses across targets.

Tuning persistence

For each neuron, we first computed the PETH as described above. Then, for each neuron i, we identified the target for which the neuron achieved the greatest peak response, then found the time point t_rmax of that peak. An anchor tuning vector v_i(t_rmax) was computed, whose entries are the PETH amplitudes across targets at time t_rmax. At every time point t, we measured the similarity of the instantaneous tuning profile vector v_i(t) to the anchor vector by the Pearson correlation

This yields a “tuning-distance” curve r_i(t) that falls off as the neuron’s PETH tuning profile deviates from that at t - t_rmax. To summarize these correlation plots into a single value, we computed the area under the tuning-distance curve

using trapezoidal integration over the sampled timepoints. Neurons with slowly decaying correlations (large AUC) achieved higher tuning persistence, whereas rapidly decaying neurons (small AUC) received low tuning persistence.

Local dimensionality

We estimated the effective dimensionality of local population activity by computing a common and simple metric, the participation ratio (PR), in the spatial neighborhood of each recorded neuron. This metric is limited to only considering the singular value spectrum of the local activity, but when applied exclusively to modulated neurons it provides at least a coarse idea of the variety of responses. For each neuron, we identified its K = 30 nearest neighbors (including the anchor neuron) by Euclidean distance in anatomical (x, y) coordinates. For each neighbor, the trial-averaged PETH was vectorized; stacking these vectors formed a data matrix X whose columns correspond to neurons. We then computed the eigenvalues λ of the covariance matrix of X. Finally, the participation ratio for that neighborhood was taken as:

PR increases as variance is distributed more evenly across orthogonal components (higher local dimensionality), approaching 1 when variability is dominated by a single component (neuron) and approaching K (number of neighbors) when contributions are comparable across all components. The procedure was repeated with each neuron serving as the anchor, yielding one PR value per neuron.

Smoothing feature distributions

In Figures 4A, 5, 8S1 and 9G, distributions of features were smoothed with Kernel Density Estimation (KDE). Specifically, for a distribution, the density at point x was computed as

Where {x_1, x_2, x_3,…, x_n} represent the distribution of interest, K is the Gaussian kernel, and the bandwidth h was set by Silverman’s rule (Läuter 1988):

with n being the number of observations and s being the sample standard deviation.

PCA-based alternative to feature-based analyses

As an alternative to our feature extraction method for summarizing population structure in trial-averaged activity, we applied PCA to the vectorized PETHs of all modulated neurons and used the first 20 scores as a feature space. For increased interpretability, we performed a VARIMAX rotation on the scores to sparsify them without distorting the space. The rotated scores’ spatial structure was visualized in Figure 5S1 and their likelihood maps in Figure 9S1.

Gradient calculation

In Figure 6, we asked whether the single neuron PETH feature values changed more over the same spatial distance in some locations than others. This corresponds to a spatial derivative of the metric field. To compute this for each map in Figure 6, we first smoothed the spatial distribution of the feature map using the irregular Gaussian smooth function as described in the section “Computing map images from single-cell properties”. We then gridded the smoothed map values into square bins as described above, and filled missing interior values with linear interpolations of their neighbors. At each pixel, we then computed the 2D gradient magnitude using MATLAB’s centered finite-difference method. Specifically, if M(x, y) denotes the gridded feature map, then the gradient magnitude |∇M(x, y)| was calculated as:

This produced a scalar map of local changes in the corresponding feature map. For the version in Figure 9S4, the two partial derivatives (x and y) were computed for each of the four subpopulation likelihood maps. The quadratic mean was then calculated as above over the set of individual feature (or likelihood) gradient magnitude maps.

Random boundary analysis

In Figure 6G,H we generated random boundaries through rejection sampling. Potential midpoint coordinates were randomly selected uniformly within the bounding box of the map, and an orientation was selected from a uniform distribution of 0-360°. We then computed the endpoints for a 1 mm line segment with that center and orientation. This boundary candidate was accepted if two criteria were met. First, both endpoints must have been within a non-empty pixel. Second, at least 200 neurons must have been in the valid zones on each side of the boundary (between 50 µm and 500 µm from the boundary). To account for non-uniform distributions of neurons, we downsampled neurons on the denser side of each random boundary to match the count on the sparser side.

To assess how well we could classify neurons as above or below a random boundary using the neurons’ feature vectors, we trained a binary SVM with a radial basis function kernel and parameter optimization via Bayesian hyperparameter search. Decoding performance was evaluated using 5-fold cross-validation.

Nonlinear dimensionality reduction

We used t-SNE to visualize the structure of our neural population’s feature space. Each neuron was represented by a five-dimensional feature vector comprising the PETH feature values described above. To improve the representation of large-scale structure, we seeded the t-SNE optimization with the top two principal components of the feature matrix (Kobak and Linderman 2021). The affinity metric was chosen to be the Pearson correlation. t-SNE perplexity hyperparameter was chosen according to the N/100 rule of thumb (Kobak and Berens 2019). The t-SNE space produced a two-dimensional projection of our neural data with coordinates of Y_i = (y_i,1, y_i,2) for each neuron i.

To produce contours in the t-SNE space for each area, we first selected the points for a given area and performed two-dimensional Gaussian KDE on these points to produce an estimate of probability density. We then drew contours (contour in Matlab) at the 90th, 95th, and 98th percentiles of the peak. To verify that our results were not specific to t-SNE, we also applied two-dimensional Uniform Manifold Approximation and Projection (UMAP; McInnes et al. (2020)) and obtained qualitatively similar results.

Gaussian Mixture Modeling and cluster selection

For the analyses in Figures 8, 9 and 10, we modeled the PETH feature vectors from each area using Gaussian Mixture Models (GMMs). For each candidate number of clusters k, we fit a GMM with full covariance matrices using the Expectation-Maximization algorithm (fitgmdist in Matlab).To select the optimal number of mixture components, we used the Integrated Completed Likelihood (ICL) criterion (Bertoletti et al. 2015), computed as:

where P_ij denotes the posterior probability that neuron i belongs to cluster j, BIC represents the Bayesian Information Criterion (Schwarz 1978), and ε was chosen as 10^-4. For each area, we fitted GMMs with one to ten components, and selected the fit that produced the maximal ICL.

Likelihood computations

Following model selection for each region, we computed the log likelihood of every neuron according to the GMM of each anatomical region. The log-likelihood ℒ_i of neuron i given the full model for a single area was defined as:

where K is the number of components in the GMM and π_j is the mixing weight for component j. In addition to the full-mixture model, log-likelihoods were calculated for each component within an anatomical area’s full mixture model to quantify the spatial footprint of each Gaussian component and to compare this spatial map across components of different anatomical areas. These log likelihoods provide a computational approach for understanding the structure of single neuron feature properties across the sensorimotor cortex.

GMM component distances and clustering

In Figure 9B, we computed the distance between each pair of GMM components found in different areas. To do so, we computed pairwise Bhattacharyya distances between the Gaussians. This distance captures differences in the means (as a Mahalanobis-like distance), and differences in the variances and covariances (differences in the orientations, shapes, and sizes). Each Gaussian k was represented by its mean vector μ_k ∊ ℝ^N and covariance matrix Σ_k ∊ ℝ^N×N, where N is the feature dimensionality (5). The Bhattacharyya distance D_B between multivariate Gaussians i and j was defined as

For display, the distance matrix produced by the Bhattacharyya distance was ordered using hierarchical clustering with Ward’s method (Fig. 9B). The number of clusters (‘subpopulations’) was determined by simple examination of this matrix. Each cluster comprised GMM components with heavily overlapping distributions in the feature space. Clusters were further validated by comparing the spatial distributions (in cortex) of the member neurons (Fig. 10E), which were very similar for the components grouped in each cluster.

As the simplest way to fairly represent the components that were grouped to form subpopulations, for further analysis we modeled each subpopulation as arising from a GMM comprising the Gaussian components that were grouped together by clustering, weighted equally.

Acknowledgements

The authors thank P. Ravishankar for animal care assistance and L. Pinto for helpful comments on the manuscript. This work was funded by The University of Chicago, NIH-NINDS R01 NS121535 (MK), the Simons Foundation (MK), the Whitehall Foundation (MK), the NSF-Simons National Institute for Theory and Mathematics in Biology via grants NSF DMS-2235451 and Simons Foundation MP-TMPS-00005320 (MK), and NIH T32 NS121763 (HG).

Additional information

Data availability

Data will be made available on FigShare upon publication.

Code availability

Code will be made available on GitHub upon publication.