Analysis of signaling during locomotion

Recordings included for locomotion analysis come from 12 Vglut2⁺ mice (four males, seven females, one unrecorded), six Calb1⁺ mice (five males, one femle), nine Anxa1⁺ mice (four males, five females), 14 Aldh1a1⁺ mice (five males, eight females, one unrecorded) and 14 DAT⁺ mice (three males, seven females, two unrecorded).

All times within a 5-s window after any stimulus was delivered (reward, air puff) were excluded from all locomotion analysis described below; this excluded reward/air puff movement reactions and consumptive licking behavior from the locomotion analysis.

Only locomotion time bins were included for locomotion analysis in Fig. ​Fig.22 and Extended Data Figs. ​Figs.1,1, ​,6,6, ​,77 and ​and9.9. Locomotion versus rest bins were selected using a double threshold on the velocity trace in both positive and negative directions (thresh1 = ±0.024 m s⁻¹ and thresh2 = ±0.010 m s⁻¹). Isolated one-bin-long locomotion periods (no other movement within two bins on either side) were excluded as well as rest periods shorter than 0.5 s. Time bins were considered as locomotion periods only if they lasted longer than 0.5 s and had an average velocity greater than 0.2 m s⁻¹. For a recording to be included in the locomotion analysis, the recording needed to include a total of at least 100 s of locomotion.

Acceleration was calculated from the velocity traces as the difference between consecutive treadmill velocity time bins (first smoothed over six bins, 0.06 s) and then multiplied by the sampling frequency (100 Hz) for proper m s⁻² units.

Cross-correlations between ΔF/F and acceleration (Figs. 2d,f and ​and3a3a and Extended Data Figs. 1e,f, 6c,h,i and 9h,g) were calculated for locomotion periods only (defined above) using MATLAB’s crosscorr function over a 1-s lag window (100 time bins). The same process was used to calculate the cross-correlation between corresponding 405-nm ΔF/F traces and acceleration, and any recording with a peak cross-correlation (between 405-nm ΔF/F trace and acceleration) above 0.1 was excluded from all locomotion analysis. This same strategy was also used to calculate cross-correlations between different variables (licking versus velocity, licking versus acceleration, licking versus ΔF/F and velocity versus ΔF/F) as shown in Extended Data Fig. ​Fig.7h.7h. Averages per mouse (Extended Data Fig. ​Fig.6i)6i) were obtained by averaging together the cross-correlations for all recordings made from the same mouse.

For triggered averages of ΔF/F on accelerations and decelerations (Figs. ​(Figs.2g2g and ​and5a5a and Extended Data Figs. ​Figs.6b,6b, ​,7d7d and ​and9j),9j), isolated large accelerations and decelerations were selected by first locating the zero-crossings on the acceleration trace (points where the acceleration trace crosses zero, going from negative to positive or vice versa), considering individual accelerations/decelerations the interval between two zero-crossings of the trace. Accelerations/decelerations were included if they had a duration of at least 50 ms (0.05 s) and a peak greater than 2 m s⁻² (accelerations) or lower than −2 m s⁻² (decelerations), but only if they were not surrounded by other large accelerations or decelerations (no acceleration >2 m s⁻² or <−2 m s⁻² in a window of 0.25 s on either side). Triggered averages were the result of averaging 135 ± 77 accelerations and 138 ± 101 decelerations (mean ± s.d.) per recording (Fig. ​(Fig.2g).2g). For an example recording of each subtype showing that recording’s triggered average with a heat map of each individual event that contributes to the recording’s average, see Extended Data Fig. ​Fig.7d.7d. The probability that ΔF/F transients follow accelerations is 57.5% for Anxa1, and the probability that transients follow decelerations is 62.4% for Vglut2⁺ and 62.3% for Calb1⁺, as shown in Extended Data Fig. ​Fig.7f.7f. This was obtained by first calculating for each event the integral of the ΔF/F trace within a 0.75-s window from the start of the acceleration/deceleration (t = 0 s), after subtracting the ΔF/F value at t = 0 s. Histograms were then obtained for the percent of accelerations/decelerations per recording with positive values for this calculation.

Conversely, for triggered averages of acceleration on ΔF/F transient peaks (Fig. ​(Fig.2h2h and Extended Data Figs. ​Figs.6d,6d, ​,7e7e and 9i), we selected well-isolated transients from non-normalized ΔF/F traces, as defined by having a large, fast rise (30 ΔF/F s⁻¹) immediately followed by a decay (as used in the calculation of signal-to-noise ratio above). Triggered averages were the result of averaging 423 ± 271 transients (mean ± s.d.) per recording (Fig. ​(Fig.2h).2h). For an example recording of each subtype showing that recording’s triggered average with a heat map of each individual event that contributes to the recording’s average, see Extended Data Fig. ​Fig.7e7e.

For triggered averages on movement onsets and offsets (Extended Data Fig. ​Fig.7g),7g), we started with the transitions between locomotion and rest bins as selected above. Often mice moved backwards or jittered before starting to run or after stopping, so to select only clean onsets we only included transitions that reached a velocity of 0.4 m s⁻¹ within 0.75 s of starting to move, with an initial acceleration peak of at least 1 m s⁻² and with no negative velocities below –0.05 m s⁻¹ before this strong acceleration. For offsets, the symmetric conditions were required (stopping from a velocity of at least 0.4 m s⁻¹ within 0.75 s, with a final deceleration of at least –1 m s⁻² and no negative velocities below –0.05 m s⁻¹ at the end of the offset). For plotting of cross-correlation and triggered averages above, traces were smoothed over five time-lag bins (0.05 s). Shaded areas represent the mean ± s.e.m. across recordings, and accompanying heat maps show cross-correlations/triggered averages for all individual recordings. Heat maps in Extended Data Fig. ​Fig.1f1f were sorted by the integral of the ΔF/F-acceleration cross-correlation at positive lags (see below), whereas heat maps in Fig. 2f–h and Extended Data Figs. 6b–d and 7g–i were sorted by PC1/PC2 angle (see PCA subsection below)—other than that, the data plotted in Extended Data Fig. ​Fig.1f1f and in Fig. ​Fig.2f2f and Extended Data Fig. ​Fig.6c6c are the same (with the addition of Anxa1⁺).

SNc recordings (Extended Data Fig. 9h–k) were analyzed in the same manner as striatal recordings. Recordings included for locomotion analysis in SNc come from 11 Vglut2⁺ mice (five males, six females), three Calb1⁺ mice (three males), eight Anxa1⁺ mice (four males, four females), 13 Aldh1a1⁺ mice (five males, seven females, one unrecorded) and eight DAT⁺ mice (two males, three females, two unrecorded).

For the initial functional characterization shown in Extended Data Fig. 1h,i, differences in locomotion signaling were quantified by calculating the integral of the cross-correlation between ΔF/F and acceleration at positive lags (0–1 s), where positive values indicate a peak in the cross-correlation and, thus, ΔF/F transients after accelerations, whereas negative values indicate a trough and, thus, ΔF/F transients after decelerations. For the quantification of acceleration/deceleration signaling across depths in striatum shown in Extended Data Fig. ​Fig.1h,1h, depth from surface was defined as the depth at which the fiber tip was located from the brain surface, as measured by the micromanipulator used to move the fiber during photometry. To reduce overlap between data points at the same depth plotted, a random amount between +0.1 mm and −0.1 mm was added to each depth. This measure of locomotion signaling was also used to plot the relationship between locomotion signaling and reward responses in Extended Data Fig. ​Fig.1i1i (for reward response calculation, see ‘Analysis of responses to rewards and air puffs’ subsection below) and to sort the ΔF/F-acceleration correlation plots in Extended Data Fig. ​Fig.1f1f.

For analysis of timing differences between Calb1⁺ and Vglut2⁺ deceleration signaling shown in Fig. ​Fig.2i2i and Extended Data Fig. 6f,g, the lag between ΔF/F transient peaks and deceleration peaks was quantified by locating in time the minimum cross-correlated value between 0 s and 1 s for the ΔF/F-acceleration cross-correlations for each recording (Fig. ​(Fig.2i),2i), the maximum ΔF/F value between 0 s and 1 s for the triggered average on deceleration (Extended Data Fig. ​Fig.6f)6f) or the minimum acceleration value between −1 s and 0 s for the triggered average on transient peaks (Extended Data Fig. ​Fig.6g6g).

For calculating the relationship between velocity and ΔF/F as shown in Extended Data Fig. ​Fig.7c,7c, we divided the velocity and ΔF/F traces based on the velocity at each timepoint into bins of 0.1 m s⁻¹ ((−0.05:0.1:0.75 inf)) and averaged the ΔF/F for each subtype and bin.

For checking whether the locomotion signaling observed in DAT mice across depths could be explained by mixtures of the Anxa1⁺ and Calb1⁺ subtypes (Extended Data Fig. 6h,h′), we first divided DAT recordings made in the anterior striatum (anterior to bregma) based on the depth from the brain surface at which they were made, from 1.5 mm to 4 mm in 0.5-mm bins, and obtained the average cross-correlation between ΔF/F and acceleration for each subset (H), as explained above. We then calculated weighted averages between the average cross-correlations for the Calb1+ and Anxa1⁺ subtypes in different ratios to match the approximate relative abundance of each subtype’s axons across depths: from 100% Anxa1⁺, 0% Calb1⁺ for dorsal striatum to 0% Anxa1⁺, 100% Calb1⁺ for ventral striatum (H′).

For determining whether the size of the ΔF/F response scaled with the size of the acceleration/deceleration (Fig. 5a,b), for each recording we divided all the accelerations/decelerations that fulfilled the requirements described above (for triggered averages on accelerations and decelerations) into five quartiles per recording based on the peak acceleration/deceleration and calculated the average acceleration (Fig. ​(Fig.5a,5a, left) and ΔF/F (Fig. ​(Fig.5a,5a, right) triggered on accelerations/decelerations within each of these five quartiles. For plotting the transient amplitude for each of these acceleration/deceleration quantiles (Fig. ​(Fig.5b),5b), we calculated the difference between the ΔF/F value at t = 0 (trigger point, start of the acceleration/deceleration) and the maximum ΔF/F value within the 1-s window following it. The fold increase as reported in the legend was calculated by dividing the transient amplitude for the largest deceleration/acceleration quintile by the transient amplitude for the smallest deceleration/acceleration quintile. Statistical significance was calculated using a paired Wilcoxon signed-rank test with Bonferroni correction (P values multiplied by 3) comparing the transient amplitude for the smallest versus largest acceleration/deceleration quintiles for each recording.

PCA was applied to the matrix of all cross-correlation traces from striatal recordings (shown in Fig. ​Fig.2f),2f), from all functionally homogeneous subtypes (Vglut2⁺, Calb1⁺ and Anxa1⁺), using MATLAB’s pca function without centering: ‘centered’, ‘off’. Centering was not used so as to maintain the cross-correlation values’ relationship to 0 and to avoid biasing the results based on the relative number of recordings from different subtypes; however, equivalent results were obtained when we repeated the PCA analysis with centering (Extended Data Fig. ​Fig.7a).7a). This pca function outputs the PCs (loadings and eigenvectors), the scores for each recording’s cross-correlation along each PC (matrix of all SNc cross-correlation traces multiplied by the loadings matrix) and the variance explained by each PC across all recordings. For the representation of combinations of the first two PCs (PC1 and PC2) shown in Fig. ​Fig.2j,2j, PC1 and PC2 were weighted by the s.d. of their scores across recordings (~1 for PC1 and ~0.7 for PC2), to accurately represent each quadrant in Fig. 2k,l and Extended Data Figs. ​Figs.6e6e and 9k. Figure ​Figure2k2k and Extended Data Fig. ​Fig.6e6e show the PC1 and PC2 scores for each recording of each subtype. In Extended Data Fig. ​Fig.6e,6e, recordings were color-coded based on the depth from brain surface at which they were recorded, as measured by the micromanipulator used to move the fiber during photometry. For Extended Data Fig. ​Fig.7b,7b, the cross-correlation traces were first normalized by dividing the trace by its absolute maximum value with its sign, so that its lowest value was –1 or its maximum value was +1 while maintaining 0, before PCA analysis. Because PC1 and PC2 explain most of the variance, this results in data points being pushed to a ring around the origin.

For SNc recordings, the cross-correlations between ΔF/F and acceleration for all recordings of all subtypes, as shown in Extended Data Fig. ​Fig.9k,9k, were decomposed using the same PCs calculated above from the striatal cross-correlations. Scores for SNc cross-correlations (Extended Data Fig. ​Fig.9h)9h) were calculated by multiplying the matrix of all SNc cross-correlation traces by the striatal loadings matrix (PCs). The percent of SNc variance explained by each PC (PC1 = 53.2% of variance, PC2 = 24.3%) was calculated as the variance without the mean subtracted (not centered).

In the PC1/PC2 space shown in Fig. ​Fig.2k,2k, the angle of each point from the origin represents the shape of the cross-correlation between acceleration and ΔF/F and, thus, the different relationships between subtypes’ signaling and acceleration, whereas the distance from the origin represents the amplitude of the cross-correlation. To quantify the shape of the cross-correlation across subtypes, we calculated the angle of each recording in the PC1/PC2 space (with each PC weighted by its s.d.) and plotted it in a radial histogram (Fig. ​(Fig.2l).2l). This angle was also used for plotting of subtypes in Fig. 6a,b and Extended Data Fig. 9f,g. All angles in this paper are reported as standard with 0° set between quadrants I and IV and angles increasing in the counterclockwise direction (that is, up is 90°). P values for reporting statistical significance for the difference between subtypes across this angle PC1/PC2 space were calculating by opening the angular space at 45° (the region where the least recordings from Calb1⁺/Vglut2⁺/Anxa1⁺ fall) and using a Wilcoxon rank-sum test with Bonferroni correction (multiply P values by 3) to compare subtypes. This angle was also used to sort cross-correlation and triggered average heat maps in Fig. 2f–h and Extended Data Figs. 6b–d and 9h–j, starting by the middle of the quadrant opposite to the center of mass for each subtype (315° for Vglut2⁺, 45° for Calb1⁺ and DAT⁺ and 135° for Anxa1⁺ and Aldh1a1⁺) and going counterclockwise. Figure ​Figure3b3b and Extended Data Fig. ​Fig.6a6a show the anatomical location of each recording color-coded based on the PC1/PC2 angle and distance from the origin for that recording. The colormap was defined by assigning a different color to the middle of each quadrant (45°, 135°, 225° and 315°), where the center of mass of each subtype approximately falls at, and then fading that color to white as the values of PC1 and PC2 decrease to 0. In Fig. ​Fig.3b3b (but not in Extended Data Fig. ​Fig.6a),6a), recording locations were collapsed into a single brain slice for anterior striatum and another for posterior striatum, and locations were shifted a random amount between ±0.4 mm mediolaterally for visibility. For details on how the x–y–z coordinates for each recording were calculated, see the ‘Fiber placement localization’ subsection.

To calculate the difference in locomotion signaling between pairs of recordings based on the distance between them, as shown in Extended Data Fig. ​Fig.6j,6j, we used the difference between the PC1/PC2 angles of each pair of recordings calculated as above (maximum angle difference is 180°). For the distance between each pair’s recording locations, we used the Euclidian distance between the x–y–z coordinates of the recordings, obtained as described in the ‘Fiber placement localization’ subsection of these methods below. For within-subtype comparisons (Calb1⁺, Anxa1⁺ and Vglut2⁺) and for DAT, all recordings for that subtype/DAT were compared with all other recordings from that same subtype/DAT. For the mismatch–subtype comparisons, each recording from Calb1⁺, Anxa1⁺ and Vglut2⁺ was compared to all recordings from the other two subtypes (for example, each Calb1⁺ recording was compared to each Anxa1⁺ and Vglut2⁺ recording). Statistical significance was calculated using a Mann–Whitney U-test with Bonferroni correction (P values multiplied by 21, the total number of comparisons performed) comparing all the mismatch–subtype pairs with all the Vglut2⁺, Calb1⁺ or Anxa1⁺ pairs.