Higher-order thalamic input to cortex selectively conveys state information

Cell counts and fluorescence intensity quantification in histological sections

Fluorescently labeled cells were counted manually in each histological section. Sections were aligned to bregma based on the detection of a small volume (<50 nL) of red retrobeads (Lumafluor) injected into the brain at the anterior-posterior coordinate of bregma immediately before transcardial perfusions as a fiducial marker. Fluorescently labeled cells were assigned to specific brain regions based on alignment with coronal sections in the Paxinos-Franklin mouse brain atlas and previously reported coordinates.^96,97 To calculate the total cell density in each brain region, the total number of cells counted in that region was divided by an estimate of its total volume, computed by numerical integration of the cross-sectional areas of the portions of each section corresponding to the brain region. Normalized cell count as a function of anterior-posterior distance from bregma was calculated by dividing the cell count associated with a given histological section by the highest cell count among all the sections from a given animal. To calculate the relative proportions of projections to PM for each mouse, we put the total number of tdTomato-positive cells counted across all areas in the denominator and the total number of tdTomato-positive cells counted in V1 and LM in the numerator and multiplied by 100 to yield the percentage that V1 and LM contribute to the feedforward projections to PM. Fluorescence intensity of axonal or dendritic labeling was measured in 10- or 5-μm bins across the laminar depth of the cortex relative to the pia for whole-cortex and layer 1 measurements, respectively. Mean normalized laminar fluorescence intensity for a given animal was calculated by dividing all fluorescence intensities by the highest fluorescence intensity for that animal and then averaging resulting normalized fluorescence vs. depth data from all sections in that animal.

Wheel position and change-points

Wheel position was determined from the output of the linear angle detector. The circular wheel position variable was first transformed to the [−π, π] interval. The phases were then circularly unwrapped to yield running distance as a linear variable, and locomotion speed was computed as a differential of distance (cm/s). A change-point detection algorithm detected locomotion onset/offset times based on changes in standard deviation of speed. Locomotion onset or offset times were estimated from periods when the moving standard deviations, as determined in a 0.5-s window, exceeded or fell below an empirical threshold of 0.1. Locomotion trials were required to have average speed exceeding 0.5 cm/s and last longer than 1 s. Quiescence trials were required to last longer than 2 s and have an average speed <0.5 cm/s.

Quantification of calcium signals

Analysis of imaging data was performed using ImageJ and custom routines in MATLAB (Mathworks). Motion artifacts and drifts in the Ca²⁺ signal were corrected with the moco plug-in in ImageJ,⁹⁸ and regions of interest (ROIs) were selected as previously described.^{31,72,99–101} All pixels in each ROI were averaged as a measure of fluorescence, and the neuropil signal was subtracted. ΔF/F was calculated as (F-F₀)/F₀, where F₀ was the lowest 10% of values from the neuropil-subtracted trace for each session.

For some analyses (Figure 4; Figure S6), ΔF/F signals were deconvolved to estimate the times of underlying Ca²⁺ events. The AR1 FOOPSI algorithm in the Oasis toolbox¹⁰² was used to find the optimal convolution kernel, baseline fluorescence, and noise distribution. Ca²⁺ events were then detected as times at which the deconvolved signal increased beyond 3 standard deviations. Time-varying Ca²⁺ event rates were then computed as

where $t_j$ indicates the time at which a Ca²⁺ event occurred for $n$ events in a recording, and $w$ is a sliding Gaussian window,

where ${\sigma }_{\mathrm{w}}=100\mathrm{m}\mathrm{s}$.

Cellular ROI selection

Image series collected from 2-photon calcium imaging were corrected for motion artifacts using the moco plug-in on ImageJ.⁹⁸ ROIs corresponding to cell bodies or axon terminals within the imaging field of view were chosen with a semi-automated procedure using mean- or max-projections of the imaging series.^99–101,103 All pixels within an ROI were averaged as a measure of fluorescence. To correct for fluorescence originating from neuropil potentially contaminating the ROI signal, a group of pixels corresponding to a ~20 μm region around each ROI (excluding all selected ROIs) was averaged and subtracted from the ROI signal at each time point:

From these neuropil-corrected raw fluorescence traces for each ROI, calcium activity was quantified as $\mathrm{\Delta }F/F$:

where $F_0(t)$ is a sliding 10^th percentile of the raw fluorescence distribution, with a sliding window 60 s in width.

For ROIs corresponding to axon terminals, additional efforts were made to minimize inclusion of redundant ROIs that putatively belong to the same neuron in the presynaptic source region. Whereas redundancies in cell body ROIs are highly unlikely (i.e., 2 or more ROIs corresponding to the same neuronal cell body), ROIs for axon terminals could suffer from redundancy due to axonal branching. In order to minimize the number of redundant axon terminal ROIs in our dataset, we used a procedure based on activity correlations. For this procedure, we first generated a dataset derived from axonal calcium imaging sessions consisting of axon segments that were unambiguously part of the same stretch of axon (“sub-ROIs”). For all such sub-ROI pairs within an imaging field of view, we then calculated the partial correlation of their activity throughout the entire imaging session, subtracting the mean value of all pairwise correlation values from each raw pairwise correlation. The resulting distribution of partial correlation values corresponds to the case in which all axon terminal ROIs in a dataset are actually from the same neuron. We then compared this distribution of sub-ROI partial correlation values with that of the partial correlation values between axon terminal ROIs that were considered to belong to different neurons based upon structural images (i.e., terminals that were non-contiguous within the field of view). Based on an analysis of the distributions of correlations for ROIs and sub-ROIs (Figure S1), a partial correlation value of 0.3 was the most conservative threshold for separating these distributions in our datasets. Thus, in our larger dataset of axon terminal ROIs, if any 2 ROIs exhibited partial correlations ≥0.3, one of the ROIs was discarded as putatively redundant and not included in further analysis. This method was also applied to eliminate putatively redundant ROIs from experiments in which we imaged apical dendrites from pyramidal neurons arborizing in layer 1 (Figure S1; Figure S4).

Quantification of visual responses and visual feature selectivity

Prior to imaging sessions used for characterizing neuronal visual responses, the population receptive field for the ROIs in a field of view was mapped by partitioning the LCD monitor into a 3×3 grid and presenting binarized Gaussian noise stimuli¹⁰⁴ within each patch in the grid. In subsequent imaging experiments using a given field of view, the monitor was positioned such that it was centered at the location within the visual field that elicited the largest visual responses from the neuron ROIs during these initial receptive field mapping sessions. With the exception of size tuning and motion coherence protocols (see below), visual stimuli were ~40° in diameter.

For ROIs corresponding to cell bodies or axon terminals, visually evoked response magnitude was quantified as the mean z-scored $\mathrm{\Delta }F/F$ during a 2-s visual stimulus, with the mean and standard deviation of the $\mathrm{\Delta }F/F$ signal during a 1-s baseline period prior to visual stimulus onset used for z-scoring. The overall visual response magnitude for a given ROI was the mean visually evoked z-scored $\mathrm{\Delta }F/F$ averaged across all visual stimulus presentations during an imaging session. ROIs were considered significantly visually responsive if $\mathrm{\Delta }F/F$ values during visual stimulation periods were significantly larger than $\mathrm{\Delta }F/F$ values during 1-s periods before visual stimulus onset (p < 0.05, Student’s t-test). Calculation of selectivity for different visual features is described below. Visual trials were not separated by behavioral state.

Quantification of modulation of neuronal activity by locomotion

For calculation of locomotion modulation index (LMI), only locomotion bouts occurring after and before at least 15 s of quiescence were used. LMI for each ROI was calculated as

where ${Ca}_{L\mathit{-on}}$ is the mean $\mathrm{\Delta }F/F$ during the 4 s following locomotion onset and ${Ca}_Q$ is the mean $\mathrm{\Delta }F/F$ during the 10–15 s after locomotion offset.

Quantification of modulation of neuronal activity by pupil size

From video of the mouse’s face, pupil size of the non-visually-stimulated eye was computed from individual video frames. Frames were first cropped to a region around the eye and image pixels were binarized such that the majority of the pixels within the pupil were white and the majority of the pixels outside of the pupil were black (the IR light from the Ti-sapphire laser that escapes through the mouses eye causes the pupil to appear white when captured by the camera). Canny edge detection was then used on the binarized image frames to detect edges of the pupil. Pupil size was then estimated as the authalic diameter of the area enclosed by the detected edges (i.e., the diameter of a circle with the same area).⁴⁴

All analyses relating neuronal activity to pupil fluctuations were restricted to quiescent (i.e., non-locomotion) periods of imaging sessions. Visual stimuli were randomly presented throughout the session and these periods were not excluded from the overall analysis. To relate pupil dilation and constriction to fluctuations in neuronal activity, the time series of pupil size was bandpass-filtered between 0.1 and 1 Hz. The bandpass-filtered signal of pupil size was converted to phase of the pupil dilation-constriction cycle through a Hilbert transform, which transforms the signal into a time series of complex numbers. The inverse tangent of the ratio between the imaginary and real parts of the complex number at each time point yields the phase of the bandpass-filtered pupil signal. To quantify neuronal activity as a function of the phase of the dilation-constriction cycle, $\mathrm{\Delta }F/F$ time series for each ROI were also band-pass-filtered between 0.1 and 1 Hz. These bandpass-filtered $\mathrm{\Delta }F/F$ signals were then z-scored using baseline periods associated with minimal changes in pupil size and small baseline pupil diameter. Specifically, these baseline periods needed to be associated with a mean of the bandpass-filtered pupil signal that was less than 1.5x the interquartile range of this signal and a mean pupil size that was within the bottom 25^th percentile of the distribution of pupil sizes observed during quiescent periods during a given session. The mean and standard deviation of the bandpass-filtered $\mathrm{\Delta }F/F$ signal during these baseline periods was used for z-scoring the $\mathrm{\Delta }F/F$ signal. For each dilation-constriction cycle, z-scored $\mathrm{\Delta }F/F$ for each ROI was then averaged in bins of width π/32 from −π to π. Only dilation-constriction periods larger than a threshold of 1.5x the interquartile range of the bandpass-filtered pupil signal were analyzed.

To calculate cross-correlations between neuronal activity and pupil fluctuations, both the time series of pupil size and the time series of $\mathrm{\Delta }F/F$ for each ROI were lowpass-filtered at 10 Hz. Correlations between these signals were then computed at different lags in 100-ms increments in a time window of 8 s, using the pupil signal as the reference signal. The strength of the modulation of $\mathrm{\Delta }F/F$ by pupil size for each ROI was quantified as the value of the peak cross-correlation between $\mathrm{\Delta }F/F$ and pupil size.

Quantification of modulation of neuronal activity by facial motion/whisking

From video of the mouse’s face, a region consisting of the nose, mouth, and whisker pads was cropped from each frame. A region of interest was then drawn around the whisker pad and facial motion energy was calculated as the absolute value of the difference in average pixel intensity of the whisker pad between successive image frames. Onset times of facial motion bouts were computed by first z-scoring the entire time series of facial motion energy. High and low thresholds of facial motion energy were then defined as the 60% and 40% quantiles of the z-scored facial motion energy distribution during quiescent periods. After smoothing the time series of z-scored facial motion energy with a 1-s moving-average filter, high and low facial movement states were defined as periods in which the smoothed signal remained above or below the high and low threshold levels for at least 500 ms. Onset of a facial motion bout was then defined as the time at which a high facial motion state was preceded by at least 4 s of a low facial motion state. ROI calcium activity aligned to these onsets was z-scored using the mean and standard deviation of the $\mathrm{\Delta }F/F$ during low facial motion periods. As for pupil diameter, facial motion data were measured during imaging sessions where visual stimuli were randomly presented.

To calculate cross-correlations between neuronal activity and facial motion, both the time series of facial motion energy and the time series of $\mathrm{\Delta }F/F$ for each ROI were lowpass-filtered at 10 Hz. Correlations between these signals were then computed at different lags in 100-ms increments in a time window of 8 s, using the facial motion energy as the reference signal. The strength of the modulation of $\mathrm{\Delta }F/F$ by facial motion for each ROI was quantified as the value of the peak cross-correlation between $\mathrm{\Delta }F/F$ and facial motion energy.

Calculation of Ca²⁺ event-triggered averages

For data shown in Figures 4E and 4F and Figures S6B and S6C, Ca²⁺ activity aligned to cell body Ca²⁺ events was computed as follows. For each riboGCaMP6m-expressing cell body identified in a field of view, the times of its Ca²⁺ events were estimated from the deconvolved ΔF/F signal. In an 8-s window around each Ca²⁺ event from the cell body, the population average Ca²⁺ event rate from all other relevant ROIs (i.e., axon ROIs or ROIs from other cell bodies) was calculated. For each cell body whose event-triggered average was calculated this way, the peak event rate was taken as the maximum event rate in the 8-s window.

Statistical tests

In both the text and figures “N” refers to number of mice and “n” refers to number of neuron ROIs (cell bodies or axon terminals). Unless otherwise noted, data are summarized by either box-and-whisker charts (with “boxes” marking the median value and the 25% and 75% quantiles, and “whiskers” marking the minimum and maximum values in the dataset) or summary curves indicating the mean value and shading indicating the standard error of the mean.

For statistical comparisons among groups, we accounted for the nested structure of the data (i.e., cellular ROIs nested within mice). Without taking nesting into consideration, false positive rates are inflated since cellular ROIs from the same mouse cannot be considered independent measurements.¹⁰⁷ We used semi-weighted error estimators to account for the nesting of cellular ROIs within mice.⁷² For this procedure, let $y_i$ be the mean for a given parameter (e.g., locomotion modulation index) for the $i$-th mouse in a group, where $x_j$ is the parameter value for the $j$-th ROI in a set if $M$ ROIs for mouse $i$:

The unweighted estimator of the mean is defined as the mean of $y_i$ for all mice in a group of $N$ mice, where $N$ is the number of degrees of freedom. This analysis is suboptimal, as some animals have many cells and yield more reliable estimates, whereas other animals have few cells and yield less reliable estimates. At the other extreme, the weighted estimator is defined by pooling parameter values from all ROIs across the $N$ mice in a group and using this number as the degrees of freedom. While the weighted estimator is commonly used in statistical comparisons, the false positive rate is improperly controlled in this case. The semi-weighted estimator for a given parameter measured from multiple ROIs from multiple mice in a group is defined as:

where the weight for the $i$-th mouse is defined as

and $s_i^2$ is the within-mouse variance across ROIs and ${\tau }^2$ is the estimated variance across all mice within the group using maximum likelihood estimation¹⁰⁸:

Equations 12 and 13 were solved iteratively and checked for convergence of ${\tau }^2$ to <0.001. If there is minimal across-mouse variance, then the contribution of each ROI is weighted by the number of mice (the weighted estimator), whereas if across-mouse variance is large (i.e., parameter measurements for ROIs within a mouse are strongly dependent), then each mouse is given the same weight.

The semi-weighted estimators derived for different cell types (e.g., V1 axons vs. LP axons) or experimental conditions (e.g., control vs. eOPN3) were then used to calculate the statistical significance of pairwise comparisons using a Student’s t-test. The Benjamini-Hochberg procedure was used to control the false discovery rate from multiple pairwise comparisons. p-values <0.05 after multiple comparisons correction were considered statistically significant. Central markers on box-whisker plots indicate median values, boxes delineate the 25^th and 75^th percentiles, and upper/lower whiskers indicate maximum/minimum values.