Dynamic reorganization of the cortico-basal ganglia-thalamo-cortical network during task learning.

Adaptive behavior is coordinated by neuronal networks that are distributed across multiple brain regions such as in the cortico-basal ganglia-thalamo-cortical (CBGTC) network. Here, we ask how cross-regional interactions within such mesoscale circuits reorganize when an animal learns a new task. We apply multi-fiber photometry to chronically record simultaneous activity in 12 or 48 brain regions of mice trained in a tactile discrimination task. With improving task performance, most regions shift their peak activity from the time of reward-related action to the reward-predicting stimulus. By estimating cross-regional interactions using transfer entropy, we reveal that functional networks encompassing basal ganglia, thalamus, neocortex, and hippocampus grow and stabilize upon learning, especially at stimulus presentation time. The internal globus pallidus, ventromedial thalamus, and several regions in the frontal cortex emerge as salient hub regions. Our results highlight the learning-related dynamic reorganization that brain networks undergo when task-appropriate mesoscale network dynamics are established for goal-oriented behavior.

Introduction

Neuronal activity governing behavior is distributed across the brain-wide network of interconnected regions (Kaplan and Zimmer, 2020). For example, goal-directed perceptual tasks require computations engaging various regions—across sensory, motor, memory, and reward circuits—to generate complex apt behaviors (Salkoff et al., 2020; Steinmetz et al., 2019; Chen et al., 2015). Recent advances in electrophysiology (Steinmetz et al., 2018) and optical imaging (Cardin et al., 2020; Najafi and Churchland, 2018) enable direct measurement of large-scale, distributed brain activity in mice performing various tasks, including discrimination of visual (Musall et al., 2019; Steinmetz et al., 2019), tactile (Gilad et al., 2018; Gallero-Salas et al., 2021), auditory (Gallero-Salas et al., 2021; Musall et al., 2019), and olfactory (Allen et al., 2019) stimuli as well as memory-guided tasks (Pinto et al., 2019). These studies showed that in expert animals, task-related behavioral variables are represented in multi-regional activity patterns, indicating distributed and parallel sensorimotor processing. Such coordinated network dynamics, appropriate for solving a specific task, likely emerges during learning and involves various computational aspects (Makino et al., 2016). For example, neural ensembles in sensory and higher-order cortical areas learn to identify and discriminate relevant stimulus features in specific contexts (Li et al., 2008; Wiest et al., 2010; Chen et al., 2013; Poort et al., 2015; Kim et al., 2020) and to associate these with particular outcomes such as reward or punishment; thalamic regions process and channel parallel streams of sensory (Bennett et al., 2019) and motor (Roth et al., 2016; Sauerbrei et al., 2020) information to the cortex; basal ganglia nuclei receive reward prediction error signals (Schultz, 2006) and transform contextual information into goal-directed behaviors (Gerraty et al., 2018; Gremel and Costa, 2013; Samejima et al., 2005); and finally, motor circuits encompassing subcortical and cortical regions undergo refinement to implement temporally precise execution of apt actions (Huber et al., 2012; Peters et al., 2014).

How large-scale brain activity reorganizes during task learning—linking stimulus-related sensory computations to task-relevant motor representations and thereby instantiating behavioral changes—remains largely unknown. In human research, neuroimaging experiments can reveal changes in functional brain networks upon task learning (Finc et al., 2020; Mohr et al., 2016). However, they remain limited in spatiotemporal resolution. In animal research, new opportunities have emerged for previously challenging longitudinal studies of large-scale brain dynamics. One promising approach is repeated imaging of neuronal activity across the entire training schedule, from naive to expert. Such chronic imaging has been applied with cellular resolution in individual or multiple regions (Komiyama et al., 2010; Poort et al., 2015; Chen et al., 2015; Wagner et al., 2019) and with coarse resolution in multiple areas across the dorsal cortex using wide-field calcium imaging (Clancy and Mrsic-Flogel, 2021; Gilad and Helmchen, 2020). The latter experiments are limited, however, to the most accessible cortical areas. While opto-fMRI allows cell-type-specific interrogation of brain-wide circuits in mice and primates (Gerits et al., 2012; Lee et al., 2010), it does not enable task-related measurements with high temporal fidelity.

To extend chronic optical measurements of brain network activity to subcortical regions, we recently introduced a high-density multi-fiber photometry method for simultaneous recordings from many cortical and subcortical regions in the mouse brain (Sych et al., 2019). Here, we apply multi-fiber photometry to track large-scale brain activity changes when mice learn a texture discrimination task. We particularly focus on the cortico-basal ganglia-thalamo-cortical (CBGTC) network (Foster et al., 2021), which is presumed key for forming stimulus-response associations. Beyond analyzing changes in individual regions, we also apply the framework of network connectivity (Bassett and Sporns, 2017) and information theoretic metrics (Vicente et al., 2011) to reveal cross-regional interactions and their stabilization within the CBGTC network during task learning.

Results

Multi-fiber photometry of distributed brain activity during texture discrimination learning

We trained 14 mice in a go/no-go texture discrimination task (Chen et al., 2013). Mice had to lick in response to a “go” texture (rough sandpaper, grit size P100; “Hit” if correct, “Miss” if not licking) and to suppress licking for a “no-go” texture (smooth sandpaper, P1200; “correct rejection” [CR], if correct; “false alarm” [FA], if licking by mistake; Figure 1A; STAR Methods). Using an expert criterion of 70% correct performance, mice learned the task within 1–15 days (average 6.2 ± 4.7 sessions; mean ± SD; one session per day; one mouse learned the task in the first session; the total number of trials ranged from 1,999 to 13,115). We aligned data to the first session after reaching expert criterion and divided them into the naive phase (sessions before first expert session) and the subsequent expert phase (Figure 1B). During task learning, mice refined a set of goal-directed actions. They developed anticipatory whisking of increased envelope amplitude (Figure 1C; 8.1° ± 2.9° and 15.2° ± 8.7° in naive versus expert phase in an early texture presentation window, 3–3.5 s after trial start; p = 0.015; n = 11 mice), likely facilitating early texture discrimination (Chen et al., 2015). In addition, expert mice reported their go decisions faster than did naive animals (licking onset after texture stop 1.5 ± 0.5 s versus 2.1 ± 0.4 s; p = 0.01, Mann-Whitney U test; n = 14 mice). In general, behavior became more stereotypical throughout task learning (Figure S1).

Figure 1

Learning-related changes in behavior and multiregional brain activity

(A) Go/no-go texture discrimination task.

(B) Left: Aligned task performance across 14 mice (means ± SEMs).

(C) Whisker-to-texture touch probability, whisking envelope amplitude, and licking rate for naive (gray, dashed) and expert (blue, solid) mice (n = 11; means ± SEMs; black bars indicate p < 0.05; frame-by-frame Mann-Whitney U test).

(D) Top: Schematic of the CBGTC network with reference to the direct pathway. Bottom: Hypothesized learning-related recruitment of cross-regional interactions in the CBGTC network. S1BF, primary somatosensory barrel field; M1, primary motor cortex; CPu, dorsolateral striatum; iGP, globus pallidus internal capsule; VM, ventral medial thalamus; VL, ventral lateral thalamus.

(E) Schematic of implanted 12-fiber array (left) and of multi-fiber front piece, connector, and fluorescence recording setup (right; 473-nm and 561-nm laser light to excite GCaMP6m and R-CaMP1.07, respectively; complementary metal oxide semiconductor (CMOS) camera sensor to detect all relay fiber signals at 20 Hz frame rate.

(F) Left: Fiber tip positions pooled across all mice (coronal view, anterior-posterior bregma, 1.34 mm). Right: Post hoc histology examples of GCaMP6m expression in CPu, iGP, VM, and RT (reticular nucleus of thalamus). Fiber tip positions indicated by dashed lines.

(G) Example GCaMP6m ΔF/F traces for the 12 brain regions for naive and expert Hit trials. In addition to the regions shown in (D), we targeted LD (lateral dorsal thalamus), RT, and hippocampal subfields: CA1Py (dorsal pyramidal layer), dorsal CA1Mol (dorsal molecular layer), and DG (dentate gyrus).

(H) Pearson correlation coefficients between ΔF/F calcium traces and whisker-to-texture touch (left, n = 11 mice), whisking envelope (center, n = 11), and licking rate (right, n = 14) for Hit trials. Marker size indicates magnitude of correlation coefficient, blue color highlights regions with significant changes from naive to expert phase (p < 0.05, Mann-Whitney U test).

To reveal how cross-regional activity in the CBGTC loop adapts during learning (Figure 1D), we applied multi-fiber photometry with chronically implanted multi-fiber arrays (Sych et al., 2019). In 10 mice, we targeted regions of the CBGTC network plus hippocampal regions with 12-fiber implants; in 4 additional mice, we implanted 48 fibers targeting the same set of regions and many additional regions in more posterior and anterior parts of the brain (Figures 1E and 1F; STAR Methods; for a list of regions and acronyms, see Table S1). Fiber implants were combined with the viral expression of a genetically encoded calcium indicator (GCaMP6m or R-CaMP1.07), enabling simultaneous recording of bulk calcium signals (expressed as ΔF/F) in all targeted brain regions across many trials from naive to expert phase (Figure 1G). For the 12 regions common to all mice, we analyzed the correlation of calcium signals with behavioral variables and the changes upon learning. Nearly all regional ΔF/F signals showed increased correlation with the texture touch in expert mice, most prominently in cortical areas S1BF (primary somatosensory barrel field) and M1 (primary motor cortex), in the ventromedial (VM) nucleus of thalamus, and in basal ganglia regions CPu (dorsolateral striatum) and iGP (globus pallidus internal capsule) (Figure 1H; p ≤ 0.015 for S1BF, M1, VM, and CPu; p = 0.03 for iGP). Correlation with whisking also tended to increase in most regions, although not reaching significance. Correlation to licking decreased in basal ganglia CPu and iGP, cortical S1BF, thalamic VM, LD (lateral dorsal thalamus), RT (reticular nucleus of thalamus), and hippocampal DG (dentate gyrus), CA1Py (dorsal pyramidal layer) regions (p ≤ 0.01 for CPu, LD, DG, and CA1Py; p ≤ 0.05 for iGP, VM, RT, and S1BF). Overall, the representation of the reward-predicting texture stimulus became stronger during learning, whereas the representation of licking action and reward collection decreased in parallel in several brain regions (Figure 1H). These results indicate that the mesoscale network spanning cortex, basal ganglia, thalamus, and hippocampus undergoes a major functional reorganization during learning.

Widespread increase of activity and discrimination power for the reward-predicting stimulus

We next asked how enhanced behavioral performance is reflected in trial-related activity of different brain regions and their ability to discriminate Hit versus CR trials. In naive mice, ΔF/F signals did not significantly differ for Hit and CR trials at early trial times, including the initial whisker-texture touch period, but diverged later when the animal either initiated movements or remained quiet and refrained from licking (Figure 2A). The late action-reward-related activity was particularly prominent in Hit trials. In expert mice, the major calcium signal peak displayed a shift to the earlier texture presentation period in most brain regions (Figure 2B). In general, two prominent activity peaks were discernible: an early peak associated with the texture touch and a late one related to reward collection (in addition, several regions showed activation upon the initial auditory trial-start cue). We quantified changes in regional activity upon learning from the mean ΔF/F signals in an early “stimulus window” and a late “action-reward window”. ΔF/F signals during the stimulus presentation significantly increased in 7 of the 12 regions, whereas action-reward-related signals tended to decrease (Figure 2C). Thus, upon learning the discrimination task, a widespread shift of neural activity occurs from the action-reward period to the time when the reward-predicting stimulus is presented.

Figure 2

ΔF/F responses and discrimination power for the reward-predicting stimulus increase during learning

(A) Example ΔF/F signals (means ± SEMs) in S1BF, RT, and VL from 1 mouse for Hit and CR trials (blue/red, n = 324/326) for naive versus expert sessions (dashed/solid lines). Brackets indicate stimulus window (3–3.5 s) and action-reward window (6–6.5 s).

(B) Example of trial-averaged ΔF/F signals in iGP, CPu, and CA1Py in naive and expert sessions for Hit (left) and CR (right) trials.

(C) Hit-trial-averaged ΔF/F signals in stimulus and action-reward window (top versus bottom), pooled across naive (gray) and expert (blue) sessions (each data point represents a session in n = 13 mice). False discovery rate (FDR) correction was performed with Bonferroni adjustment.

(D) Example of the Hit/CR classification accuracy across trial time for naive (dashed) and expert (blue) trials (same mouse/sessions as in A and B). Classification with shuffled trial labels was at 50% level (solid black lines).

(E) Hit/CR classification accuracy (n = 14 mice). Mann-Whitney U test for naive versus expert sessions, ±10 sessions before (gray) and after (blue) the first expert session. Vertical axis ranges from 50% to 60% classification accuracy. Color bar shows SEM across mice.

(F and G) Same as (D) and (E), but for Hit/CR classification accuracy estimated simultaneously from all brain regions.

(H) Schematic of regions increasing peak calcium signals during texture presentation.

(I) Schematic of regions increasing task-related classification accuracy during learning.

(C–E) ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001; Mann-Whitney U test, FDR corrected.

To assess how well and when regional activity can discriminate trial types during trial time, we next quantified Hit/CR classification accuracy using linear discriminant analysis (LDA). We trained the classifier on each time bin and for each brain region separately (random 80% training data; 20% test data; 5-fold cross-validated). In the naive phase, classification accuracy reached 60%–80% only late in the trials, specifically in the action-reward window (Figure 2D). In expert mice, all of the regions showed an earlier increase of classification accuracy—in several regions as early as the texture touch. To evaluate stimulus-related classification accuracy across learning while avoiding confounding effects of action-related licking, we trained the LDA classifier on the mean ΔF/F signal in the stimulus window and aligned all of the sessions to the first expert session. Task-related Hit/CR classification accuracy increased from naive to expert phase in nearly all of the regions, reaching significance in S1BF, iGP, CPu, VL (ventral lateral thalamus), RT, LD, and DG (Figure 2E; comparing 10 sessions before and after first expert session). Hit/CR classification reached higher accuracy (∼70%) for late expert sessions when we applied LDA across all regions simultaneously (Figures 2F and 2G). For further verification that Hit/CR classification in experts reflects discrimination of stimuli rather than actions, we also classified CR versus FA trials (same stimulus, different actions). Action-related CR/FA classification accuracy transiently increased during the learning period (±1 sessions around time zero) in several regions, but gradually decreased again in the late expert phase, indicating that motor actions alone cannot fully explain Hit/CR discrimination power (Figure S2). We conclude that the widespread shift of activity to the stimulus window (Figure 2H) is paralleled by enhanced go/no-go (Hit/CR) classification in partly overlapping sets of regions (Figure 2I). The improvement of Hit/CR classification when all of the regions are used for LDA prompted us to further study learning-related adaptations on the network level by quantifying cross-regional interactions.

Cross-regional interactions increase during learning

To evaluate cross-regional interactions, we analyzed effective connectivity at the two relevant timescales of trials and of learning. We estimated effective connectivity using transfer entropy (TE) (Novelli et al., 2019; Schreiber, 2000; Vicente et al., 2011). TE is a directed measure of effective connectivity based on information-theoretic measures of regional activity distributions across trials and time steps (Battaglia et al., 2012). For every pair of regions, an effective connection is assigned if a signal of a source region in a past time bin t-Δt is predictive of the signal of a target region in the present time bin t, but only if the latter cannot be explained by its own past. We used mainly a multivariate TE approach (Wollstadt et al., 2019), which applies an optimal information transfer criterion to prune weaker connections that can be fully explained by stronger ones. From the experimental point of view, we used multi-regional viral labeling, which may result in both neuronal cell bodies and long-range projecting axons (neuropil) to express GCaMP6m. Therefore, part of the signal from a source brain region could contaminate the signal in a target brain region via axonal cross-talk. We evaluated the sensitivity of our TE analysis to this possibility in simulations (Figure S3). Whereas inter-areal axonal contamination of calcium signals (cross-talk) increases zero-lag correlations between brain regions, the TE approach is fairly resistant to this effect because it statistically tests the dependence of present signals on past dynamics. Consequently, moderate levels of axonal cross-talk minimally affect estimates of cross-regional connectivity obtained with TE (Figure S3).

To analyze effective connectivity changes across learning, we estimated TE in 200-ms time bins in the stimulus window (3.0–3.6 s; 3 time bins; Figure 3B). Estimates were made separately for each session from all within-session trials. A connection between two regions was labeled as existent if TE was significant in at least one of the three stimulus time bins (p < 0.01; compared to shuffled data; STAR Methods). For each session, the total number of detected network connections was calculated (the maximum number of possible connections is k · (k − 1); i.e., 132 for k = 12 regions). After aligning sessions to the first expert session and pooling across animals, we found that the total number of stimulus-related connections increases from naive to expert phase, for both Hit and CR trials and using either bivariate or multivariate analysis (Figures 3C and 3D; we mainly use multivariate analysis as the more conservative approach; STAR Methods). Moreover, we defined connection strength as the mean frequency of significant TE bins, under false discovery rate (FDR) correction, in a given period (e.g., for the 3 time bins in the stimulus window, frequencies of 0.0, 0.33, 0.67, and 1.0 are possible for an individual session). The distribution of connection strength (averaged across sessions) shifted toward higher values from naive to expert phase for both Hit and CR trials (p < 0.001; two-sample Kolmogorov-Smirnov test), with a subset of network connections clearly gaining strengths (Figures 3E and 3F for Hit and CR trials, respectively). In expert mice, the 8 strongest connections encompassed cortical, thalamic, and hippocampal regions, with substantial overlap between Hit and CR trials (Figure 3E and 3F). Cross-regional interactions were not correlated with the magnitude of stimulus-evoked activity (Figure S4), although both increased in parallel during task learning. We conclude that cross-regional mesoscale interactions are enhanced during task learning at stimulus presentation time, consistent with the observed widespread increase of ΔF/F activity at this moment during trials.

Figure 3

Learning is associated with enhanced network connectivity

(A) Illustration of transfer entropy analysis. Top: Signal B depends on past dynamics of A. Bottom: Projection of the 2 time series (blue) onto the subspace. A(t) × B(t + Δt) shows a clear correlation (green), whereas projection onto B(t) × B(t + Δt) shows no correlation (red).

(B) Example of TE estimates across trial time for all connections from VM to other target regions for expert Hit trials (same mouse as Figures 2A and 2B). Black bars indicate frames with significant TE; FDR corrected.

(C) Total number of connections in stimulus window (3–3.6 s) for Hit (blue) and CR (red) trials across sessions (means ± SEMs, n = 13 mice) with bivariate (left) and multivariate (right) analysis.

(D) Total number of connections for naive and expert phase, estimated with bivariate (left) and multivariate (right) analysis (means ± SEMs, n = 13 mice; gray lines, individual mice; ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001; Mann-Whitney U test).

(E) Top left: Histogram of connection strength, averaged across mice, during Hit trials for naive (gray) and expert (blue) phase. Top right: The 8 strongest connections in expert sessions shown as network of the respective regions. Bottom: Matrices of connection strength averaged for naive (left) and expert (right) phase (averaged over n = 13 mice).

(F) Same as (E) for CR trials.

Stabilization of cross-regional interactions during learning

We next asked whether such effective connectivity changes also occur during other specific trial periods and how they compare for different trial types. Expert mice displayed higher effective connectivity primarily during the texture presentation period when mice gathered relevant information to discriminate between textures (Figure 4A). During later trial times, the total number of connections decreased, converging to the naive level. This observation, together with the increased numbers of connections in both Hit and CR expert trials, renders it unlikely that effective connectivity simply relates to animal movements. Consistent with this notion, subdividing sessions according to high and low values for specific behavioral variables revealed no dependence of the total number of effective connections on whisking amplitude (albeit on touch probability in experts) and even a higher number for lower licking rates in experts (Figure S5). We used variance partitioning to estimate whether behavioral variables or task performance better explain connectivity emerging during learning. Task performance uniquely explained more variance in connectivity as compared to behavioral parameters (Figure S6).

Figure 4

Stabilization of trial-related network dynamics during learning

(A) Total number of connections throughout trial time for Hit (blue) and CR (red) trials (means ± SEMs; n = 13 mice; solid/dashed lines expert/naive phase; SEMs, gray-shaded areas).

(B) Bar plots of naive and expert mean clustering coefficient (± SEMs) for those regions with significant increases (Hit and CR expert conditions, blue and red; naive conditions, gray; FDR corrected; gray lines indicate individual mice).

(C) Number of connections shared across sessions for naive and expert conditions (means ± SEMs; Mann-Whitney U test; gray lines indicate individual mice).

(D) Left: Scatterplot shows the likelihood-to-be-random (log p value) of the observed number of shared connections in the stimulus window versus task performance (top, Hit trials; bottom, CR trials; solid lines, moving averages; each dot = 1 session; horizontal dashed line, p = 0.01; vertical line, 70% performance threshold); >, above upper bound, <, below lower bound, of random distribution. Right: Violin plot of log p distributions for naive and expert phase (Mann-Whitney U test).

(E) Same as (D), but for action-reward in Hit trials (top) and reward-omission in CR trials (bottom; outcome window 6–6.6 s).

(F) Strength for all connections across trial time for expert phase (mice randomly split into train [n = 7] and test [n = 6] sets). Across-mice average connection strengths for Hit and CR sessions were sorted from early- to late-active according to the peak time. Sorting of the train set was applied to the test mice. Average connection strength for connections 1–40 (group 1), 41–80 (group 2), and 81–120 (group 3) are shown to the right for the train (black) and test (magenta) data.

(G) Rank distribution of the difference between Hit and CR connection strength for all connections (positive Hit-dominated ranks, blue; negative CR-dominated ranks, red; null distribution of ranks, gray; Mann-Whitney U test, n = 13 mice).

(H) Most frequent connections for Hit (blue) and CR (red) trials in expert phase (Wilcoxon signed-rank test p < 0.05, not FDR corrected; n = 13 mice). ∗p < 0.05, ∗∗p < 0.01, ∗∗∗, p < 0.001, ∗∗∗∗p < 0.0001.

To assess whether enhanced connectivity may imply that some regions selectively increase their outgoing connectivity, we evaluated naive-to-expert changes in the clustering coefficient (Barrat et al., 2004; Rubinov and Sporns, 2010) for the stimulus window (STAR Methods). The clustering coefficient consistently increased across mice for subcortical regions iGP, VM, and RT in Hit trials, and for CA1Py and RT in CR trials (Figure 4B). This result indicates that the brain network adapts in a reliable way and stabilizes in an appropriate reorganized state. Consistently, the number of connections shared across consecutive sessions increased from naive to expert phase (Figure 4C).

To further assess network stabilization, we compared session-to-session changes in estimated connectivity with random network dynamics. Specifically, we calculated the likelihood (log p value) of obtaining the same number of session-to-session shared connections, assuming that the exact connections were randomly permuted for each session (null model). As task performance improved, this likelihood-to-be-random markedly decreased in the stimulus window, whereas in the naive phase, randomized connectivity could explain the observed number of shared connections (Figure 4D; pooled across sessions; p < 0.01). Session-to-session network dynamics during the action-reward window for Hit trials and reward-omission for CR trials also showed a stabilization trend (Figure 4E; p < 0.01). For both time windows, these results were comparable for Hit and CR trials. We conclude that mesoscale effective connectivity stabilizes throughout learning, especially during the early texture presentation period, with a subset of regions establishing strong network clustering motifs.

We also analyzed trial-related dynamics of connection strength for all possible connections in the 12-region network. For a subset of expert mice (“train,” n = 7) we sorted the connections according to peak time of their connection strength, revealing a sequential activation of connections during trial time (Figure 4F; STAR Methods). To test for consistency across mice, we applied the same sorting to both naive and expert phase of the remaining “test” subset of mice (n = 6). To visualize the sequential recruitment of connections, we divided connections into three groups of equal size. Group 1 showed high activity early in trial time, especially between the auditory cue and texture arrival; group 2 displayed the highest activity during texture touch; and connections in group 3 were most active during licking and reward time. The group assignment from the train mice was partially conserved in the test mice (Figure 4F). In experts, approximately one-third of connections engaged strongly after the auditory cue and peaked during the texture presentation period, displaying a temporal profile similar to the ΔF/F calcium signals (Figure 4F). To test for differences in effective connectivity between Hit and CR trials, we ranked Hit minus CR strength differences individually for each connection (averaged from 1 to 8 s of trial time) across all mice, summed the ranks over all connections, and compared the resulting test statistic to a null model (STAR Methods). The resulting test statistic was significantly larger than zero, demonstrating that overall connection strength was higher for Hit compared to CR trials (Figure 4G). To identify differences in Hit and CR networks, we tested separately for each connection—and in each time frame—whether connection strength differed between Hit and CR trials. Only a few connections, all between thalamic regions VM/VL and S1BF or iGP, consistently showed Hit/CR differences across mice at specific trial moments (Figure 4H). We conclude that connections are activated in sequential order during trials, with the rough order according to the most salient events being preserved across mice. In addition, networks engaged in Hit and CR trials share a substantial amount of connections, with only a few trial-type specific connections.

Internal globus pallidus coordinates effective connectivity in thalamus and cortex

To reveal which regions strongly contribute to functional network formation, we quantified how Hit/CR differences of the clustering coefficient emerge during learning. In general, we found higher clustering coefficients in Hit trials. In naive mice, Hit/CR differences mainly occurred around the time of reward, and they were significant for several regions in the action-reward window (CPu, iGP, VL, S1BF, M1, and CA1Mol [dorsal molecular layer]; p < 0.05; Figures 5A and 5B). In experts, such differences occurred at earlier times during texture presentation, with a significant increase in the stimulus window only for iGP (p = 0.002). The basal ganglia, including its output region iGP, thus may learn to integrate the reward-predicting signal and coordinate appropriate motor commands in the downstream thalamic targets.

Figure 5

Dynamic changes of regional clustering during learning

(A) Hit-CR differences (Δ) of clustering coefficient (ΔCC) versus trial time for all 12 regions (n = 13 mice). Naive phase, left; expert phase, right. Dashed boxes, stimulus, and action-reward window.

(B) Clustering coefficient during Hit and CR trials for regions showing significant differences in either stimulus or action-reward window (∗p < 0.05, ∗∗p < 0.01; Mann-Whitney U test; means ± SEMs; n = 13 mice; gray lines are individual mice). No significant differences were found in stimulus window for naive mice and in action-reward window for expert mice.

(C) Regions targeted by iGP during cue (1–1.5 s), stimulus (3–3.5 s), and action-reward (6–6.5 s) windows in Hit (top row) and CR (bottom row) trials. Line thickness indicates overall number of observed targets pooled across all mice (sorted counterclockwise from highest to low).

(D) Probability to be targeted by iGP (mean + 3 SD; Hit, blue; CR, red) comparing action-reward window for naive and stimulus window for expert phase. Open bars show mean + 3 SD for shuffled connections. Only regions above shuffle mean + 3 SD are shown.

(E) Left: Schematic of optogenetic excitation of iGP neurons. Center: Proportion of Miss and FA trials without and with laser illumination. Right: Overall task performance without and with laser illumination (9 independent iGP sites in 7 mice pooled; means ± SEMs; Wilcoxon signed-rank test).

(F) Violin plots for the distribution of mean ΔF/F R-CaMP1.07 signal during a late-reward period (6.5–7 s) in VL (n = 4 mice; 336 Hit, 68 Miss/LaserOff, and 73 Miss/LaserOn trials), LD (n = 4; 640 Hit, 144 Miss/LaserOff and 162 Miss/LaserOn trials), and S1BF (n = 6, 8 sites pooled; 339 Hit, 95 Miss/LaserOff, and 103 Miss/LaserOn trials). These regions are downstream of iGP, consistent with increased targeting of VL, LD, and S1BF in (C) and an increased clustering coefficient. Red line, mean; black line, median. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001; Mann-Whitney U test across trials.

In expert mice, we further analyzed which regions contribute to the formation of triangular network motifs as captured by the iGP clustering coefficient. Most frequently, iGP connected to thalamic regions VL, VM, and LD and to cortical S1BF and M1. The connectivity pattern varied across trial time, with VL and VM the most prominent targets before texture presentation and during the action-reward window, whereas S1BF, LD, and VL were targeted during the stimulus window (Figure 5C). Compared to random networks, S1BF and LD received above-chance iGP connections in Hit trials and VL and LD in CR trials, whereas in naive animals, only VM was above chance in the action-reward window (Figure 5D).

To establish a link between network connectivity and mouse behavior, we used optogenetics to activate GABAergic neurons in VGAT-ChR2 EYFP mice (n = 7 mice across 9 independent sites; STAR Methods) (Sych et al., 2019). Guided by the magnitude of the estimated clustering coefficient, we targeted iGP through one of the 120-μm optical fibers and applied laser illumination (473 nm) during the entire texture presentation period (2–6 s of trial time; modulated at 20 Hz with 20-ms on and 30-ms off periods). This manipulation increased the proportion of Miss trials from 7.8% to 25% (Figure 5E; n = 7 mice, p = 0.003), consistent with the notion that iGP suppresses motor actions via thalamo-cortical inhibition and thus lowers task performance (the average number of FA trials remained unchanged). To quantify the effect of activating GABAergic neurons in iGP on the downstream targets VL, LD, and S1BF, we compared their mean ΔF/F signals (measured with R-CaMP1.07) with and without laser illumination directly after texture presentation (during 6.5–7 s of trial time; STAR Methods). We found a significant decrease in both VL and S1BF signals (Figure 5F), consistent with the expected mainly inhibitory effect on thalamus (Parent et al., 2001; Wallace et al., 2017; Kim et al., 2017). These findings highlight the impact of iGP on CBGTC circuit dynamics and the reorganization of brain networks during learning (Schechtman et al., 2016), as well as the strong pallidal influence on downstream thalamic regions. While the effective connectivity to thalamus can be explained by the direct anatomical projections (Oh et al., 2014), connections to S1BF and M1 cortical regions require TE transfer via an intermediate region involving at least two synapses.

Prefrontal regions integrate input from a wide range of brain regions

In 4 mice we implanted 48-fiber arrays (Sych et al., 2019) to investigate learning-related network dynamics on an even broader scale, targeting multiple additional brain regions in frontal and posterior cortex and in the amygdala (Figure 6A; Table S1). Our recordings from this larger brain network (2,256 possible connections) substantiated our main findings regarding recruitment and stabilization of mesoscale networks during learning. As for the 12-region network, the total number of effective connections increased from naive to expert phase in the 48-region network, particularly during the texture presentation period (Figure 6B). Assessment of network stabilization (comparing session-to-session dynamics with random networks as in Figure 4) again revealed that the likelihood of random session-to-session variability decreased with increasing task performance. This decrease was significant in the stimulus window for both Hit and CR trials when pooled across all of the sessions and mice (Figure 6C; p < 0.001). For both Hit and CR trials we observed no significant gain in stability in the outcome window (Figure 6D), possibly due to the higher variability of motor behavior in this trial period. These results highlight the importance of a brain-wide functional network attaining stable and robust dynamics, specifically during the decisive texture touch period when sensory information is pre-processed and relayed to the motor system to generate appropriate motor commands.

Figure 6

Recruitment of prefrontal cortical circuits during learning

(A) Top view of 3 implanted multi-fiber arrays: one 24-fiber array at −1.4 mm bregma, targeting among other medial (m) regions the same set of regions as in the 12-fiber experiments; two additional 12-fiber arrays at +2.4 mm and −4.1 mm bregma targeting additional anterior (a) and posterior (p) regions.

(B) Changes in total number of connections across trial time (means ± SEMs; n = 4 mice; Hit, blue, CR, red; dashed lines and shaded gray areas, naive phase).

(C) Left: Scatterplot shows the likelihood-to-be-random (log p value) of the observed number of shared connections in the stimulus window versus task performance (top, Hit trials; bottom, CR trials; solid lines, moving averages; each dot = 1 session, horizontal dashed line: p = 0.01; vertical line, 70% performance threshold); >, above upper bound, <, below lower bound, of random distribution. Right: Violin plots of log p-distributions for naive and expert phase (Mann-Whitney U test, ∗∗∗p < 0.001).

(D) Same as (C), but for action-reward in Hit trials and reward-omission in CR trials in the outcome window.

(E) Strength of all connections across trial time (4 mice average) for expert phase for Hit and CR trials. Left: Connections sorted from early to late activity peak times. Right: Change of connection strength (ΔCS) sorted in descending order according to mean value in the stimulus window.

(F) Mean connectivity matrices of ΔCS during the stimulus window.

(G) Network of regions for the top 1% connection strength for Hit (blue) and CR (red) trials in expert phase.

We further aimed to identify in this larger network those connections that clearly adapted their strength during learning. Connection strength dynamically changed during trial time (Figure 6E). As in Figure 4F, connections were sequentially activated throughout the trial when sorted according to their peak activation times. We calculated the connection strength change ΔCS in the stimulus window relative to the pre-stimulus period (1–1.5 s). Approximately one-third of connections increased their strength compared to baseline, whereas another third displayed suppression. To identify which connections adapted their strength during the stimulus window, we plotted the mean ΔCS values for the stimulus window as matrices for naive and expert phases (Figure 6F). In the naive phase, target regions in the frontal cortex (e.g., PrL [prelimbic cortex], Cg1 [cingulate cortex, area 1], MO [medial orbitofrontal cortex], AI [agranular insular cortex]) received strong input from the source regions of the CBGTC network (labeled m-subsection in Figure 6A). The subcortical deep mesencephalic nucleus (DpMe) and thalamic VM also received broad input. In the expert phase, Cg1, MO, PrL, DpMe, and VM maintained their strong input connectivity acting as integrative hub regions (Figure 6F), whereas AI instead received a reduced number of inputs. Other regions (e.g., S2) also showed marked changes in their input patterns. These findings demonstrate how stimulus-associated effective connectivity in the CBGTC loop and beyond reorganizes during learning, comprising both strengthening and weakening of specific connections.

To visualize the most salient stimulus-associated connectivity motifs in experts, we plotted the strongest 1% of connections in the stimulus window (Figure 6G). PrL, MO, and thalamic VM appeared as hub regions. Particularly, MO received feedforward inputs from S1BF during Hit trials and from S1BF and M1 during CR trials, corroborating the idea that signals necessary to update the decision strategy should converge to MO (Sul et al., 2010). We observed strong connections between MO and PrL, supporting the hypothesis that these regions interact as decision centers, possibly containing a map of learned stimulus-value associations (Wilson et al., 2014). Taken together, these results demonstrate that, in addition to the engagement of the CBGTC loop, prefrontal cortical circuits are recruited during task learning (Le Merre et al., 2018; Nakajima et al., 2019).

Discussion

Multi-fiber photometry provides opportunities to study distributed brain activity and its relationship to behavior (Kim et al., 2016; Sych et al., 2019). Here, we demonstrated how this method can be chronically applied to reveal functional reorganization of brain networks during learning. We found widely distributed learning-related changes across brain regions. A large fraction of connections in the subnetworks we recorded from with either 12- or 48-fiber arrays were recruited into the task-related mesoscale network. Upon learning, network connectivity increased throughout most of the trial time, starting from the initial cue signaling the texture approach, up to and throughout the decision-relevant texture presentation period. Multiple regions, especially in the frontal cortex, increased their activity during the texture approach period, suggesting a role in attentive and anticipatory circuit components. Consistent with previous studies, sensorimotor processing in expert mice engaged functional circuits in the cortex (Poort et al., 2015; Wiest et al., 2010; Yang et al., 2016), striatum (Yin and Knowlton, 2006), thalamus (Schmitt et al., 2017; Wimmer et al., 2015), and hippocampus (Le Merre et al., 2018; Pereira et al., 2007; Itskov et al., 2011). A salient characteristic observed in many regions was the temporal shift of activity and discrimination power from the unexpected reward to the reward-predicting stimulus upon learning. This shift is reminiscent of the temporal shift of dopaminergic neuron activity in reinforcement learning studies (Cox and Witten, 2019; Watabe-Uchida et al., 2017). In experts, touch-related activity may not only reflect refined sensory representations, particularly for the rewarded stimulus, but may also comprise neural signals involved in task-specific movements (Engelhard et al., 2019; Howe and Dombeck, 2016), initiation or suppression of a learned motor program, or prediction of reward outcome. The temporal shift from reward period to early texture presentation period was particularly apparent in the local network connectivity of iGP, as estimated by the clustering coefficient. We hypothesize that this dynamic shift may be established by the local striatal circuits integrating reward expectation signals with the relevant motor commands. In this scenario, iGP acts as a control region coordinating downstream activity in thalamic, cortical, and habenula circuits (Stephenson-Jones et al., 2016). iGP may act as a “brake” for early licking responses, potentially increasing the sensory integration time and, hence, the confidence of a behavioral choice. Consistent with this notion, optogenetic activation of GABAergic iGP neurons increased the proportion of Miss trials and reduced mesoscale brain activity primarily in anatomically distinct thalamic and cortical regions. Our results encourage one to think of the task learning not only in terms of sensory selection but also appropriate inhibitory control of motor behavior.

By studying functional changes on the mesoscale level, our results provide direct evidence of recurrent interactions between CBGTC network regions and the hippocampus. Multiple previous studies observed task-related brain-wide activity (Allen et al., 2017, 2019; Musall et al., 2019; Pinto et al., 2019; Salkoff et al., 2020; Steinmetz et al., 2019), but it remained unclear how it is established. Here, we show that learning engages signal flow across the CBGTC network, which shapes cortical dynamics throughout the task training. Several mechanisms could underlie the recruitment of cross-regional interactions and their stabilization during learning, including synaptic plasticity (Rioult-Pedotti et al., 2000), modulation of inhibition (Khan et al., 2018), cross-regional synchronization that could facilitate recurrent interactions (Vezoli et al., 2021; Hipp et al., 2011; Antzoulatos and Miller, 2014; Grion et al., 2016), computations of prediction error from midbrain nuclei (Cox and Witten, 2019; Hosp et al., 2011), and common top-down feedback from prefrontal cortical regions (Le Merre et al., 2018; Paneri and Gregoriou, 2017; Ferenczi et al., 2016). Resolving the mechanisms underlying the temporal shift of activity from reward to reward-predicting stimulus will be particularly important, given its fundamental role for credit assignment in reinforcement learning. Chemogenetic or optogenetic manipulations of specific pathways or subsets of pathways could help to verify their causal role in learning-related brain-wide network adaptations. Several further implications can be drawn from the observed network transformations. For example, increased connectivity implies shorter path lengths to connect any pair of regions, facilitating the information flow across the entire network (Bassett and Sporns, 2017). In addition, as networks become more robust, deletion or perturbation of single regions may be compensated by alternative pathways sustaining a large network motif (Albert et al., 2000).

Limitations of the study

We applied two major approaches to analyze learning-related changes in brain activity. First, we focused on stimulus-induced activity for each brain region and its correlation to behavior. Second, we applied effective connectivity measures to estimate cross-regional interactions, an approach previously used mainly for the analysis of fMRI data (Bassett et al., 2011, 2015) and more recently also of mesoscale calcium imaging data (Kuroki et al., 2018). We wish to stress that the terms used here to describe functional relations, such as cross-regional interactions, connectivity, and connection strength should be interpreted in a purely correlative and not causal manner. In the future, both types of analysis will benefit from refined photometry approaches that should overcome several limitations. Depending on labeling approach and brain region, fiber-photometric signals may represent a mixture of cellular signal components. Here, we expressed calcium indicators using rather generic adeno-associated virus (AAV) promoters and serotypes. We therefore cannot exclude that some regional signals comprised multiple neuronal subpopulations with distinct response properties. In addition, multiregional viral injections also cause labeling of long-range axonal projections, which may also contribute to regional signals. In simulations, we found that TE analysis is more tolerant to such axonal contributions as compared to cross-regional correlation (Figure S3), presumably because TE infers information transfer from time-lagged activities both within a region and across regions. Cell-type or pathway-specific expression of calcium indicators, especially using the recently introduced soma-targeted calcium indicators (Chen et al., 2020; Shemesh et al., 2020), should help to overcome these limitations and thus enable the dissection of signal components and a more refined analysis of regional activities and interactions. Moreover, sensitive genetically encoded voltage indicators could be combined with photometry to overcome the limitations of calcium sensors in temporal resolution (Marshall et al., 2016).

In summary, large-scale calcium recordings open new avenues to reveal the dynamic reorganization of cross-regional interactions during learning. Identification of the pathways recruited during learning and of key features of brain network adaptations will help us to better understand the principles of brain dynamics.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Fritjof Helmchen (helmchen@hifo.uzh.ch).

Materials availability

Calcium sensors are readily available on Addgene https://www.addgene.org/viral-service/aav-prep/sensors/.

Experimental model and subject details

All animal experiment procedures were carried out according to the guidelines of the Veterinary Office of Switzerland and approved by the Cantonal Veterinary Office in Zurich.

C57BL/6 mice

A total of 20 male C57BL/6 mice (2–6 months age) were used in this study. Mice were injected with the viral construct and implanted with the high-density multi-fiber arrays. Mice were housed in the animal facility in a room on reversed light cycle. After the surgery, mice were single housed in T 2L (IVC) cages for the recovery and further behavioral training. Handling mice by the experimenter started in the first week post-surgery and included gentle lifting on the experimenter’s palm. After the second week post-surgery, mice were gradually accommodated to the head fixation. Head fixations were done in the experimental room in the open T 2L (IVC) cage after 5 min exploration period. Fixation periods initially were brief, lasting ∼10 s. Over a few days the duration of head fixation was gradually increased from seconds to 5 min. After the third week post-surgery mice were water-scheduled according to the standard experimental procedures and trained in the texture discrimination task. Every day prior to training the mouse’s weight was recorded. Generally, pre-training weight ranged between 25-35 g.

VGAT-ChR2 EYFP transgenic mice

A total of 7 male VGAT-ChR2 EYFP transgenic mice (2–6 months age) were used in this study. We used the same standard experimental procedures for VGAT-ChR2 EYFP transgenic mice as described above for C57BL/6 mice.

Method details

Surgical procedures

GCaMP6m or R-CaMP1.07 expression was induced by stereotaxic injections of AAV2.9-hSyn-GCaMP6m and AAV1-EF1a-R-CaMP1.07, respectively. GCaMP6m was used for 12-fiber array implants (n = 7) and R-CaMP1.07 for 12-fiber arrays (n = 3) and 48-fiber arrays (n = 4). Several expert sessions out of 6 mice with the GCaMP expression overlap with the previous study (Sych et al., 2019). Multi-fiber implantation was performed in a separate surgical session, typically 2–5 days after virus injections (Schlegel et al., 2018). Briefly, we anesthetized mice with 2% isoflurane (mixed in pure oxygen) and maintained body temperature at 37°C with a heating pad. To prevent inflammation and for analgesia, we applied Metacam (s.c., 0.1 μL/g bw). After removal of connective tissue, we polished and dried the skull and applied iBond (Kulzer, Total Etch) to ensure best adhesion of the skull to the connective dental cement. To further stabilize the implant, we formed a thin ring of Charisma (Kulzer, A1) on the skull rim. Both iBond and Charisma were light-cured. Small slit-like craniotomies were made to allow for virus injections and fiber-array implantations. First, ∼120 nL of virus-containing solution were pressure-injected into all areas of interest at a rate of about 20 nL/min. In order to allow for local diffusion and to minimize potential reflux, the glass injection pipette (10–15 μm diameter) was kept in place for 10 min after each injection. Second, the fiber-array front piece(s) was (were) implanted with the help of a stereotaxic manipulator. The 12-fiber array was implanted at 0.2-mm distance from midline, tilted at an angle of 15 degrees. We oriented the fiber array such that the most lateral fiber efficiently targeted CPu (−1.06 mm from bregma) and the most medial fiber targeted hippocampal areas and posterior motor cortex (CA1, DG, M1; −1.46 mm). Prior to fiber implantation we slightly scratched the dura surface. After, we sealed the craniotomy with vaseline, which melts at body temperature and completely covers the craniotomy. For the 48-fiber experiments, a 24-fiber array was implanted similarly and in an equivalent position as in 12-fiber experiments (−1.46 mm from bregma). Two additional 12-fiber arrays were implanted in frontal and posterior regions (+2.4 mm and −4.1 mm from bregma, respectively). We applied dental cement (Tetric EvoFlow A1) on the skull around the implants followed by UV light curing. A light-weight metal head post was additionally cemented to the skull, allowing for head-fixation during the behavioral experiments.

Training procedures for the texture discrimination task

Mice recovered for 2 weeks after fiber-array implantation. We habituated mice to head-fixation through a series of short-duration head-fixations in the experimental holder. A day after water scheduling started (for 5 days a week), a mouse was placed in the behavioral setup and first trained to lick following texture presentation. During this shaping period, we presented mostly the ‘go’ texture (rough sandpaper, grit size P100). The auto-reward function automatically opened the water valve at 6 s of trial time conditioning mice to lick on the water spout. The trial structure during the shaping period was the same as for the expert performance, i.e. both auditory cues and the texture were presented. After mice robustly licked, we also presented the ‘no-go’ texture (smooth sandpaper, P1200) and trained mice to discriminate between the two texture types. Mice had to lick for the ‘go’ texture (‘Hit’ if correct, ‘Miss’ if not licking) and to suppress licking for the ‘no-go’ texture (‘correct rejection’, CR, if correct; ‘false alarm’, FA, if licking wrongly). We permitted early licks on go trials whereas licking on FA no-go trials was punished with a prolonged white noise sound. The lick detector was reachable throughout the entire session. Textures were presented pseudo-randomly with no more than three consecutive presentations of the same texture type. In each trial, a first auditory cue (2-kHz tone, two 100-ms pulses at an interval of 50 ms) was given at 1 s after trial start, signaling that the texture was starting to move towards the whiskers on the right side of the snout. The whisker-to-texture touch typically occurred between seconds 2 and 3 of trial time before the texture reached its end position at 3 s. The texture stayed in place for 2 s and was then retracted, also indicated by a second auditory cue (4-kHz tone; four 50-ms pulses at 25-ms interval). A small water reward was delivered in Hit trials (total reward during a full expert session amounted to 1.5 - 3 mL). Mice reported their decision by licking. Typically, mice picked a reward between 5.8 – 6.5 s as can be seen in the increased licking rate starting from 6 s (Figure 1C).

We trained 14 mice. Mice showed different learning rates and on average learned the task within 6 days from start of discrimination training (6.2 ± 4.7 sessions; mean ± s.d.; 1 session per day; the slowest mouse required 15 sessions; one mouse learned the task in the first session and therefore could not be included in session-based analysis; total number of trials ranged from 1999 to 13115). We aligned all data to the first session after reaching expert criterion (70% performance) and divided them in a naïve phase (sessions before the first expert session) and an expert phase (Figure 1B). The average number of trials in both naive and expert phase was around 3500 (3419 ± 2950 and 3556 ± 1638 trials, respectively; mean ± s.d.).

Tracking mouse behavior throughout learning

Several behavioral variables were measured throughout task learning. Whisker movement and whisker-to-texture touch times were quantified with markerless pose estimation using DeepLabCut (DLC) (Mathis et al., 2018). Videos of the mouse’s snout and whiskers were recorded from the top at 40 frames per second with a camera (A504k; Basler) and 940-nm infra-red LED illumination. We trained DLC to track five body features: (1) the nose tip; (2) the base of the presumed C2 whisker and (3) the whisker shaft position 10 mm away from the base; (4) the bottom corner of the approaching texture; and (5) the intersection of the whisker tip and the texture aiming to quantify the first touch moment. Tracked across frames, the extracted coordinates were used only if DLC confidence was high (p ≤ 0.05). If only one frame had a low confidence the missing value was interpolated as the average of the coordinates from the previous and the following frame. Longer periods of low confidence were assigned NaNs. We calculated the whisker angle from the dot product of two Euclidean vectors (vector nose to whisker-base) and (vector from whisker-base to the 10-mm shaft position). Next, we transformed the time-dependent whisker angle into whisking amplitude envelope by calculating the maximum minus minimum amplitudes in a 250-ms sliding window. Whisker-to-texture touch was expressed as a binary value for every frame during the trial (0 for p > 0.05, and 1 for p ≤ 0.05). The average of this vector across trials represented the touch probability. The touch probability increased as the texture came into reach of the whiskers and was close to 1 when the texture was standing next to the whiskers. Licking was quantified as event rate from the piezoelectric lick sensor. All behavioral variables were resampled to 20 Hz, the rate of photometry recordings. Pearson correlation coefficient was calculated as the covariance of two inputs normalized by the product of their standard deviations. First input represented calcium signals from all Hit trials concatenated across the session and the second input was similarly structured respective behavioral variable. To convert whisker-to-texture touch into continuous variable for the trial-based correlation analysis, binary values of whisker-to-texture touch were convolved with Calcium-like kernel function of 500ms decay time. Only correlation coefficients with p ≤ 0.05 were used in further comparison for naïve and expert phase (correlation coefficients with p > 0.05 were set to NaN) and later averaged within naïve and expert sessions for every mouse and tested with a Mann–Whitney U test.

Multi-fiber photometry setup

We used an Omicron LuxX 473-nm laser for excitation of GCaMP6m and a Coherent OBIS LS 561-nm laser to excite R-CaMP1.07. To achieve stable CW operation, lasers were run at 80% of maximal output power. A variable neutral density filter (NDC-25C-4M; Thorlabs) was used to reduce fluorescence excitation power to ∼1.3 mW/mm2 per fiber channel at the fiber tip. One or two cylindrical lenses shaped the excitation beam in an appropriate illumination pattern at the object plane of the objective. First, we expanded and collimated the circular beam with an achromatic Galilean beam expander (GBE05-A; Thorlabs). Second, to create a line illumination pattern matching the 12-fiber array we used a 75-mm focal length cylindrical lens (LJ1703RM-A; Thorlabs) placed at ∼145 mm distance from the objective (TL4x-SAP; Thorlabs). For the 48-fiber array a rectangular illumination pattern was required. For this purpose, we added a cylindrical lens with 150-mm focal length (LJ1629RM-A; Thorlabs) oriented 90° with respect to the 75-mm focal length cylindrical lens. A dichroic beamsplitter (F58-486 dual line, AHF) coupled the excitation light for both wavelengths (473 nm or 561 nm) into the objective (TL4x-SAP, Thorlabs) and transmitted the fluorescence in the emission spectral windows for GCaMP6m and R-CaMP1.07. To separate fluorescence signals from residuals of the excitation light and to minimize auto-fluorescence generated in a broader spectral range we used emission filters (525/50 nm, F37-516, AHF, for GCaMP6m; and 605/70 nm, F47-605, AHF, for R-CaMP1.07) and a triple-line notch filter at 425, 473 and 561 nm (ZET405/473/561 for Omicron LuxX 473nm and Coherent OBIS LS 561nm lasers, or, alternatively, for a combination of ZET405/488 ZET488/561 Chroma Technology Corp. for Coherent OBIS LX 488nm and Coherent OBIS LS 561nm lasers). To image the end face of the fiber array onto the camera we used a 200-mm focal tube lens (Proximity Series InfiniTube) with internal focusing. The image was created at the back focal plane of the tube lens on the CMOS sensor (ORCA Flash4.0, Hamamatsu camera). Calcium signals were expressed as percentage ΔF/F relative to the fluorescence baseline (average within the 0–1 s period before the first auditory cue).

Optogenetic manipulation of iGP

We used VGAT-ChR2 EYFP transgenic mice to activate GABAergic neurons in iGP. A 473-nm laser (Omicron LuxX) was used to excite GABAergic neurons. The beam was expanded (5× GBE05-A, Thorlabs), collimated and directed to the objective. The beam was coupled into the selected fiber channel corresponding to the implanted iGP target region by steering it with the mirror placed close to the microscope objective’s back aperture (Sych et al., 2019). We applied a 4-s long period of 473-nm stimulation from 2-6 s after trial start using 20-ms pulses at 20 Hz produced by a waveform generator (Agilent 33500B; TTL-triggered from the DAQ board for behavior control). To avoid interference of the optogenetic laser illumination with the photometric R-CaMP1.07 measurements of brain activity at 561-nm illumination, we quantified the mean ΔF/F activity in Hit, Miss/LaserOff and Miss/LaserOn trials (Figure 5) after the texture presentation period (during 6.5 - 7 s of trial time).

Post-hoc immunohistochemistry

We followed protocols for brain slice preparation as previously described (Sych et al., 2019). Mice were anaesthetized (100 mg /kg bw ketamine and 20 mg /kg bw xylazine) and perfused transcardially with 4% paraformaldehyde in phosphate buffer, pH 7.4. After perfusion, tissue was removed from the skull and the head including the multi-fiber implant was additionally fixated in 4% paraformaldehyde for one week. Then, the ventral (bottom) side of the skull bone was removed, and the brain was carefully extruded. Coronal sections (75 – 100 μm thickness) were cut with a vibratome (VT100, Leica). Stained sections were mounted onto glass slides and confocal images were acquired with an Olympus FV1000. Sequential imaging of coronal slices allowed us to track fiber shafts and to verify the 3D positions of the fiber tips. The final assignment of each fiber channel was completed after the fiber tip positions were tracked on the histological slices and aligned to the Mouse Brain atlas.

Quantification and statistical analysis

We used Matlab R2021b (MathWorks) for the statistical analyses. Mann–Whitney U test was used throughout the paper with non-matched groups. Wilcoxon signed-rank test was used with matched groups. Bar plots and plots with shaded errors express mean ± SEM. All of the statistical details of experiments can be found in the figure legends. To compute TE for our data we have used the IDTxl Python package (Wollstadt et al., 2019).

Effective connectivity estimated by transfer entropy

In order to compute connectivity (the amount of cross-regional interaction) we used transfer entropy (TE) as a directed information-theoretic metric for estimating the statistical dependence between the past activity state of a source region and the present activity state of a target region (Schreiber, 2000; Wollstadt et al., 2019; Novelli et al., 2019). In its basic form, the Bivariate TE (BTE) from a source variable to a target variable is defined as the Mutual Information (MI) between the current value of , and the past value of , conditioned on the immediate past values of .where and . Given the source and target time windows , TE considers the multiple past time-delayed values of the source and target (from to ) on the TE. One can also consider a time delay u between the source and target variables (i.e., instead of ) and select the value of u that maximizes the TE. We will refer to this value as the time lag, which approximates the timescale of the effective connection. The true probability distributions of the variables and (e.g., neuronal signals) are not known a priori and must be estimated from repeated sampling of and over multiple trials, or over several consecutive time bins if the signal is stationary in the corresponding time window. In general, robust estimation of information-theoretic quantities such as MI and TE from sample data is not trivial, for example, the dimension of the space increases exponentially with L. We have used the IDTxl Python package (Wollstadt et al., 2019) to compute TE for our data. The BTE algorithm implemented in IDTxl selects the most informative subset of past variables to maximize the TE and thus forms an optimal embedding for each variable (independently of the others). BTE has two major advantages over simpler measures of functional connectivity such as cross-correlation. First, it can handle arbitrary nonlinear relationships provided sufficient data are available and, second, it explicitly eliminates any spurious connections that are already well-explained by the past of the target variable. However, it only considers two variables at a time, and thus it is overly liberal when addressing higher-order interactions such as causal chains or common causes. To control false positive rate due to high-order interactions, we used the Multivariate TE (MTE) algorithm provided by IDTxl, which builds a minimal network incrementally by adding one connection at a time, as long as adding a new connection further improves the prediction with respect to the previously inferred connections (Novelli et al., 2019). Due to redundancy in recorded neural data, the true effective connectivity is likely a compromise between the minimal network provided by MTE and the denser one inferred by BTE. To capture time-dependent changes in effective connectivity, we computed BTE and MTE across trial time in a sliding window of six 200-ms time bins (L = 5).

Connection strength

The connectivity matrix for every session was calculated with the IDTxl (Wollstadt et al., 2019) library. We defined a binary connection between each brain region based on a significance level, i.e. a connection was set to 1 (‘present’) if p ≤ 0.01 and to zero otherwise (‘absent’). Connection strength was defined as the mean frequency of significant TE bins in a given analysis window (stimulus or action-reward), averaged across all mice and all expert sessions (n = 13 mice for Figure 4 and n = 4 mice for Figure 6). To estimate differences in effective connectivity between Hit and CR trials (Figure 4G), we used a connectivity matrix inferred from MTE. For every mouse we averaged across all expert sessions. The resulting distributions of connection strength across mice for Hit and CR trials were tested with the Wilcoxon signed-rank test at every bin of trial time. Connections with p ≤ 0.05 were shown as significantly different for Hit vs CR trials.

Analysis of functional network stability across sessions

In order to evaluate the stability of emerging effective connections, we designed a test based on the probability distribution for the number of the effective connections shared between subsequent sessions. The number of connections observed in any given session is computed by summing up the entries of the corresponding binary connectivity matrix M.

Note that the diagonal values of M are zero because, by definition, only connections between different channels are considered. Thus, the highest possible number of connections is given by the number of off-diagonal entries of the matrix M, namelywhere is the number of matrix rows (that is, the number of recorded brain regions).

The number of shared connections N is defined as the sum of entries of the overlap between the matrices of two subsequent sessions. Denoting the "previous" session with x and the "next" session with y,

It can be easily seen that the maximum number of shared connections is given by the number of connections of the session that has the smallest of them

Given prior knowledge of and , the constraints on possible values of can be summarized in the following: (1) is non-negative; (2) does not exceed the maximum number of possible connections; (3) does not exceed the number of connections lowest among the two matrices; (4) given sufficiently large and , there is a guaranteed non-zero overlap.

We construct the following null hypothesis : For all sessions the effective connectivity is completely random, except for the total number of connections . Namely, under all permutations of connections among possible connections are equiprobable. We proceed to calculate the probability distribution of the number of shared connections under , which will be used to test the significance of shared connections observed in real data.

The denominator is given by the number of possible permutations of and . The numerator is constructed by first choosing from , then choosing from , and finally distributing the non-shared connections of among the remaining zeros of .

We use the above distribution to calculate the p-values of the following two-sided tests: (1) , the probability of getting as few or less shared connections as observed; (2) , the probability of getting as many or more shared connections as observed. Calculating the above p-values for observed data demonstrates that the number of the observed shared connections is consistently very unlikely under , proving that does not hold for the observed data. This implies that the overlap between subsequent sessions is significantly different from a purely random one, mostly because there are many more shared connections than expected under .

Further, we use as a measure of stability of the network. The larger the absolute value of , the more shared connections there are, and the more stable is the network. We plotted as a scatter plot of mouse performance in order to demonstrate the relationship between the network stability and learning (Figures 4D, 4E, 6C, and 6D). The scatter plots combined the results of all consecutive session pairs from all mice. We also performed a Mann–Whitney U test to determine whether the networks were more stable for expert or naive mice. We visualized the results of this test using violin plots. Empirical distributions of are highly asymmetric and thus are difficult to accurately represent using bar plots. Finally, we used a sliding window to construct a trendline for the scatter plot to show the dynamic changes of stability with learning. For any given window, the likelihood of observing all data points in that window is given byand the associated log-likelihoodwhere i was used to index the data points within the sliding window. The log-likelihood is a linear function of the number of data points in the window , so we normalized the log-likelihood with that number to account for the fact that the number of data points in the window may vary. Thus, we defined the trend line as the arithmetic mean of the logarithms of individual probabilities, which translates to the geometric mean of the probabilities themselves.

In order to smoothen the trend line, we replaced the window-based estimation of running average with kernel smoothing. We used a Gaussian kernel with units of performance squared.

Clustering coefficient

We used a local clustering coefficient as a measure of local change in connectivity during learning. For every region of the connectivity matrix estimated from the MTE analysis we calculated the number of outgoing connections nout. If nout ≥ 2 and both target regions were also connected, we calculated the clustering coefficient aswhere we summed across all outgoing connections for the source to pairs of targets and normalize it by the node degree and number of targeted regions minus one. If we compared the changes in clustering coefficient for the stimulus window and the action-reward window, the connectivity matrix is first averaged across all connections within the given time interval and then the weighted clustering coefficient was calculated, where represent averaged across the time interval connection strength. Next, clustering coefficient within each time interval is averaged across respective groups for naïve and expert sessions (Figure 4B) or for the groups for Hit and CR trials (Figure 5B). Group differences are estimated with a Mann–Whitney U test across mice.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Bacterial and virus strains
AAV2.9-hSyn-GCaMP6m	Addgene	Plasmid #100841 (https://www.addgene.org/100841/)
AAV1-EF1a-R-CaMP1.07	Viral Vector Facility UZH/ETH Zurich	AAV vector plasmid p78 (vvf.uzh.ch)
Experimental models: Organisms/strains
C57BL/6 mice	The Jackson Laboratory	Bred in house but originally from Jackson LaboratoryStrain #014548, RRID: IMSR_JAX:000664 (https://www.jax.org/strain/000664)
VGAT-ChR2 EYFP transgenic mice	The Jackson Laboratory	Bred in house but originally from Jackson LaboratoryStrain #014548, RRID: IMSR_JAX:014548 (https://www.jax.org/strain/014548)
Software and Algorithms
DOI	https://doi.org/10.5281/zenodo.7016970
URL	https://github.com/YaroslavSych/CBTC_Network/tree/v1.1