Thalamus-driven functional populations in frontal cortex support decision-making

Mice

This study was based on data from 186 mice (age > postnatal day 60, both male and female mice). Sixteen VGAT-ChR2-EYFP mice (JAX 014548), four PV-ires-cre (JAX 008069)59 crossed with a red-shifted channelrhodopsin (ReaChR) reporter mice (Rosa26-LSL-ReaChR-mCitrine, JAX 026294)60 and one PV-ires-cre mouse crossed with Ai32 (Rosa26-LSL-ChR2-EYFP, JAX 012569) were used for electrophysiology and photoinhibition during behavior experiments. Seven GtACR1 reporter mice (R26-LNL-GtACR1-Fred-Kv2.1, JAX 033089) were used for ThalALM photoinhibition44. Seventeen Ai32 mice were used for in vivo whole-cell recording to characterize connectivity from vS1 and secondary motor cortex (M2) to vM1. Twenty-eight ReaChR mice were used for in vivo whole-cell recording to characterize S1/S2, cALM and ThalALM connectivity to ALM in untrained passive mice. Sixteen ReaChR and 11 Ai32 mice were used for in vivo whole-cell recording of S1/S2, cALM and ThalALM connectivity to ALM during tactile decision-making behavior. Ten ReaChR mice were used for silicon probe recording and ChR2-tagging to characterize S1/S2, cALM and ThalALM connectivity to ALM during behavior experiments. Finally, five wildtype mice were used for anatomical tracing.

We analyzed three silicon probe recording datasets previously collected in the same task conditions (from refs. 14,24,31). Combined, the primary extracellular recording dataset (1.3-s delay epoch) included new data from 31 mice (described above, 4,967 units) and reused data from 42 mice (4,659 units)14,24. The second extracellular recording dataset (1.7-s delay epoch) included reused data from refs. 24,31 (29 mice, 10,420 neurons).

All procedures were in accordance with protocols approved by the Institutional Animal Care and Use Committees at Baylor College of Medicine. Mice were housed at a constant temperature (22 ± 1 °C) and humidity (30–55%) under a 12:12 reversed light/dark cycle and tested during the dark phase. On days not tested, mice received 0.5–1 ml of water. On other days, mice were tested in experimental sessions lasting 1–2 h where they received all their water (0.3–1 ml). If mice did not maintain a stable body weight, they received supplementary water61. All surgical procedures were carried out aseptically under 1–2% isoflurane anesthesia. Buprenorphine Sustained Release (1 mg kg−1) and Meloxicam Sustained Release (4 mg kg−1) were used for preoperative and postoperative analgesia. A mixture of bupivacaine and lidocaine was administered topically before scalp removal. After surgery, mice were allowed to recover for at least 3 d with free access to water before water restriction.

Surgery

Mice were prepared with a clear skull-cap and a headpost26,61. The scalp and periosteum over the dorsal skull were removed. A layer of cyanoacrylate adhesive (Krazy glue, Elmer) was directly applied to the intact skull. A custom headpost was placed on the skull over the visual cortex and cemented in place with clear dental acrylic (Lang Dental Jet Repair Acrylic; Part no. 1223-clear). A thin layer of clear dental acrylic was applied over the cyanoacrylate adhesive covering the entire exposed skull, followed by a thin layer of clear nail polish (Electron Microscopy Sciences, Part no. 72180). For ThalALM photoinhibition, a 5-mm optic fiber (Thorlabs, Part no. CFMLC12L05) was implanted above the left ThalALM25. For cALM or S1/S2 photoinhibition, a plastic fitting was glued onto the clear skull implant above the right ALM or left S1/S2 for attachment of the optic fiber.

Behavior

The behavioral task and training have been described61,62. The stimulus was a metal pin (0.9 mm in diameter), presented at one of two possible positions (Fig. 1a). The two pole positions were 5 mm apart along the anterior–posterior axis. The posterior pole position was 5 mm from the whisker pad. A two-spout lickport (4.5 mm between spouts) was used to deliver water rewards and record licks.

At the beginning of each trial, the vertical pole moved into reach of the whiskers (0.2-s travel time), where it remained for 1 s, after which it was retracted (retraction time 0.2 s). The sample epoch was defined as the time between the pole movement onset to 0.1 s after the pole retraction onset (sample epoch, 1.3 s; Fig. 1a). Mice touched the object at both pole positions, typically with a different set of whiskers. The delay epoch (1.3 s for primary dataset, 1.7 s for second datasets) followed the sample epoch. An auditory ‘go’ cue indicated the end of the delay epoch (pure tone, 3.4 kHz, 0.1 s). Licking early during the trial was punished by a loud alarm sound (0.05 s) and a brief timeout (1–1.2 s). Licking the correct lickport after the ‘go’ cue led to a water reward (2–3 µl). Licking the incorrect lickport triggered a timeout (2–6 s). Trials in which mice did not lick within a 1.5-s window after the ‘go’ cue (‘ignore’) were rare and typically occurred at the end of a session. Reaction time was from the ‘go’ cue onset to the first lickport contact.

Viral injection and histology

Glass pipettes (20–30-µm-diameter tip and beveled) were back-filled with mineral oil and front-loaded with viral suspension immediately before injection.

For anatomical tracing, AAV2.CAG.GFP (Addgene, 37825) was injected in the left ALM (anterior 2.5 mm from lambda, lateral 1.5 mm, depth 0.5 and 0.8 mm, 100 nl at each depth). At 14 d post injection, WGA (Thermo Fisher Scientific, WGA-Alexa594, 2% in PBS, 200 nl) was injected in the same location and incubated for 24 h. Mice were perfused transcardially with PBS followed by 4% paraformaldehyde (PFA)/0.1 M PBS. The brains were fixed overnight in 4% PFA and transferred to 30% sucrose before sectioning on a cryostat (Thermo Scientific, Cryostar NX70). Coronal 30-µm free-floating sections were mounted with mounting medium with DAPI (Vector Laboratories, H-1500-10), imaged on a fluorescence macroscope (Olympus MVX10) and processed in ImageJ.

To characterize long-range input connectivity in ALM, we injected AAV9.CamKII.HI.eGFP-Cre.WPRE.SV40 in ReaChR or Ai32 mice in the right ALM (anterior 2.5 mm from bregma, lateral 1.5 mm, depth 0.5 and 0.8 mm, 100 nl at each depth), left ThalALM (posterior 1.5 mm, lateral 0.8 mm, depth 4.1 mm, 150 nl) or left S1/S2. To target a region spanning vS1 and S2, the left hemisphere was tilted down by 5° from the horizontal plane and the injection pipettes entered the brain vertically at posterior 1.5 mm and lateral 4 mm from bregma. Virus was injected at depths 0.8, 1.2 and 2 mm (100 nl at each depth). To characterize long-range input connectivity in vM1 (Extended Data Fig. 6), we injected AAV9.CamKII.HI.eGFP-Cre.WPRE.SV40 (Penn Vector Core, University of Pennsylvania) in Ai32 mice in vS1 (posterior 1.0 mm from bregma, lateral 3.1 mm, depth 0.5 and 0.8 mm, 100 nl at each depth) or M2 (anterior 2.7 mm, lateral 0.9 mm, depth 0.5 and 0.8 mm, 100 nl at each depth).

To quantify the fraction of anterogradely labeled neurons in ALM due to potential tropism of the Cre viruses (Extended Data Fig. 6e), we collected 30-µm free-floating coronal sections around ALM (three mice each for S1/S2, cALM and ThalALM injections). Sections were stained with NeuN. Regions covering ALM layers 2/3 and 5 were imaged with an LSM710 (Zeiss) and processed with ImageJ. Cell counting was performed manually (Extended Data Fig. 6e).

PhotostimulationChR2-tagging and ChR2-assisted circuit mapping

Photostimulation and electrophysiology recordings were performed in the left ALM to photostimulate ChR2- or ReaChR-expressing axons from left S1/S2, right ALM or left ThalALM. Light from a 473-nm (UltraLasers, MBL-FN-473-300 mW) or 593-nm laser (UltraLasers, MGL-N-593.5-200 mW) was controlled by an acousto-optical modulator (Quanta Tech, MTS110-A3-VIS), and focused onto the brain surface through a craniotomy (beam diameter: 400 µm at 4σ). For whole-cell recordings, photostimulation consisted of four powers (1, 5, 10, 20 mW) and four pulse conditions (1, 3, 5, 10 pulses; 2-ms pulses at 5-ms interval). For silicon probe recordings, photostimulation consisted of three powers (10, 20, 30 mW) and three light pulses (1-ms pulses at 200-ms interval). Photostimulation was tested outside of the behavioral task.

Photoinhibition of S1/S2, cALM and ThalALM

For photoinhibition of S1/S2 and cALM, we photostimulated GABAergic neurons in VGAT-ChR2-EYFP, PV-ires-cre × Ai32 or PV-ires-cre × ReaChR mice. These methods resulted in similar photoinhibition44. Light was delivered through an optic fiber placed on the clear skull implant (S1/S2, bregma posterior 1.5 mm, lateral 4 mm; cALM, anterior 2.5 mm, lateral 1.5 mm). We used 40-Hz photostimulation with a sinusoidal temporal profile. The duration was 1.3 s including a linear ramp during laser offset (100 or 200 ms). The average power was 4 mW. In a subset of S1/S2 photoinhibition, 8 mW was used. At this power, the photoinhibition silenced activity in a cortical area of 2-mm radius (at half-max) through all cortical layers44. For photoinhibition of ThalALM, we photostimulated the thalamic reticular nucleus axons in VGAT-ChR2-EYFP mice25. Photostimulation was through an optic fiber (Thorlabs, Part no. CFMLC22L05) implanted above ventral-medial nucleus/VAL (bregma posterior 1.5 mm, lateral 0.8 mm, depth 4.1 mm). The average power was 3 mW measured at the tip of the optic fiber. In these experiments, photostimulation occurred during either the sample or delay epoch randomly in 25% of trials.

We also directly photoinhibited ThalALM using soma-targeted GtACR1 (ref. 46). AAV9.CamKII.HI.eGFP-Cre.WPRE.SV40 (University of Pennsylvania Vector Core) was injected in ventral-medial nucleus/VAL (posterior 1.5 mm from bregma, lateral 0.8 mm, depth 4.1 mm, 120 nl) of R26-LNL-GtACR1-Fred-Kv2.1 mice44. In these experiments, photostimulation occurred during the sample, delay or response epoch randomly in 25% of trials. The average power was 1–3 mW.

To prevent mice from distinguishing photostimulation trials from control trials using visual cues, a masking flash (1-ms pulses at 10 Hz) was delivered using 470-nm or 590.6-nm light-emitting diodes (Luxeon Star) near the eyes of the mice. The masking flash began as the pole started to move and continued through the end of the epoch in which photostimulation could occur.

ElectrophysiologySilicon probe recordings

A craniotomy (diameter < 1 mm) was made over the left ALM. A silicon probe was acutely inserted 0.9–1.11 mm below the brain surface. To minimize brain movement, a drop of silicone gel (3-4680, Dow Corning) was applied over the craniotomy after the electrode was in the tissue. The tissue was allowed to settle for 15 min before the recording started. Extracellular spikes were recorded using 64-channel Cambridge NeuroTech silicon probes (H2 acute probe, 25-µm spacing, 2 shanks). The voltage signals were amplified and digitized on an Intan RHD2164 64-Channel Amplifier Board (Intan Technology) at 16 bits, recorded on an Intan RHD2000-Series Amplifier Evaluation System (sampling at 20,000 Hz) and stored for offline analysis. Two to eight recordings were made from each craniotomy. DiI was applied to the tip of the silicon probe in the last session to label the recording tracks.

Whole-cell recording

A craniotomy (diameter, 100–200 µm) was made in vM1 or ALM. Recordings were obtained using a glass pipette (tip resistance, 7–11 MΩ) and MultiClamp 700B amplifier (Molecular Devices). The signal was sampled at 20 kHz using Wavesurfer (http://wavesurfer.janelia.org/). Membrane potential was not corrected for liquid junction potential. The intracellular solution contained (in mM): 128 potassium gluconate, 4 MgCl2, 10 HEPES, 1 EGTA, 4 Na2ATP, 0.4 Na2GTP, 10 sodium phosphocreatine (pH 7.23; 283 mOsm). Aliquoted ATP/GTP was added to the internal solution on the day of recording. Positive pressure (200 mBar) was applied before insertion to reduce pipette tip contamination. Then, 1.5% agar (Sigma, A1296) in artificial cerebrospinal fluid (aCSF) was applied over the craniotomy after the pipette tip reached pia surface. Recording depth was based on manipulator reading. The series resistance was monitored through a current pulse (100 ms, −0.2 nA) injection. Only neurons with GΩ-seal were included for analysis. Once a GΩ-seal was achieved, an increasing negative pressure was applied slowly until break-in was established. A family of step currents (500 ms or 750 ms, in 40-pA steps) were injected in current-clamp mode (Extended Data Fig. 7). Each craniotomy was used for 1–2 recording sessions.

In calibration recordings, 1 mM TTX (Tocris Bioscience) and 100 µM 4-AP (Acros Organics) were applied topically over the recording craniotomy to verify synaptic connection with long-range input axons. TTX was a sodium channel blocker that prevented local action potential transmission and excitation in unconnected neurons (Fig. 5e,f and Extended Data Fig. 6a,b). Unlike ChR2-assisted circuit mapping in vitro, we found that 4-AP was not required to elicit light-induced EPSPs in connected neurons during application of TTX. To confirm the EPSPs arose from glutamatergic transmission, 20 mM NBQX (Tocris Bioscience) and 30 mM AP5 (Tocris Bioscience) were applied topically to block glutamate receptors. This abolished the EPSPs (Extended Data Fig. 6b).

Behavioral data analysis

We separately computed performance for ‘lick right’ and ‘lick left’ trials as the fraction of correct choices, excluding lick early trials and ignore trials (Fig. 4b). Significance of the performance change in each photostimulation condition was determined using a nested bootstrap to account for variability across mice, sessions and trials26. We tested against the null hypothesis that the performance change caused by photostimulation was due to normal behavioral variability. In each round of bootstrap, we replaced the original behavioral dataset with a resampled dataset in which we resampled with replacement from: (1) mice, (2) sessions performed by each mouse, (3) the trials within each session. We then computed the performance change on the resampled dataset. Repeating this procedure 10,000 times produced a distribution of performance changes that reflected the behavioral variability. The P value of the observed performance change was computed as the fraction of times the bootstrap produced an inconsistent performance change (for example, if a performance decrease was observed during photostimulation, the P value was the fraction of times a performance increase was observed during bootstrap).

Electrophysiology data analysisSilicon probe recording preprocessing

The extracellular recording traces were band-pass filtered (300–6 kHz). Events that exceeded an amplitude threshold (4 s.d. of the background) were subjected to manual spike sorting to extract single units26. The primary dataset consisted of 9,626 single units, 73 mice, 347 sessions. The second dataset consisted of 10,420 neurons, 29 mice, 110 sessions.

Spike width was the trough-to-peak interval in the mean spike waveform. Units with spike width < 0.35 ms were defined as fast-spiking neurons (1,045 of 20,046) and units with spike widths > 0.45 ms as putative pyramidal neurons (18,266 of 20,046). Units with intermediate values (0.35–0.45 ms, 735 of 20,046) were excluded from analyses. This classification was previously verified by optogenetic tagging of GABAergic neurons26,44. Unless stated otherwise, we concentrated our analyses on the putative pyramidal neurons.

t-SNE and clustering analysis of individual neuron response profiles

We computed each neuron’s average PSTHs for ‘lick left’ and ‘lick right’ trials (correct trials only) and concatenated the PSTHs. Each PSTH was baseline-subtracted and magnitude-normalized by dividing by the norm of the PSTH. We excluded neurons that did not exhibit consistent PSTHs. Specifically, we split each neuron’s trial data in half and computed PSTHs twice using the split data. We then computed Pearson’s correlation between the PSTHs. Neurons with correlation coefficient less than 0.5 were excluded (3,970 of 20,046). In whole-cell recordings, some neurons did not produce enough spikes to calculate PSTH. However, we found that the PSTH calculated from trial-averaged membrane potential (Vm) closely matched the spiking activity PSTH27. We therefore used PSTHs calculated from Vm for the whole-cell data. Vm was downsampled in time to match the PSTHs from spiking activity.

The input data for the t-SNE were an n × t matrix, where the rows contain the PSTHs of individual neurons. We tested a range of parameters for t-SNE, including perplexity (30 to 1,600), distance metrics (correlation, cosine or Euclidean distance) and the number of principal components (20–100). Only the perplexity affected the embedding outcome, but the results were similar for perplexity 30–100. We therefore used perplexity of 50, 50 principal components and cosine distance for the embedding. We computed t-SNE ten times and picked the outcome with the lowest Kullback–Leibler divergence. We performed t-SNE separately on the primary and second datasets.

ePAIRS test for clustering of response profiles

To test if ALM neurons exhibited clusters of prototypical response profiles or a uniform continuum of response profiles, we used the projection angle index of response similarity (PAIRS) test first presented by ref. 9. We used a modified version of the PAIRS test presented by ref. 12 which accounted for the variance structure of the data, that is, ePAIRS test.

The input data were an n × t population response matrix, where the rows contain the PSTHs of individual neurons (‘lick left’ and ‘lick right’ trials concatenated). The PSTHs were baseline-subtracted and magnitude-normalized. We used PCA to reduce the dimensionality from n to 26, capturing 98% of the activity variance over time. We then examined the loadings matrix (n × 26), which represented the weights of individual neurons for the 26 principal components. Each neuron’s response profile over time was thus represented by a 26-element vector. For each neuron, we computed the average vector angle between the neuron and its k nearest neighbors. Across the population of neurons, we obtained a distribution of average angles. The median of this distribution should be small if neurons exhibited similar response profiles.

For comparison, we generated null distributions that exhibited no clustering. We drew n samples from a 26-dimensional multivariate Gaussian distribution using the MATLAB function mvnrnd(). As in the ePAIRS test presented by ref. 12, we drew samples from a multivariate distribution with zero mean and the variance along each dimension was matched to the neural data. We then computed the average vector angles of nearest neighbors for this simulated dataset. This yielded null distributions for neuronal populations with a uniform continuum of PSTH shapes. To statistically compare the null distributions and the empirical distribution, we simulated the null distribution 10,000 times and calculated the fraction of times the median of the empirical distribution exceeded the median of the null distribution. This corresponded to a P value.

The average vector angle depended on the number of dimensions considered. Here we used 26 principal components, but we also tried as few as eight principal components (the same as in ref. 9) and reached the same statistical outcome. The average vector angle also depended on the parameter k. We tested a range of k (1–10) and arrived at the same statistical outcome. Following ref. 12, we used simulated data to validate the ePAIRS test. We drew samples from either a single multivariate normal distribution (that is, no clustering) or from multiple multivariate normal distributions (that is, multiple clusters). We found that the ePAIRS test appropriately captured only cases where samples originated from multiple distributions.

Clustering of response profiles

Density peak clustering63 was performed in the two-dimensional t-SNE representation. Clustering using the top principal components also produced similar results, but clustering in the t-SNE gave slightly more consistent PSTHs within clusters. Density peak clustering required manual selection of clusters based on local density. We evaluated the robustness of cluster number across a range of population size. Subpopulations were created by subsampling neurons in the dataset and clusters were selected blind to the population size. The number of clusters saturated at ~100 (Extended Data Fig. 1e). To correct for over-clustering, we manually examined the PSTHs of each cluster and combined a small number of clusters (<10%) with very similar PSTHs. The primary dataset yielded 94 clusters. t-SNE and clustering were performed independently on the second dataset, resulting in 86 clusters.

To examine the consistency of response profiles between the primary and second datasets (Fig. 1f,g and Extended Data Fig. 1f,g), we matched clusters across the two datasets. Because the second dataset has a longer delay epoch (1.7 s), we downsampled the PSTHs in the delay epoch using the MATLAB function resample(). For each cluster from the second dataset, we computed Pearson’s correlation for its PSTHs (‘lick left’ and ‘lick right’ trials concatenated) with all the clusters from the primary dataset. The clusters were matched based on the highest correlation coefficient (Extended Data Fig. 1f). In some cases, the cluster with the highest correlation coefficient had already been matched to another cluster. The matched clusters were then defined as the next best match based on Pearson’s correlation coefficient and visual inspection of the PSTHs. We did not exhaustively match all clusters of the two datasets. Rather, we focused on a subset of the clusters from the second dataset (48 of 86) that could be easily matched to a cluster from the primary dataset.

To visualize the response profiles of all clusters (Fig. 1f and Extended Data Fig. 1h), the clusters were first sorted into three groups: lick right preferring, lick left preferring and nonselective. Within each group, clusters were further sorted by activity onset time. For the nonselective clusters, clusters were further subdivided into excitatory and suppressive responses before sorting by activity onset time.

Cluster reproducibility

To evaluate the robustness of clusters from density peak clustering, we also performed Louvain–Jaccard clustering on the primary dataset. We calculated the matrix for co-clustering of every pair of neurons in each method (Extended Data Fig. 1c). We sorted the neurons based on co-clustering in density peak clustering. The block structure along the diagonal of the cell–cell co-clustering matrix was preserved in Louvain–Jaccard clustering, which indicates that if two cells belonged to the same cluster in density peak clustering, then their co-clustering probability was high for Louvain–Jaccard clustering.

To define reproducible clusters, for each cluster in density peak clustering, we found its matching cluster in Louvain–Jaccard clustering. A matching cluster must have >50% of its units also present in the original cluster. By this criterion, 70 of 94 clusters from density peak clustering could be matched to a cluster in Louvain–Jaccard clustering. For the matched clusters, 59 of 70 clusters defined by density peak clustering had >50% of their units captured by their matching clusters defined by Louvain–Jaccard clustering. We considered these clusters to be reproducible. The irreproducible clusters tend to be small in size (Extended Data Fig. 1d), representing 25.8% of the neurons in the dataset.

Noise correlation

We calculated noise correlation between all simultaneously recorded neuron pairs with three constraints (Fig. 1h and Extended Data Fig. 1i). First, neuron pairs must be >100 µm apart. This avoided contamination of spikes from spike sorting. Second, each neuron must be recorded for >10 trials for each trial type. Third, each neuron was only used once in neuron pairs. This avoided using the same neuron in multiple neuron pairs, making neuron pairs nonindependent from each other in statistical tests. In total, we obtained 1,060 within-cluster pairs out of 2,658 possible neuron pairs, and 1,598 across-cluster pairs instead of 107,136 possible pairs. Noise correlation was calculated separately for ‘lick right’ and ‘lick left’ trials and separately during baseline, sample, delay and response epochs. The baseline epoch was 500 ms before the sample epoch. Only correct trials were used. For each trial, we calculated spike counts within 100-ms windows. For each time window, we subtracted the mean spike count calculated from all trials of the same trial type. Noise correlation was Pearson’s correlation of the mean-subtracted spike counts across trials and time windows.

Decoding

Decoding was performed independently at each time point (200-ms windows in 50-ms steps). Decoding was performed using a linear support vector machine on pseudo-populations that combined neurons from different recordings. We concatenated the single-trial spike counts of individual neurons to generate population response vectors. Because neurons were not recorded simultaneously, trials from different neurons were randomly matched. This approach ignored any trial-to-trial correlation between neurons. Response vectors for testing were built using trials that were not used for training. Decoding was repeated 20 times using population responses generated from different combinations of single neuron trial data. Standard errors of mean performance were calculated as the standard deviations of performances across different runs. We used pseudo-populations because most recording sessions did not yield many neurons simultaneously recorded for enough error trials.

To decode trial type, we matched the number of correct and error ‘lick left’ and ‘lick right’ trials. Trials could be classified by stimulus, choice or outcome while holding the other two variables constant29,31,32. To decode trial epochs, correct ‘lick left’ and ‘lick right’ trials were combined. Population response vectors were subjected to four-way classification (baseline, sample, delay or response epochs). To decode reaction time, correct ‘lick left’ and ‘lick right’ trials were combined and sorted by reaction time (interval between ‘go’ cue onset and the first lick). The top and bottom one-third of the sorted trials were used for binary classification. To decode ignore trials, all ‘lick left’ and ‘lick right’ trials were combined and classified by whether mice generated a licking response within a 1.5-s time window following the ‘go’ cue.

Population selectivity vectors

Across n neurons, we concatenated the trial-averaged responses within specific time windows into n × 1 response vectors that described the population response for each trial type: correct ‘lick right’ trials, \(\overline }}}}\); correct ‘lick left’ trials, \(\overline }}}}\); error ‘lick right’ trials, \(\overline }}}}\); error ‘lick left’ trials, \(\overline }}}}\). The population selectivity vectors were calculated as:

$$\begin}\,}\,} = \frac}}}} - \overline }}}} } \right) + \left( }}}} - \overline }}}} } \right)}}\\ }\,}\,} = \frac}}}} - \overline }}}} } \right) + \left( }}}} - \overline }}}} } \right)}}\\ }\,}\,} = \frac}}}} - \overline }}}} } \right) + \left( }}}} - \overline }}}} } \right)}}\end$$

Stimulus selectivity represents the response in posterior object location trials (CR, ER) relative to anterior trials (CL, EL). Choice selectivity represents the response when mice licked right (CR, EL) relative to licking left (CL, ER). Outcome selectivity represents the response in correct trials (CR, CL) relative to error trials (ER, EL).

Selectivity vectors were calculated in analysis windows of 100 ms at different time points (1-ms steps). We quantified the similarity of selectivity vectors across time using Pearson’s correlation (Fig. 2b, off diagonal). Within each analysis window, we calculated Pearson’s correlation between the selectivity vectors calculated from split-half trials (Fig. 2b, diagonal).

Activity modes

Across n neurons, we defined a set of orthogonal directions in activity space (Mode, n × 1 vectors) that captured components of population activity (Fig. 2). We defined the activity modes using a portion of the trials. Separate trials were used for activity projections. At each time point, we calculated the trial-averaged population response vectors (r, n × 1) for specific trial types. Activity projections were calculated as ModeTr. To obtain standard errors, we bootstrapped the neurons in the dataset. Standard error was the standard deviation of the activity projections calculated on the resampled datasets.

We calculated stimulus, choice and outcome modes from selectivity vectors (see Population selectivity vectors). Stimulus selectivity vectors were similar during the sample and early delay epochs (Fig. 2b). We averaged the stimulus selectivity vectors in the sample epoch to obtain the stimulus mode, \(}}}}}_}} = \frac}}}} - \overline }}}} } \right) + \left( }}}} - \overline }}}} } \right)}}\). Choice selectivity vectors developed during the late sample epoch and were stable during the delay epoch (Fig. 2b). We averaged the choice selectivity vectors in the delay epoch to obtain the choice mode, \(}}}}}_}} = \frac}}}} - \overline }}}} } \right) + \left( }}}} - \overline }}}} } \right)}}\). We defined the action mode based on choice selectivity during movement initiation (0.1 s < t < 0.3 s after the ‘go’ cue), \(}}}}}_}} = \overline }}}} - \overline }}}}\). We defined the outcome mode by averaging the outcome selectivity vectors during the response epoch (0 s < t < 1.3 s after the ‘go’ cue), \(}}}}}_}} = \frac}}}} - \overline }}}} } \right) + \left( }}}} - \overline }}}} } \right)}}\).

We additionally defined two non-trial-type-selective activity modes previously shown to play roles in decision-making and motor response30,35. We defined a ramping mode as \(}}}}}_}} = }}}_ - }}}_\), where \(}}}_\) represents the population response vector 500 ms before the sample epoch and \(}}}_\) represents the population response vector during the last 500 ms of the delay epoch. We defined a go mode that captured the phasic activity after the ‘go’ cue35, \(}}}}}_}} = }}}_ - }}}_\), where \(}}}_\) and \(}}}_\) represent the population response vectors 100 ms before and after the ‘go’ cue. The ramping and go modes were calculated using the combined responses from correct ‘lick left’ and ‘lick right’ trials.

Finally, we calculated an activity mode that captured most of the remaining activity variance. We calculated eigenvectors of the population response using singular value decomposition (SVD). The data for the SVD were an n × t population response matrix containing the baseline-subtracted PSTHs of n neurons (‘lick right’ and ‘lick left’ trials concatenated). Consistent with previous analyses of frontal cortex activity7,41, the eigenvector carrying the most variance showed non-trial-type-selective modulation during the response epoch, which we defined as the response mode (Moderesponse).

The seven activity modes were near orthogonal to each other (Extended Data Fig. 2a). For all analyses in the paper, we further rotated the activity modes using the Gram–Schmidt process to be fully orthogonal to each other. Together, the seven activity modes captured 69% of ALM activity variance, 71% of the stimulus selectivity, 92% of the choice selectivity and 93% of the outcome selectivity (Extended Data Fig. 2b). Activity variance was the root mean square (r.m.s.) of the baseline-subtracted activity over the sample, delay and response epochs. The population stimulus and choice selectivity were the r.m.s. values of the stimulus and choice selectivity over the sample and delay epochs. The population outcome selectivity was the r.m.s. of the outcome selectivity during the response epoch.

Our primary analyses were based on neurons combined from different recordings. We also performed the same analysis on simultaneously recorded populations (Extended Data Fig. 2c). We restricted the analysis to sessions with at least ten neurons recorded simultaneously for at least ten trials of each trial type (33 sessions, 10–57 neurons recorded simultaneously, 24 neurons on average). To average activity projections across multiple sessions, we offset the activity projections of each session by subtracting the global mean of activity projections across all trials and time points. This removed session-to-session fluctuations in mean activity. The offsets were computed using the trials that were used to construct the activity modes. Independent trials were used for activity projections.

Activity modes from demixed PCA

Demixed PCA was performed using the dPCA package, https://github.com/machenslab/dPCA (v.1.0.5)7 (Extended Data Fig. 2c,d). The input to the dPCA was an n × s × d × t × k matrix where each entry was the spike rates of individual neurons in individual trials (calculated in 200-ms windows). n corresponds to neurons, s corresponds to trial types instructed by object location (‘stimulus’, anterior versus posterior), d corresponds to lick directions (‘choice’, left versus right), t corresponds to time steps (1 ms) and k corresponds to individual trials. Neurons were combined from different recordings. For conditions with fewer numbers of trials, the empty entries in the response matrix were filled in with nan’s. dPCA was performed using the MATLAB function dpca() in the package with the default parameters. We found that ALM population activity reorganized dramatically after the ‘go’ cue (Fig. 2b). We therefore used the timeSplits option in the dpca() function to split the time periods at the ‘go’ cue for separate marginalization.

ePAIRS test for mixed selectivity

We tested for a notion of mixed selectivity where a shared neuronal population encoded random mixtures of the seven activity modes defined above (Fig. 2). Each activity mode corresponds to a weighted sum of individual neuron activities. Each neuron’s weights for the activity modes constituted a seven-dimensional coding vector. We calculated the angles between each neuron’s coding vector and its nearest neighbors. We then tested if the distribution of the nearest-neighbor angles differed from null distributions expected from random distribution of coding vectors. For null distribution, we drew vectors from a seven-dimensional multivariate Gaussian distribution with variance along each dimension matched to the neural data. The procedures for the ePAIRS test were the same as above (see ePAIRS test for clustering of response profiles), but here for the seven-dimensional vectors (Fig. 3c).

Joint coding of specific activity modes

We characterized individual neurons’ weights for pairs of activity modes as two-dimensional vectors (Fig. 3d). If a pair of activity modes were randomly mixed across a shared neuronal population, the angles of the coding vectors would be uniformly distributed. If neurons coded individual activity modes, the distribution would exhibit peaks at 0° and 90°. Because weights can take on positive or negative values, we limited the angles to 0° and 90° by taking the absolute value of the weights. Finally, we limited the analysis to the top 20% of the neurons ranked by the length of their coding vectors. If neurons were not selective for either activity mode, their coding vectors would appear uniformly distributed even though the vector lengths were very small. Importantly, this selection did not affect the coding vector distribution because neurons coding both activity modes also exhibited large vector length. Using a more inclusive criterion (for example, the top 50% of the neurons) yielded similar results.

Synthetic neuronal population

We generated a synthetic population coding random mixture of activity modes (Fig. 3c,d and Extended Data Fig. 4e–g). We first computed the seven activity modes as described above (see Activity modes). We then found eigenvectors of the population response matrix (n × t) using SVD. The population response matrix contained the baseline-subtracted PSTHs of n neurons (‘lick right’ and ‘lick left’ trials concatenated). We rotated the eigenvectors using the Gram–Schmidt process to be orthogonal to the seven activity modes. Individual neuron PSTHs could be reconstructed from linear combinations of the activity modes and eigenvectors. We constructed synthetic neuron responses using random combinations of the activity modes and eigenvectors (Extended Data Fig. 4e,f). The weights were drawn from a Gaussian distribution with zero mean. The variance of the Gaussian distribution was scaled so each activity mode and eigenvector in the synthetic population carried the same amount of activity variance as in the original population. This procedure thus preserved the activity modes, but randomly redistributed them across the synthetic population. We calculated t-SNE of the synthetic responses (Extended Data Fig. 4g). We also recalculated the activity modes using the synthetic responses and obtained the weights of individual neurons (Extended Data Fig. 4g). We carried out the ePAIRS test as described above on the synthetic population (Fig. 3c and Extended Data Fig. 4b).

Functional populations

We performed k-means clustering on the activity mode weights (n × 7 matrix for a population of n neurons) to divide neurons into functional populations (Fig.

留言 (0)

沒有登入
gif