Striatal ensemble activity in an innate naturalistic behavior

Self-grooming is an innate, naturalistic behavior found in a wide variety of organisms. The control of rodent grooming has been shown to be mediated by the dorsolateral striatum through lesion studies and in-vivo extracellular recordings. Yet, it is unclear how populations of neurons in the striatum encode grooming. We recorded single-unit extracellular activity from populations of neurons in freely moving mice and developed a semi-automated approach to detect self-grooming events from 117 hours of simultaneous multi-camera video recordings of mouse behavior. We first characterized the grooming transition-aligned response profiles of striatal projection neuron and fast spiking interneuron single units. We identified striatal ensembles whose units were more strongly correlated during grooming than during the entire session. These ensembles display varied grooming responses, including transient changes around grooming transitions or sustained changes in activity throughout the duration of grooming. Neural trajectories computed from the identified ensembles retain the grooming related dynamics present in trajectories computed from all units in the session. These results elaborate striatal function in rodent self-grooming and demonstrate that striatal grooming-related activity is organized within functional ensembles, improving our understanding of how the striatum guides action selection in a naturalistic behavior.

The selection and initiation of actions is controlled by the striatum [23][24][25][26][27][28][29] . The striatum has been shown to encode action space 30 and to flexibly combine behavioral motifs into actions 25 . Cells within the striatum, predominantly consisting of striatal spiny projection neurons [31][32][33] , are organized into functional clusters of co-active units; this clustering is considered to be important for striatal network dynamics 34,35 and behavioral control 30,[36][37][38][39][40] . The striatum has been implicated in the production of self-grooming behavior, and lesions of the striatum disrupt grooming bouts 41 . Within the striatum, neurons of the dorsolateral striatum 42,43 and nearby central striatal regions 44 display grooming-related activity. Yet, how populations or ensembles of neurons in the striatum-defined here as sets of neurons that are more likely than chance to be co-active-encode self-grooming remains unclear. Given the likely importance of ensemble activity in the striatum, elaborating whether grooming-associated neural activity maps onto striatal ensembles is significant, and it remains to be elucidated.
Here we recorded simultaneous activity of populations of neurons in the dorsolateral striatum of freely moving mice during spontaneous grooming using extracellular probes. Because grooming can be a relatively rare behavior, many hours of data must be acquired to capture an adequate sample of its neural correlates. To overcome this obstacle, we developed a semi-automated approach to detect grooming bouts from behavioral videos using 3D pose estimation and postural heuristics. We found striatal projection neurons (SPN) and fast spiking interneurons (FSI) with temporally diverse grooming-related activity. Furthermore, we identified striatal ensembles that encode core parameters of grooming bouts, including the transitions in and out of individual bouts, as well as bout duration.  . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 23, 2023. ; https://doi.org/10.1101/2023.02.23.529669 doi: bioRxiv preprint To record striatal activity during mouse spontaneous grooming, we implanted 64-channel electrodes into the dorsolateral striatum of adult male and female mice (Figure 2-figure supplement 1A-D, N=6 mice). We classified well-isolated units as putative striatal spiny projection neurons (SPN) and putative fast spiking interneurons (FSI) on the basis of their firing rates and spike waveform, following established criteria 45-47. Units that did not fit the criteria for SPNs or FSIs were labelled as 'other' and excluded from further analyses. (Figure 2A, SPN: 88.2% 2,755 units, FSI: 3.2% 100 units, other: 8.6% 269 units).
To identify example individual units with grooming-related activity, we isolated units whose activity within a two-second window around the start or end of grooming was two standard deviations above mean activity during a grooming-free baseline period. All subsequent analyses were performed on the full dataset. Example SPNs and FSIs with grooming related activity are shown in Figure 2B-E. We observed units displaying diverse grooming-related responses, including increased activity at the start of grooming, end of grooming, and both start and end of grooming, as well as units that showed elevated activity for the duration of grooming, and units that showed reduced activity for the duration of grooming. Left: units were categorized as either SPNs, FSIs, or 'other' by their firing rates, spike waveform peak width, and duration between spike waveform peak to valley (SPN: 88.2% 2,755 units, FSI: 3.2% 100 units, other: 8.6% 269 units). Right: average +/-SEM spike waveforms for units in each category. B. Activity of 4 example SPNs during grooming. Neural activity are aligned to the start of grooming bouts, denoted by t = 0. Top: spike raster plots for the example neuron during each grooming bout in the given session. Grooming bouts are sorted by grooming bout duration, denoted by the green rectangles. Bottom: unit average firing rate aligned to grooming start and normalized to [0, 1]. C. Activity of 4 example SPNs as in B, but for aligned to groom end. D. Activity of 4 example units as in B, but for FSIs aligned to groom start. E. Activity of 4 example units as in B, but for FSIs aligned to groom end.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  To characterize different grooming-associated responses in the recorded striatal units, we first performed principal components analysis (PCA) 48 on the grooming transition-triggered-average activity of all recorded SPNs and FSIs. We performed PCA separately on SPN activity aligned to groom start, SPN activity aligned to groom end, FSI activity aligned to groom start, and FSI activity aligned to groom end. The first principal component (PC) from each group reflects units that undergo a step-like increase or decrease in their activity at the grooming transition ( Figure 3A, B, F, G, K, L, P, Q, SPNs aligned to groom start: 336 units with positive contributions to the first PC (denoted "> 0"), 640 < 0, SPNs aligned to groom end: 526 > 0, 424 < 0, FSIs aligned to groom start: 66 > 0, 1 < 0, FSIs aligned to groom end: 22 > 0, 44 < 0). The second PC from each group reflects units whose activity transiently peaks or decreases around the grooming transition ( Figure 3C, D, H, I, M, N, R, S, SPNs aligned to groom start: 364 > 0, 292 < 0, SPNs aligned to groom end: 164 > 0, 518 < 0, FSIs aligned to groom start: 4 > 0, 19 < 0, FSIs aligned to groom end: 5 > 0, 19 < 0). Taken together, these results show that single unit SPNs and FSIs encode transitions into and out of an innate naturalistic behavior.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The event-triggered average responses of striatal neurons reveal several patterns of firing rate modulation aligned with the start or end of grooming. However, because this analysis relies on trial-averaging, it cannot be used to determine whether the striatum contains ensembles of neurons that are consistently coactive. To address this question, we looked for functional clusters, or 'ensembles' of co-active striatal neurons with grooming related activity. To achieve this, we first constructed a matrix containing the grooming activity of all units in a given session, where the activity of each unit plus 5 seconds before and after grooming was concatenated. This matrix is shown for an example session in Figure 4A; dendrogrambased sorting of recorded units makes evident the presence of synchronous activity in groups of units. Synchrony among units is also suggested by the block-diagonal structure of the correlation matrix for this example session ( Figure 4B). This synchrony is further suggested by comparison to a time-shuffled version of the data ( To test for the presence of ensembles and to estimate the number of potential ensembles within each session, we used an eigenvalue-based statistical method 49,50 , where we identify dimensions of the neural covariance matrix that capture more variance than expected by chance and found each session to have three to ten ensembles (Figure 4C, 6.2 ± 0.3 ensembles, 33 sessions, see Methods for details). Next, using the meta-k-means clustering algorithm 51 , we identified the ensembles present in each session. Meta-kmeans is an extension of k-means where one performs N iterations of k-means, keeping track of how many times each pair of units gets clustered together. The proportion of times that each pair gets clustered together is used to assign intermediate clusters before a final merging step. We chose to use meta-k-means for two reasons. First, meta-k-means does not always result in cluster assignment for all units, which is preferable because our recording configuration only samples a small subset of neurons in striatum and we do not expect all recorded units to belong to an identified striatal ensemble. Second, meta-k-means allows for the final number of identified ensembles to be greater or less than the initial choice of k, which is preferable because we do not have a priori knowledge of ensemble numbers present in each session.
We found that 80% of recorded SPNs (1,074 units) and 70% of recorded FSIs (41 units) were clustered into ensembles of two or more units (Figure 4-figure supplement 1C). The majority of clusters were . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 23, 2023. ; https://doi.org/10.1101/2023.02.23.529669 doi: bioRxiv preprint composed of only SPNs, and FSIs were almost always assigned clusters together with SPNs (Figure 4figure supplement 1D; 88% (214 clusters) SPN only, 1% (1 cluster) FSIs only, 11% (28 clusters) SPN and FSIs). Within our identified ensembles, units in the same ensemble tended to be spatially closer than units that are not within the same ensemble ( Figure 4D), consistent with striatal ensembles detected using imaging techniques 30,38,39 . During grooming, the pairwise correlation between units within the same cluster is higher than for the correlation between units not in the same cluster, whether those units are unclustered or belong to a different cluster ( Figure 4E). The distributions of pairwise correlations computed using data from the whole session are shown in Figure 4F. This in-versus-out of cluster difference is weaker when comparing unit activity during the whole session, rather than during grooming ( Figure 4G, grooming AUC: 0.81, whole session AUC: 0.69). (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 23, 2023. ; https://doi.org/10.1101/2023.02.23.529669 doi: bioRxiv preprint (blue). Plot depicts the mean ± SEM of the cumulative distribution across mice (4 mice). E. Distribution of pairwise unit correlations during grooming for each pair of units that are within the same cluster (pink 2,747 pairs) and each pair of units that are not within the same cluster (blue 27,421 pairs). Left: histogram. Right: cumulative distribution. F. Same as in E, but for pairwise correlations computed from activity during the entire session including grooming times. G. Comparison between the difference in pairwise unit correlations for units within and not within the same cluster computed during grooming (green) and during the whole session (orange) (grooming AUC: 0.81, whole session AUC: 0.69).
The distribution of the number of ensembles identified in each session is shown in  Arrows are in the same position as in Figure 4A highlighting the absence of synchrony after shuffling. Each unit's activity is normalized to the range from zero to one. B. Correlation matrix for the shuffled activity shown in A. C. Total number of SPNs and FSIs that were clustered (pink) and unclustered (blue) (4 mice, 33 sessions, 243 clusters). D. Total number of clusters that comprised of only SPNs (purple), only FSIs (green), and of both SPNs and FSIs . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  Having demonstrated the existence of ensemble activity during grooming, we next asked whether ensemble activity was enriched during specific timepoints during grooming. A heatmap displaying the average grooming activity for all significant ensembles is shown in Figure 5A (112 ensembles, 4 mice, 33 sessions). To visualize the average grooming activity of ensembles across grooming bouts of varying duration, we first linearly time-warped ensemble activity during grooming. We found that ensemble peak activity is enriched around the transitions into and out of grooming ( Figure 5B), but that the dynamics of grooming ensemble activation are diverse. For example, we found ensembles with higher activity during grooming, lower activity during grooming, and with peak activity at the start or end of grooming, with individual cluster examples of these patterns shown in . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made To contrast the encoding of grooming-related activity in ensembles vs unclustered units, we used nonnegative matrix factorization to perform dimensionality reduction on the neural activity during grooming within a single session. The trajectories traced out by the population activity of all units during all grooming bouts in an example session are shown in Figure 5G (65 units). During most grooming bouts the population activity moves out along factor 1, then up along factor 2, and finally returns toward the origin along factor 3, corresponding to population encoding of grooming onset, maintenance, and termination. The neural trajectories computed from unclustered units do not retain the dynamics present in the trajectories from all units or the clustered ones (Figure 5-figure supplement 1D, 16 unclustered units; Figure 5H, 49 clustered units). Notably, the ensemble trajectories computed from the activity of ensembles retain most of the dynamics seen in the trajectories computed from the activity of single units (Figure 5I, 13 ensembles). Heatmap depicting the percent of units within each ensemble that are active during grooming averaged across all grooming bouts in a given session. Prior to averaging across grooming bouts, ensemble activity during each grooming bout was interpolated to a fixed duration, excluding the 5 seconds before and after grooming. Ensembles are in the same order as in Figure 5A and data are normalized to the range from zero to one (112 ensembles, 4 mice, 33 sessions). B. Distribution of times when the percent of active units in an ensemble peaked (data as in A). C. Representative example of an ensemble with decreased average activity throughout the duration of grooming (grey region as in Figure 5C). D. Neural trajectories traced out in factor space by the population of unclustered units during all grooming bouts in an example session (16 units, 22 grooming bouts). Visualization elements as in Figure 5G.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

Discussion
Here, we developed a semi-automated approach to identify grooming events in mouse behavioral videos. We recorded population activity in the dorsolateral striatum of freely moving mice and found SPNs and FSIs with activity that encodes the start and end of spontaneous grooming bouts. Previous studies reported changes in striatal single unit activity around the start and end of naturalistic 42,43,52 and learned behaviors 24,53,54 . Our single unit data elaborate the changes in striatal unit activity around the transitions of a naturalistic behavior. Single unit responses were heterogenous and included units with changes in activity at the start or end of grooming, as well as units that were active or silenced throughout the duration of grooming. We identified striatal ensembles with units that were more correlated during grooming than during the entire session. These striatal ensembles encoded grooming start time, end time, and bout duration. Single session trajectories computed from ensemble activity retain most of the dynamics present in trajectories computed from all units in the session.
The distribution of pairwise distances between units within an ensemble is shifted towards smaller distances, compared to the distribution of distances for units that do not belong to the same ensemble. This is consistent with striatal ensembles detected using imaging techniques 30,38,39 although our identified ensembles are less spatially compact than those previously reported. Notably, however, striatal 2-photon imaging approaches and microendoscopic imaging of genetically encoded calcium indicators yield a recording area that extends along the medial-lateral and anterior posterior axes (i.e., horizontally), whereas our four-shank electrophysiological recordings yield a recording field that extends along the anteriorposterior and dorsal-ventral axes (i.e., vertically), potentially accounting for this difference in findings. Together, the ensemble distances recorded via both methods suggest the possibility that striatal ensembles are isotropically organized, which is well-aligned with studies of cortical and thalamic projection patterns to the striatum [55][56][57] .
While the ensembles we characterized were identified during grooming behavior, it is possible that they are not grooming-specific. One possibility for the functional organization of striatal ensembles is that ensemble membership is stable across different behaviors. Then, ensemble activation encodes some common movement motif or state that is present across multiple behavioral settings. However, the units within our identified striatal ensembles were more correlated during grooming than during the entire recording session. Our data therefore support an alternative to the stable membership hypothesis, wherein individual units can be members of multiple ensembles, such that when a given unit is active, the behavior the animal is performing can only be determined by looking at the population level. For example, a set of units might be members of a single ensemble that encodes the start of grooming, but each unit could also be part of different ensembles that encode aspects of other natural behaviors, such as eating or walking. This 'mixed selectivity' is broadly consistent with previously reported data on striatal activity patterns and cluster memberships that are not conserved across divergent behavioral motifs 30,39 .
Striatal FSIs are parvalbumin-containing GABAergic interneurons that make up approximately 1% of striatal neurons 33,58,59 . FSIs receive direct cortical input, form inhibitory synapses onto SPNs 33,60,61 , and are interconnected amongst themselves via gap junctions on their dendrites 62 . Previous studies found . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 23, 2023. ; https://doi.org/10.1101/2023.02.23.529669 doi: bioRxiv preprint uncoordinated, idiosyncratic task-related changes in FSI activity 63 , as well as changes in FSI activity that correlated with movement features like velocity 64 , with one study showing grooming-related activity in a small number of recorded putative FSIs 43. We recorded units that match spike waveform features and firing rates of FSIs 45,59 and found them to have grooming related activity. These units displayed decreased or increased activity during grooming, as well as transient activity changes at the transitions into and out of grooming. Consistent with strong connectivity present between FSIs and SPNs 65 , a majority of FSIs (70%) were part of identified ensembles and almost all FSIs were part of ensembles together with SPNs.
Striatal control of movement is mediated via two parallel pathways: the direct and indirect pathway, composed of dSPNs and iSPNs [66][67][68] . Previous recordings of direct and indirect pathway striatal ensemble activity have found that d-and iSPNs exhibit a similar time course of behavior-related changes in their activity 23,30,38 and have similar ensemble organization 30 , suggesting mixed pathway membership within ensembles. Indeed, simultaneous recording of direct and indirect pathway SPNs showed that striatal ensembles are composed of d-and iSPNs in equal proportions 39 . These prior data suggest that striatal ensembles encoding self-grooming might be composed of both d-and iSPNs, in similar proportions and with similar encoding of grooming, although this remains to be confirmed.
Within bouts of rodent self-grooming, animals perform highly stereotyped grooming sequences called syntactic grooming 69,70 . This sequence is comprised of four distinctive phases and found in a variety of species [8][9][10] . A large body of work has established a function for the striatum in the sequential ordering of syntactic grooming phases 41,52,71,72 . On a longer time scale, grooming is one of many behaviors an animal can perform at a given time. Since self-grooming, although broadly important, is not often an immediate concern, grooming is considered a low priority behavior 73 . Thus, the control of grooming has been characterized by a disinhibition model, in which grooming takes place in the time left over by other higher priority behaviors when grooming would be inhibited 74 . Previous studies have demonstrated a striatal role in behavioral sequencing on this longer time scale as well 25,30 . Thus, self-grooming represents an ethologically meaningful behavioral paradigm providing a path to study neural control of an innate, conserved behavior at multiple spatiotemporal scales.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

Subjects
Animals were handled according to protocols approved by the Northwestern University Animal Care and Use Committee (protocol number: IS00009022). Adult male and female C57BL/6J mice (p57-105) were used in this study (The Jackson Laboratory, RRID:IMSR_JAX:000664). All mice were group-housed in a humidity-controlled, ambient temperature facility, with standard feeding, 12 hr. light-dark cycle, and enrichment procedures.

Electrode implantation
Mice were anesthetized with isoflurane (3% for induction, 1.5-2% for maintenance) and placed on a small animal stereotax frame (David Kopf Instruments, Tujunga, CA). A 4-shank 64-channel silicon electrode (part #A4x16-Poly2-5mm-23s-200-177-H64LP_30mm mounted on a dDrive-m, NeuroNexus Technologies, Ann Arbor, MI) was implanted into the dorsolateral striatum (0.5 mm AP, 2 mm ML, and lowered 0.25 mm) and secured with Vetbond (3M, Maplewood, MN) followed by dental cement (Micron Superior 2, Prevest DenPro, Jammu, India, or C&B Metabond, Parkell, Edgewood, NY). The electrode shanks were aligned with the brain's anterior-posterior axis. A skull screw was connected to the electrode ground wire and fastened to the skull above the ipsilateral cerebellum. The screw was secured with Vetbond followed by dental cement. Mice were monitored following the surgery to ensure a full recovery and were administered post-operative analgesics. Mice recovered for at least 5 days after implantation.

Behavioral recording
On the first day of experimentation, mice were acclimated to the experimental arena for 30 minutes. All subsequent experimental sessions were 120 minutes long. The experimental arena was an equilateral triangular arena constructed from transparent acrylic (Figure 1A. 12-inch sides and height, 1/8-inch thick. Part #8536K131, McMaster-Carr, Elmhurst, IL) fastened together with clear epoxy (Part #31345, Devcon, Solon, OH). While in the arena, mice were given DietGel (ClearH2O, Westbrook, ME) for food and hydration. One mouse was in the arena at a time. The arena was cleaned with 70% ethanol after each session. The arena was enclosed by a dark curtain and illuminated by infrared LEDs. Experiments were conducted during the animal's active phase.
Mouse behavior was captured with three side-view cameras at 125 fps using a custom fork of the campy Python package 75 (Figure 1A. image size 1440 x 608 pixels. cameras: Part #BFS-U3-16S2M-CS Teledyne FLIR, Wilsonville, OR. Lenses: Part #8595755548 Yohii, China). To synchronize recordings from all cameras, each frame grab was triggered by a TTL pulse sent from a microcontroller (Arduino, Turin, Italy). To facilitate alignment of videos and neural data, each frame grab TTL pulse was also sent to the electrophysiology acquisition board and the video and neural recordings were simultaneously initiated by pressing a button connected to the microcontroller.
Electrophysiological recording . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 23, 2023. ; https://doi.org/10.1101/2023.02.23.529669 doi: bioRxiv preprint A 64-channel headstage (part #C3325 Intan Technologies, Los Angeles, CA) was connected to the implanted electrode. In the arena, the headstage was connected to a 12-channel commutator supported by a balance arm (Part # FL-12-C-MICRO-BAL Dragonfly Inc., Ridgeley, WV). Neural recordings were acquired at 30 kHz with an Open Ephys acquisition board 76 and the Open Ephys GUI 76 (GUI version 0.5.5, Open Ephys, Atlanta, GA). To facilitate alignment of videos and neural data, the video and neural recordings were simultaneously initiated by pressing a button connected to the microcontroller.
Neural recordings were spike sorted offline using KiloSort3 77 in MATLAB (MathWorks, Natick, MA). We only considered units labeled as 'good' by KiloSort. Additionally, we manually inspected the waveforms of all units labeled as 'good' and excluded units with waveform shapes that did not resemble an action potential. Recordings were performed from each animal until the number of 'good units, as identified by KiloSort, stayed below 15 for 2 consecutive days (27.3 ± 4.3 days since implant).

Unit classification
All recorded units were classified as either putative striatal projection neurons (SPN), fast-spiking interneurons (FSI), or 'other' based on their spike waveform and mean firing rate across the entire session following established criteria [45][46][47] . Units were classified as putative SPNs (88.2 %) if they had peak width > 150 µs, peak-valley interval > 500 µs, and mean firing rate ≤ 10 Hz. Units were classified as FSIs (3.2 %) if they had peak width ≤ 150 µs, peak-valley interval ≤ 500 µs, and mean firing rate ≥ 0.1 Hz. Units that were not identified as either putative SPNs or FSIs were labeled 'other' (8.6 %) and excluded from all analyses.

Histology
After the recordings, coronal brain sections were obtained from all mice to determine electrode placement. Mice were deeply anesthetized with isoflurane and transcardially perfused with 4% paraformaldehyde (PFA) in 0.1 M phosphate buffered saline (PBS). Brains were post-fixed for 1-2 days and washed in PBS. Brains were then sectioned on a vibratome (Leica Biosystems, Wetzlar, Germany) at 60 µm from frontal cortex to posterior striatum and then dried and cover slipped under glycerol:TBS (9:1) with Hoechst 33342 (2.5 µg/ml, Thermo Fisher Scientific, Waltham, MA). Sections were imaged on an Olympus VS120 slide scanning microscope (Olympus Scientific Solutions Americas, Waltham, MA) with DAPI for localization of cell nuclei and FITC for background fluorescence to enable electrode localization.

Electrode localization
Coronal slices containing electrode tracts were processed using WholeBrain software package 78 in R 79 . Histological slices were analyzed and registered to the Allen Mouse Brain Common Coordinate Framework (CCFv3) 80 as described previously 81 . Striatal sections with electrode tracts were analyzed at 60 µm intervals and frontal lobe cortical sections, which rarely contained electrode tracts, were analyzed at 100-200 µm intervals. A total of 6-18 section images were analyzed per mouse brain. Electrode placement for each brain section was denoted at the most ventral location that the electrode tract was observed. Medial-lateral and dorsal-ventral coordinates registered to the CCFv3 were obtained for each section and plotted on a representative coronal section using the WholeBrain package.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 23, 2023. ; https://doi.org/10.1101/2023.02.23.529669 doi: bioRxiv preprint

Grooming identification
We identified mouse grooming bouts in behavioral videos using a semi-automated 4-step procedure ( Figure 1C). As detailed further below, first, we tracked the mouse limb positions in 2D; second, we triangulated to the 2D limb positions to 3D; third, the 3D limb positions were used to isolate likely grooming times; and fourth, using the likely times of grooming we manually identified grooming bout start and stop times.
For 2D pose estimation, we used DeepLabCut (version 2.2.0.5) 82,83 to detect 15 keypoint positions (snout, both paws, wrists, elbows, shoulders, eyes, ears, and hind paws) on the frames from each camera. We labeled the 15 keypoints in 1,627 frames from 14 videos, 3 animals, and all 3 cameras. We trained a ResNet-50-based neural network 84,85 with default parameters on 90% of the labeled frames for 300,000 iterations. We validated with 1 shuffle and got a train error of 3.88 pixels and test error of 6.64 pixels (image size: 1440 x 608 pixels). We used a p-cutoff of 0.6 to condition the x and y coordinates for future analyses. With the p-cutoff, the train error was 3.27 pixels and test error 5.06 pixels. This network was then used to analyze all videos from the same experimental setting.
For 3D triangulation of the 2D poses we relied on Anipose 86 . We calibrated our 3 cameras with a 3-minute video of a ChArUco board throughout the camera fields of view (125 FPS. Image size 1440 x 608 pixels). The ChArUco board was printed on paper and taped onto a stiff plastic board. The ChArUco pattern had 7x10 squares containing 4-bit markers and a dictionary size of 50 markers. Before triangulation, we applied a Viterbi filter to the 2D poses with F = 12 frames. The triangulation was performed via an optimization constrained on the distance between the two eyes and the distance between each eye and the snout. We chose scale_smooth = 5, scale_length = 2, and score_threshold = 0.3.
To isolate times when the mouse was likely grooming, we identified frames where the mouse was rearing, and its paws were near its head. To identify frames where the mouse was rearing, we set thresholds on the snout height and the distance between the midpoint of the eyes and the midpoint of the hindlimbs. To identify frames where the mouse's paws were near its head, we set thresholds on the distance between the midpoint of the paws and the snout, the distance between the midpoint of the paws and the midpoint of the eyes, and height of the midpoint of the paws. Per-frame labels within 1.2 seconds were merged to form predictions over a window. Predicted grooming bouts less than 2 seconds long were discarded. Parameter values were evaluated by generating predictions for a set of videos and viewing the predictions aligned to the video in BENTO 87 .
Grooming bout start and stop times were manually refined by 4 trained annotators, with each annotator labeling a different subset of videos. Annotators used the predictions from the previous thresholding step to navigate the videos and score grooming behaviors. The videos were annotated in VLC media player 88 using the Time v3.2 89 and Speed controller 90 extensions. To assess inter-annotator reliability, all annotators were also given the same 2-hour session to annotate. The average pairwise Jaccard index (or Intersection over Union) computed on the annotations was 0.76 ± 0.04.
Identification of units with grooming-related activity . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made To isolate units with grooming related activity, we first down sampled activity to 2 Hz, by summing spike counts within non-overlapping 500 ms bins. We defined units as having grooming-related activity if their average activity at the start or end of grooming was 2 standard deviations above their average activity during a 3 second-long, grooming-free baseline period. For activity aligned to groom start, we used the activity from 5 to 2 seconds before grooming start as the baseline and compared that to the activity from 1 second before until 1 second after the grooming start time. For activity aligned to groom end, the activity from 2 to 5 seconds after grooming was chosen as the baseline and was compared to the average activity from 1 second before until 1 second after groom end. These selection criteria were only used to identify neurons to serve as examples in Figure 2; all subsequent analyses were performed on the full dataset.

Characterization of trial-averaged grooming responses
We characterized the grooming responses exhibited by the recorded units by computing each unit's eventtriggered average response in a ± 5-second window relative to the start and (separately) end of grooming. Responses were down sampled to 4 Hz, by summing spike counts within non-overlapping 250 ms bins. To simplify interpretation, we restricted averaging to grooming bouts that were 'well-isolated.' Specifically, because we included the 5 seconds before and after each bout in this analysis, grooming bouts that started within 10 seconds of the end of the previous bout were not included in this analysis. An exception to this is instances with less than 3 seconds between the end of one grooming bout and the start of the next, in which case we merged these two annotations into a single bout.
To better visualize and interpret grooming-associated activity we next performed principal component analysis (PCA) 48 on the set of event-triggered averages. Specifically, we performed PCA separately on SPN activity aligned to groom start, SPN activity aligned to groom end, FSI activity aligned to groom start, and FSI activity aligned to groom end. For each cell type and condition combination, we concatenated the associated units to form a units-by-time matrix, relative to groom start/stop. The response of each unit in this matrix was then Z-scored, after which PCA was used to extract the first two principal components.
To group units into 'response types' for visualization, we examined the weight (or 'score') for each unit's contribution to the largest and second-largest principal component. Units for which the absolute value of this weight was greater for the first principal component than the second formed the "PC1" response group, while those with greater-magnitude weights for the second principal component formed the "PC2" response group. Finally, we sorted the units within each group according to their contribution to the first (PC1) or second (PC2) principal component of the dataset (weight * activity).

Grooming-associated striatal ensemble identification
To identify grooming-associated striatal ensembles within a given session, we took the neurons-by-time matrix of recorded spiking from the full two-hour recording session and excluded all frames that did not occur during grooming or within a 5-second window before or after grooming. Due to the low firing rate of SPNs, it was rare to find units that were co-active at high sampling rates. Therefore, to instead identify units that were frequently active at around the same time, we binned spike counts at a sampling rate of . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 23, 2023. ; https://doi.org/10.1101/2023.02.23.529669 doi: bioRxiv preprint 0.667 Hz (bin size of 1.5 seconds). The spike counts for each unit were then normalized to between zero and one, and further smoothed by convolving with a Gaussian filter with a standard deviation of 3 seconds.
Using this smoothed and coarsened estimate of cell activity, we next applied meta-k-means 51 to identify sets of cells that were often co-active. We chose to use meta-k-means for clustering, because it does not force all units in a session to be assigned to a cluster, and it allows for the number of ensembles to be greater or less than the initial choice of k. Briefly, meta-k-means employs repeated runs of standard kmeans clustering to identify groups of units that are consistently clustered together. As with standard kmeans, meta-k-means requires the user specify an initial number of clusters k, which we set to be the square root of the number of units in the given session, as in Barbera et al 38 . We then ran 1,000 repeats of the k-means clustering algorithm on the matrix of filtered grooming-related activity described above, with units as variables and time bins as observations. Cluster centroids for each repeat were initialized using the greedy k-means++ algorithm 91 . As in prior work 30,38,92 , units that were assigned to the same cluster in >80% of k-means runs were considered to be part of the same meta-k-means cluster.
When merging clusters as part of the meta-k-means algorithm, we used the silhouette score 93 to evaluate the outcome of a potential cluster merge. To ensure adequate sample size for downstream analyses, ensemble identification was restricted to sessions with at least 30 SPNs and FSIs.

Statistical analysis of striatal ensembles
To compute a statistical estimate on the number of ensembles present in each session we used an eigenvalue based statistical method from Peyrache at al 49,50 . Briefly, if two neurons have correlated firing, then we expect to observe common fluctuations in their spiking activity during recording. This can be visualized by plotting the firing rates of the two neurons against each other, in the form of a direction in neural activity space along which observations tend to be distributed. Conversely, without any correlated firing, the variance in the dataset will be roughly equally distributed in all directions, and one would not expect to see structure when plotting the activity of two cells against each other. These relationships hold for any n-dimensional set of neurons, where structure emerges from correlated activity. The eigenvectors of the covariance matrix of a dataset represent the directions of maximally shared variance, thereby capturing the correlations present in the data. Each eigenvalue quantifies the variance of the data along the axis defined by its corresponding eigenvector. Thus, in a matrix of neurons-by-time, the eigenvalues of its covariance matrix capture the correlations among neurons (i.e., ensembles) and comparing the eigenvalues to those computed from the covariance matrix of a shuffled version of the data provides a measure of the number of ensembles present.
To apply this method, we first binned spike counts at the same sampling rate as for clustering (0.667 Hz, bin size of 1.5 seconds) and Z-scored the activity of each neuron. Next, we shuffled the time bins for each neuron independently, and computed the maximal eigenvalue of the covariance matrix of the shuffled data. We repeated this 5,000 times to form a null distribution and counted the number of eigenvalues from the covariance matrix of the actual data that were above the 99 th percentile of this null distribution.
Time-normalized grooming-related activity . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made To compare ensemble average activity during grooming across grooming bouts of varying duration, we linearly temporally rescaled the activity of each unit during each grooming bout to a single fixed length. We chose the length for the interpolated values from the distribution of grooming bout durations from all clustered sessions. To avoid 'blurring' activity associated with the start and end of grooming bouts, the activity during the 5 seconds before and after each bout was not interpolated.
Bootstrap significance testing of grooming-related ensemble activity To test for periods during grooming when the grooming ensemble average activity was significantly higher or lower than chance, we computed the mean, 2.5 th , and 97.5 th percentiles of the distribution of randomly sampled ensemble activity. To generate this distribution, we computed the population-average ensemble activity during 1,000 random duration windows of time sampled throughout the 2-hour recording session. Start times were sampled uniformly from 5 seconds into the recording session until 5 seconds before the end of the session. Window durations were sampled from the distribution of observed grooming bout durations across all sessions.

Neural trajectories
To visualize neural trajectories, neural activity during grooming was decomposed onto 3 components using non-negative matrix factorization 94 . For the neural trajectories, we decomposed a matrix containing each unit's activity during all grooming bouts within that session. For each unit, we concatenated its activity during all grooming bouts including the five seconds before and after each bout. For the ensemble trajectories, we averaged the activity of all units within each ensemble and concatenated each ensemble's activity during all grooming bouts including 5 seconds before and after each bout. We initialized using Nonnegative Double Singular Value Decomposition with zeros filled with the average neural activity and minimized the Frobenius norm of the loss.

Statistical methods
All statistical tests were two-sided. Statistical significance was set to p = 0.05. Summary data in all figures are reported as mean ± SEM.

Software
All custom software was written in Python 95 unless stated otherwise. In addition to those mentioned elsewhere, we used the following Python packages: numpy 96 , scipy 97 , matplotlib 98 , scikit-learn 99 , and pandas 100,101 .
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 23, 2023. ; https://doi.org/10.1101/2023.02.23.529669 doi: bioRxiv preprint