Goal-directed behavior requires sensory information to be transformed into motor commands, yet the precise temporal structure of communication between sensory and motor cortices during this process remains poorly understood. By simultaneously recording spiking activity from primary somatosensory cortex (S1) and secondary motor cortex (M2) in mice performing a tactile delayed go/no-go task, we identified a brief, stimulus-triggered communication window that predicted behavioral success. Time-shifted canonical correlation analysis revealed a sharp, feedforward surge in S1 → M2 coupling peaking at a short latency, as expected for a potentially direct transmission pathway. This transient interaction conveyed stimulus-selective information such that the S1 → M2 communication subspace aligned with stimulus-coding axes in both regions. In contrast, M2 → S1 feedback interactions were weaker and delayed. Brief optogenetic suppression of S1 reduced both M2 stimulus selectivity and behavioral discriminability, consistent with a functional role of the transient S1 → M2 coupling. These findings demonstrate that sensorimotor transformation is supported not by sustained linear coupling along the stimulus- or choice- axis, but by a short-lived and highly selective feedforward communication channel that rapidly transmits reliable sensory representations into motor-planning circuits.

Behavioral task and dual-site recordings

We investigated interactions between primary somatosensory cortex (S1) and secondary motor cortex (M2) by simultaneously recording neuronal activity in mice performing a tactile delayed go/no-go task (Fig 1A–1E). In this task, a brief air puff was delivered unilaterally to the whiskers, followed by a 1-s delay, a 2-s auditory-cued response window, and a 5.1–7.0-s inter-trial interval. In go trials, a left-sided air puff prompted mice to extend their tongue, break an IR beam, and receive a water reward (Fig 1A and 1B). No-go trials had the same events and timing as go trials, except that the air puff was delivered to the opposite side and mice were required to withhold licking (Fig 1B and 1C). Mice were trained to stable performance before recordings (S1A–S1C Fig) and typically achieved d-prime (d′) >1.5 during electrophysiological sessions (S1D Fig).

Recordings targeted S1 and M2 subregions within the whisker sensorimotor subnetwork [49] (Figs 1D, S1E and S1F). As shown in two example sessions, both regions exhibited strong initial responses to contralateral stimuli and weaker responses to ipsilateral stimuli (Figs 1E and S2A–S2H). In S1, the contralateral stimulus (both hit and miss trials) evoked reliable initial responses within 50 ms of stimulus onset in many neurons [52–54] (Fig 2A, 2B, 2E and 2F). In M2, population activity also showed strong initial responses to contralateral stimuli (hit trials), whereas responses to ipsilateral stimuli (CR trials) were weaker (S2C, S2D, S2G and S2H Fig). To quantify stimulus preference, we calculated a selectivity index for each neuron using the first 100-ms time window following stimulus onset. In this time window, M2 populations had a much narrower range of selectivity than S1 (S3A Fig). Across 50 recording sessions, the median fraction of S1 neurons with significant selectivity within 100 ms of stimulus onset was 0.29. We used this value to classify sessions into high- or low-selectivity (S3B–S3G Fig), reflecting whether the air puff effectively engaged the S1 receptive fields.

Transient S1–M2 communication predicts behavioral performance

We employed time-shifted CCA [43] to characterize population-level communication between S1 and M2 (Fig 1F and 1G). This analysis identifies feedforward and feedback interactions as well as any task-specific synchronizations between the two regions [43,55]. To understand overall interactions, including shared tuning properties between the two regions, we first computed total correlation, which preserves all possible interactions, and compared it with “noise correlation,” which reflects trial-to-trial co-fluctuations in activity between S1 and M2 by subtracting each neuron’s trial-averaged activity [41,43,48]. We first examined correlation maps of the two example sessions introduced earlier (in S2 Fig). Total correlation revealed focal, synchronous interactions shortly after the stimulus (Fig 2A and 2B, top left). Enhanced, transient interactions were also observable when go (left air puff) and no-go (right air puff) trials were separately analyzed (Fig 2A and 2B, middle and bottom left). In noise correlation, we observed much weaker stimulus-locked interactions compared to the total correlation (Fig 2A and 2B, right).

Across 50 sessions, we observed a significant transient increase in total correlation between S1 and M2, followed by a prominent decrease in the canonical correlation value (CCV) relative to the baseline period (−0.25 to 0 s) (all trials: Fig 2C, p_up = 1.9 × 10⁻⁷, z_up = 5.2, and p_down = 6.2 × 10⁻⁷, z_down = −4.9; go trials: Fig 2D, p_up = 2.0 × 10⁻⁷, z_up = 5.2 and p_down = 3.4 × 10⁻⁷, z_down = −5.1; no-go trials: Fig 2E, p_up = 0.00015, z_up = 3.8 and p_down = 6.6 × 10⁻⁸, z_down = −5.4; upper and lower confidence interval (CI) crossing points up to 1 s post-stimulus versus baseline, Wilcoxon signed-rank test). CCV returned to baseline levels during the delay period and again dropped below the lower CI after the auditory response-cue presentation in all types of trials. Considering noise correlations, the transient CCV elevation was modest but significant when all trials were pooled and in no-go trials (all trials: Fig 2F, p_up = 0.0046, z_up = 2.8 and p_down = 5.9 × 10⁻⁷, z_down = −5.0; go trials: Fig 2G, p_down = 3.5 × 10⁻⁷, z_down = −5.1; no-go trials: Fig 2H, p_up = 0.034, z_up = 2.1 and p_down = 5.6 × 10⁻⁸, z_down = −5.4). Beyond this initial transient, temporal fluctuations in noise correlations closely mirrored those in total correlations. Because total correlation reflects both stimulus-driven shared tuning (signal) and trial-to-trial (noise) components, this similarity of temporal fluctuations indicates that the contribution of shared stimulus tuning between S1 and M2 rapidly diminished after the initial response. Thus, this implies that the two populations showed largely independent tuning once the transient sensory-evoked alignment subsided.

To further relate these interactions to behavior, we analyzed trial outcome types separately (S4 Fig). Correct trials (hit and CR) exhibited interaction patterns like those observed in go and no-go trials (S4A, S4B, S4E and S4F Fig), whereas error trials showed limited fluctuations in CCVs (S4C, S4D, S4G and S4H Fig). This contrast suggests that the observed interaction dynamics reflect goal-directed communication during successful performance, although limited sample size in the error trials precludes a fair comparison with the correct trials. The stronger and more sharply timed transient observed in total correlations relative to noise correlations suggests that the stimulus onset triggers a stimulus-gated S1–M2 coupling that exceeds spontaneous coupling and shared variability not directly attributable to stimulus drive. In other words, this early interaction reflects a rapid, stimulus-evoked alignment of population activity across the two regions, which is then followed by slower modulation of internally generated shared variability.

One possible explanation for the weak coupling during the delay period is variability in reaction times within and across sessions, which could blur together lick-evoked drops in CCV, and potential delay-period rises across trials. To account for this within-session variability, we realigned the spiking data to response onset (first lick) and recomputed CCV to examine delay-period coupling using total correlation (S4I and S4J Fig). We did not observe a significant increase in CCV before the first lick, although this could still reflect a mixture of stimulus-driven and delay-period coupling. To further isolate delay-period coupling from stimulus-evoked one, we excluded sessions with short reaction times (mRT ≤ 1.2 s), in which animals typically responded before the auditory cue. Even in sessions with long reaction times (mRT > 1.2 s, n = 33), CCV did not increase during the delay period, and it instead decreased below CI upon response onset (S4J Fig, total: p_down = 5.4 × 10⁻⁷, z_down = −5.0; noise: p_down = 7.8 × 10⁻⁷, z_down = −4.9). Notably, temporal fluctuations in CCV could not be explained by changes in firing rate. We assessed the geometric mean firing rate (gmFR) [56–58] across regions, which reflects the overall task-induced activity levels in the two populations irrespective of their trial-by-trial co-variation (Fig 2C–2H). Although gmFR across regions increased following response-cue onset, this was not accompanied by a corresponding rise in CCV, indicating that transient interareal coupling dynamics are not simply a reflection of global increases in population activity across regions.

Transient CCV increases may reflect behaviorally relevant communication between stimulus-selective populations within S1 and M2. Supporting this view, the magnitude of average CCV within 100 ms from stimulus onset was significantly correlated with behavioral discriminability (d-prime) when considering all trials (total: Fig 3A and 3D; r = 0.63, p = 2 × 10⁻⁵, noise: 3G and 3J; r = 0.53, p = 5 × 10⁻⁴, Pearson’s correlation) and when restricting analysis to go trials (total: Fig 3B and 3E; r = 0.60, p = 4.7 × 10⁻⁵, noise: 3H and 3K; r = 0.39, p = 0.015), with no-go trials showing no reliable relationship (total: Fig 3C and 3F; r = 0.14, p = 0.41, noise: 3I and 3L; r = −0.03, p = 0.86). This pattern indicates that the early, stimulus-locked interaction is closely tied to successful performance, particularly in trials that require an active decision to respond. By contrast, neither the post-stimulus decrease in CCV (0.1–0.3 s after stimulus onset) nor the later CCV increase during no-go trials (0.9–1.2 s) was significantly correlated with behavioral performance (total post-stimulus CCV: p = 0.1, r = 0.27; p = 0.052, r = 0.31 and p = 0.7, r = 0.056 for all trials, go and no-go trials; noise post-stimulus CCV: p = 0.1, r = 0.27; p = 0.053, r = 0.31 and p = 0.7, r = 0.055 for all trials, go and no-go trials; total later CCV: p = 0.06, r = 0.3; p = 0.3, r = 0.17, and p = 0.99, r = −0.002 for all trials, go and no-go respectively; noise later CCV: p = 0.06, r = 0.3; p = 0.3, r = 0.17, and p = 0.97, r = −0.005 for all trials, go and no-go respectively). These later fluctuations, therefore, likely index arousal-related, global changes in cortical state, rather than performance-sensitive communication. Especially, such decreases in correlated variability are consistent with the degradation of high-amplitude, synchronized activity often observed in quiet wakefulness [59–62].

We further examined the relationship between transient interareal interactions and behavioral performance (d-prime) while accounting for session-to-session variation in neuronal selectivity using a linear mixed-effects model. In this model, behavioral performance was predicted by three factors: the average CCV within 100 ms of stimulus onset, the fraction of significantly selective neurons in M2, and the fraction of significantly selective neurons in S1, with subjects’ identity included as a random effect. Consistent with the Pearson’s correlation analyses described above, the model showed a significant positive coefficient for average CCV (p = 0.001, t = 3.6 and df = 31 for total correlation using all trials; p = 0.004, t = 3.1 and df = 31 for noise correlation using all trials, n = 50 sessions from 18 mice). In contrast, population selectivity in either area did not significantly predict behavioral performance (total: M2; p = 0.35, t = −0.94, S1; p = 0.34, t = −0.96, noise: M2; p = 0.97, t = −0.04, S1: p = 0.69, t = 0.41, df = 31). No significant interactions were found among the three predictors for either total or noise correlation. Even in reduced single-predictor models, population selectivity in S1 or M2 alone failed to predict behavioral performance (S1: p = 0.27, t = 1.1, M2: p = 0.34, t = 0.96, df = 37). These results indicate that the transient interareal interactions between S1 and M2 are a robust predictor of behavioral performance, whereas variability in population selectivity likely arises from receptive field variations among sessions and does not account for behavioral variability.

We consolidated these results with a trial-resolved analysis to test whether moment-to-moment fluctuations in inter-areal coupling predict the behavioral outcome of each trial. To minimize the influence of fluctuations in arousal or task engagement, we first estimated the animal’s moment-to-moment performance level throughout the session using a local hit-rate measure and excluded trials occurring during periods of very low engagement (local hit rate < 0.2). This procedure preferentially removed trials in which animals were likely disengaged from the task, thereby reducing potential confounds related to arousal and task engagement. Using the remaining engaged trials, we computed trial-level coupling; For each trial, we quantified coupling strength as the average correlation between the first pair of canonical variables within 100 ms of the air puff in either the feedforward (FF: 2–20 ms) or feedback (FB: −20 to −2 ms) direction. We then applied a binomial mixed-effects model with trial outcome (correct versus incorrect) as the response variable, trial-level coupling as a fixed effect, and sessions and subjects as random effects.

In high-selectivity sessions, trial-level coupling significantly predicted the trial outcome only in the feedforward direction (FF: p = 0.0045, t = 2.84, FB: p = 0.24, t = 1.2, df = 5,242, 25 sessions from 12 mice). In low-selectivity sessions, there was no significant association between trial-level coupling and trial outcome in either direction (FF: p = 0.80, t = −0.26, FB: p = 0.53, t = 0.62, df = 5,065, 25 sessions from 13 mice). Consistent with this result, including selectivity (S1sel) as a fixed-effect predictor revealed a significant, positive interaction between trial-level feedforward coupling and the fraction of selective S1 neurons in each session (FF: p_{coupling × S1sel} = 0.017, t_{coupling × S1sel} = 2.38, p_coupling = 0.44, t_coupling = −0.76, p_S1sel = 0.35, t_S1sel = −0.93, FB: p_{coupling×S1sel} = 0.21, t_{coupling × S1sel} = 1.26, p_coupling = 0.92, t_coupling = −0.11, p_S1sel = 0.89, t_S1sel = 0.14, df = 10,307, 50 sessions, S5 Fig), while selectivity in M2 (M2sel) did not interact with the trial-level coupling (FF: p_{coupling × M2sel} = 0.90, t_{coupling × M2sel} = 0.13, FB: p_{coupling × M2sel} = 0.46, t_{coupling × M2sel} = −0.74). Together, these analyses indicate that stronger early feedforward S1–M2 coupling is associated with successful behavior within trials, specifically when sensory representations in S1 are strong.

Because the canonical weights used for this analysis were estimated from all trials within each session, we repeated it with 10-fold cross-validation, estimating weights from 90% of trials and applying them to the held-out 10%. For total coupling, the association between feedforward coupling and trial outcome was preserved (go trials of high-selectivity sessions: p = 0.006, t = 2.74, df = 2,846). For noise coupling, cross-validation reduced the estimated coupling roughly 5-fold (mean 0.184 versus 0.038) and the association with trial outcome was no longer detectable (p = 0.92, t = −0.10, df = 2,847), indicating that the canonical subspace estimated from residual activity within a single session does not generalize to held-out trials. The trial-level noise-coupling results reported here, including the delay-period analyses below, should therefore be read as a within-session characterization rather than an independently validated prediction.

To test the possibility of a link between trial-level coupling and behavioral state, we regressed the trial-level coupling with a proxy of state. Although pupil dynamics around whisker stimulus are likely too slow to be used as a predictor of the transient coupling [63], reflexive movements of the facial area in response to stimulation [64] could be used as a proxy of the state with sufficiently short latency. We therefore included an analysis incorporating the motion index [65] as a proxy of state to predict trial-level coupling in high-selectivity sessions where transient trial-level coupling predicted the outcome. Motion index quickly increased upon air puff presentation, indicating reflexive movements in response to the puff (S6A Fig) as reported elsewhere [64]. We asked whether trial-level coupling in three distinct periods can be related to the motion index as a state variable during 3 periods of: (1) baseline; −100 to 0 ms from air puff, (2) transient; 0 to 100 ms from air puff and (3) decorrelation; 120 to 220 ms from air puff. Trial-level coupling was computed as detailed previously. S6B Fig shows the relationship between motion index and trial-level coupling averaged over the first 100 ms after air puff onset (5,244 trials from 25 high-selectivity sessions). After z-scoring both variables within each session, no significant correlation was observed (r = −0.02, p = 0.145, Pearson’s correlation).

Using a linear mixed-effects model with session and subject included as random effects, we tested whether trial-level coupling strength is predictable from the motion index during each of the three periods. Across all trials from high-selectivity sessions, neither feedforward nor feedback coupling showed a significant association with motion index during any of the three periods (baseline; FF; p = 0.53, t = −0.62, FB; p = 0.87, t = 0.16, transient; FF; p = 0.19, t = 1.3, FB; p = 0.91, t = 0.11, decorrelation; FF; p = 0.87, t = −0.17, FB; p = 0.23, t = 1.2, df = 5,242). Thus, trial-to-trial fluctuations in movement state did not predict fluctuations in inter-areal coupling strength despite the rise of both variables during the first 100 ms of the air puff based on the grand average data (S6C Fig). We further examined whether the motion index during the transient period contributed to behavioral performance, either independently or through an interaction with trial-level coupling strength. Trial-level motion index alone was not associated with performance (p = 0.87, t = −0.16, linear mixed-effects, df = 5,242). When both motion index and coupling strength were included in the model, only feedforward coupling during the transient epoch remained a significant predictor of behavioral performance. Neither motion index nor its interaction with coupling contributed significantly, and no significant effects were observed for feedback coupling (FF: p_motion = 0.68, t_motion = −0.42, p_{motion x coupling} = 0.71, t_{motion × coupling} = 0.37, p_coupling = 0.016, t_coupling = 2.4, FB: p_motion = 0.24, t_motion = −1.2, p_{motion x coupling} = 0.25, t_{motion × coupling} = 1.2, p_coupling = 0.61, t_coupling = 0.52, df = 5,240). These results indicate that trial-to-trial fluctuations in movement-related state do not account for the observed variability in inter-areal coupling. Instead, they support the interpretation that the transient increase in S1–M2 coupling reflects specific functional interactions between the two regions rather than nonspecific modulation due to motion-related activity.

We found no significant rise in the delay-period coupling even when the data was aligned to the response time (S4J Fig). However, it is possible that S1–M2 coupling in the delay period is highly variable from trial to trial in a way that could predict performance or response time, while it is masked by the session-level analysis. To address this possibility, we investigated whether trial-level delay-period coupling could predict future licks (i.e., reaction time) or trial outcome (correctness). We computed trial-level coupling as previously described, except that we chose larger, non-overlapping windows to account for potentially slower coupling (window size = 250 ms, step = 250 ms, bin size = 2 ms), using trials with the reaction time between 1.25 and 6.25 s. First, we tested whether the trial-level delay period coupling (noise correlations, 0.25–1.25 s from air puff) can predict the trial outcome in a binomial mixed-effects model with session and subjects as random effects. We did not find a significant association between trial-level, delay-period coupling and trial outcome in either direction (FF: p = 0.13, t = 1.5, FB: p = 0.49, t = 0.70, df = 6,853, 36 sessions from 16 mice). Including S1 selectivity (S1sel) as a fixed effect did not change the story (FF: p_coupling = 0.54, t_coupling = 0.61, p_S1sel = 0.73, t_S1sel = −0.35, p_{coupling × S1Sel} = 0.78, t_{coupling × S1sel} = 0.28, FB: p_coupling = 0.81, t_coupling = −0.24, p_S1Sel = 0.46, t_S1sel = −0.73, p_{coupling × S1Sel} = 0.40, t_{coupling × S1sel} = 0.83, df = 6,851).

Furthermore, in a linear mixed-effects model where trial-level, delay-period coupling and selectivity were used to predict reaction time, we found no significant association (FF: p_coupling = 0.63, t_coupling = 0.49, p_S1Sel = 0.33, t_S1sel = −0.98, p_{coupling × S1Sel} = 0.93, t_{coupling × S1Sel} = 0.09, FB: p_coupling = 0.059, t_coupling = 1.9, p_S1Sel = 0.64, t_S1Sel = −0.46, p_{coupling × S1Sel} = 0.59, t_{coupling × S1Sel} = −0.53, df = 4,274, 36 sessions from 16 mice). However, trial-level delay-period coupling as the sole fixed-effect variable significantly predicted reaction time in the feedback direction (FB: p = 0.020, t = 2.3, FF: p = 0.36, t = 0.92, df = 4,276) with a positive weight (β = 0.68). This relationship was maintained when using broader feedback delay (−50 to −20 ms) instead of the original feedback delays (−20 to −2 ms) for trial-level coupling, indicating that it could arise from indirect feedback from M2 to S1. Indeed, this relationship was not limited to this period and these sessions with longer reaction times. Across all sessions, a significant positive association existed between trial-level coupling and reaction time (−100–0 ms from air puff; p = 0.043, t = 2.0, 0–200 ms from air puff; p = 0.003, t = 2.9, df = 6,213, 50 sessions). Its positive relationship with reaction time indicates that the feedback from M2 to S1 may help to maintain an interacting loop between sensory and motor regions, to delay decision until further evidence (e.g., response cue) is available, or serve as a top-down learning signal [66] which relays feedback based on the recent decisions and sensory cues.

Overall, beyond the cross-session correlation analysis, these trial-resolved analyses indicate that stronger early feedforward S1–M2 coupling predicts successful behavior within trials, specifically when sensory information is robustly represented in S1, although the noise-coupling version of this analysis relies on subspaces estimated within session and was not recovered under cross-validation. Furthermore, our analyses showed that trial-to-trial fluctuations were not associated with fluctuations in motion-related behavioral state.

S1 and M2 exchange stimulus selectivity information in brief interactions

Interareal interactions reflecting direct communications must occur within short delays between the activities of the two regions [43]. To examine the directionality of these transient interactions during the air puff, we calculated feedforward and feedback components in short delays using all trials in each session for CCA (<15 ms, Fig 4A–4C). On average, the transient rise of CCV was significant only in the feedforward direction (noise: S1 → M2, Fig 4B, p_up = 0.0003, z_up = 3.7, Wilcoxon signed-rank test on upper CI crossing points relative to the baseline window, −50 ms to 0). In the feedback direction, CCV instead dropped significantly below baseline during the air puff, followed by a transient rebound to baseline and a subsequent drop (noise: M2 → S1, Fig 4C, p_down = 5.6 × 10⁻⁵, z_down = −4.0, Wilcoxon signed-rank test on lower CI crossing points (average) relative to baseline, n = 50 sessions).

To test whether the observed interactions carried stimulus-related information, we computed a coding direction () [67,68] using all trials in each session, during the first 100 ms of the stimulus presentation by subtracting the mean right-trial activity vector across neurons from the mean left-trial vector within each area. The resulting vector captures the population pattern, distinguishing between the two types of stimuli. Next, we computed another vector, a communication subspace (), as the weight of the first canonical dimension in each area. Unlike , which is fixed once defined, is time-varying and represents the axis along which S1 and M2 activities are most correlated at each moment. We quantified the similarity between and by calculating the cosine of the angle between the two vectors (i.e., projection magnitude of unit vectors, henceforth called “similarity”, Fig 4D and 4G). In the feedforward direction, both S1 and M2 similarity magnitudes rose significantly above baseline in parallel with the CCV (Fig 4E and 4H; p_up = 1.5 × 10⁻⁵, z_up = 4.3 and p_up = 6.8 × 10⁻⁵, z_up = 4.0, respectively). However, in the feedback direction, there was no significant rise in the similarity in either area (Fig 4F and 4I). Total correlations instead revealed significant coupling in both directions (S7A–S7I Fig). However, interactions were consistently weaker in the feedback direction (total; CCV: p = 5.7 × 10⁻⁹, z = 5.8, S1 similarity: p = 3.4 × 10⁻⁵, z = 4.1, M2 similarity: p = 1.0 × 10⁻⁵, z = 4.4, noise; CCV: p = 1.3 × 10⁻⁸, z = 5.7, S1 similarity: p = 2.2 × 10⁻⁵, z = 4.2, M2 similarity: p = 3.1 × 10⁻⁶, z = 4.7, Wilcoxon signed-rank test within 100 ms of air puff).

Is the early communication pattern originating from S1 required for successful task performance, and do these patterns causally shape early stimulus selectivity in M2? To address these questions, we transiently suppressed S1 during the early sensory period using optogenetics in PV-ChR2 mice (S8A Fig). Suppression of the right S1 significantly reduced behavioral discriminability (d-prime) compared with both the control condition and left S1 suppression (S8B Fig; P = 0.015, F = 9.2, df = 2, Control versus Right S1: P = 0.019, Left S1 versus Right S1: P = 0.042, 4 mice, repeated-measure ANOVA followed by Fisher’s LSD tests). Furthermore, recordings from the right M2 identified 147 neurons exhibiting significant stimulus selectivity within the first 100 ms after the air puff. Brief S1 suppression significantly decreased this early M2 selectivity (S8C Fig, p = 0.0017, t = −3.17, df = 292, 9 sessions from 4 mice, linear mixed effects). These results demonstrate that early S1 activity is required for both the emergence of robust stimulus selectivity in M2 and successful performance of the go/no-go task.

To assess whether the inclusion of multi-unit activity influenced the coupling dynamics, we repeated this analysis using the high-quality fraction of the units in both areas. Because a high inter-spike-interval (ISI) violation rate is commonly interpreted as evidence of poor unit isolation or contamination from spikes belonging to multiple neurons, we repeated the CCA analysis after excluding units with a high ISI violation rate and retaining ~80% of the units in each area (S9A Fig). This restriction enriches the dataset for better-isolated single units and reduces the likelihood that the observed delay reflects multiunit contamination or spike-sorting artifacts. Using the high-quality fraction of the units in both areas (~80%), we found highly consistent results regarding the dynamics of coupling across two areas for both types of correlation measures (CCV_noise; S9B Fig and CCV_total; S9C Fig).

Our analysis of the coding-communication alignment utilized the same trials for pd and cd estimations and therefore is prone to overfitting (i.e., alignment may arise from large co-fluctuations in limited trials). We therefore additionally performed a split-half cross-validation to independently estimate the communication and coding axes. We reproduced our similarity results in both feedforward and feedback directions using independent trials for pd and cd. Consistent with the original results, similarity in both regions increased significantly above baseline in the feedforward direction using either noise or total correlations (noise; S10A and S10B Fig; S1 similarity: p_up = 0.003, z_up = 2.9; M2 similarity: p_up = 0.033, z_up = 2.1; total; S10D and S10E Fig; S1 similarity: p_up = 5.8 × 10⁻⁶, z_up = 4.5; M2 similarity: p_up = 5.8 × 10⁻⁶, z_up = 4.5, upper confidence interval (CI) crossing points up to 100 ms post-stimulus versus baseline, Wilcoxon signed-rank test). In the feedback direction (S10C and S10F Fig), only total correlations showed a subtle but significant increase in similarity in both areas (S10F Fig; S1 similarity: p_up = 0.03, z_up = 2.2; M2 similarity: p_up = 0.01, z_up = 2.5). This cross-validated analysis shows that the alignment of the shared and local coding subspaces is generalizable to independent trials and therefore unlikely to arise from large fluctuations in limited trials. In other words, the shared activity captures coding properties of the population regardless of trial-to-trial variability in response. Furthermore, we repeated the CCA analysis using 10-fold cross-validation to control for potential overfitting of CCA itself. In each fold, trials were randomly divided into training (90%) and test (10%) sets. CCA weights were estimated using only the training trials, and the resulting canonical coefficients were then applied to the held-out test trials to compute the CCVs. Final CCVs were obtained by averaging across folds (S11 Fig; average cross-validated CCVs across 50 sessions). This cross-validated analysis produced results highly consistent with the original non-cross-validated CCA, including communication timing and coupling magnitude across delays. Overall, these analyses demonstrate that the observed feedforward communication and its alignment with the coding subspace are robust and are unlikely to arise from overfitting.

Our similarity measure is based on absolute values. We observed considerable sign flips from one time window to another within sessions and across sessions (Example sessions in S12A Fig). Importantly, such sign changes do not necessarily indicate a loss of alignment between coding and communication subspaces. Because cosine similarity is sensitive to vector orientation, a reversal in the dominant activity pattern between interacting populations can invert the sign while still preserving a strong subspace relationship. Thus, both positive and negative cosine values can reflect structured alignment between stimulus-coding and communication axes. To directly test whether sign ambiguity influenced our conclusions, we repeated the analysis using element-wise products of the signed similarity matrices from S1 and M2. The resulting map was highly similar to those obtained using absolute cosine similarity across sessions (S12B Fig; FF: p_up = 5.9 × 10⁻⁷, z_up = 5.0, FB: p_down = 0.18, z_down = −1.3, n = 50 sessions, upper and lower confidence interval (CI) crossing points up to 100 ms post-stimulus versus baseline, Wilcoxon signed-rank test), indicating that the reported effects were not driven by arbitrary sign flips. Rather, the transient increases in similarity reflect robust alignment between the stimulus-coding subspace and the communication subspace during the brief interaction window. Notably, the feedforward component of the product of similarity was significantly correlated with the sessions’ S1 selectivity, further emphasizing the specificity of the communication-coding alignment (0–100 ms; FF: p = 6.7 × 10⁻⁵, t = 4.4, FB: p = 0.63, t = 0.48, df = 48, linear mixed-effects).

Feedback interactions were observed only in total correlations. Indeed, in noise correlations, the feedback effect emerged as a small, delayed peak (50–100 ms) superimposed on the stimulus-evoked desynchronization. Notably, this enhancement was not associated with a significant rise in - similarities. During feedforward interactions, and were consistently aligned, implying alignment of signal and noise axes [41,69]. This is evident in consistently enhanced similarity values in both total and noise correlations following the air puff. In contrast, in the feedback direction of noise interactions, the lack of alignment between and implies that signal and noise axes are independent [41,69–71]. Therefore, the enhanced feedback in total correlations and the corresponding similarities indicate dominant co-fluctuation along a signal axis that is independent of noise correlations. Nevertheless, the presence of this feedback component is consistent with previous reports of feedback signaling between somatosensory and frontal areas [22,72].

These findings reveal that S1 and M2 quickly share task-relevant sensory information within a brief, stimulus-locked window, with communication confined to the earliest stimulus-evoked population pattern, rather than sustaining prolonged or recurrent communication. Moreover, the directionality of signal and noise correlations collectively suggests that a trial-to-trial interaction occurs predominantly in the feedforward direction.

To confirm that the observed similarities reflected stimulus-related rather than response-related activity, we tested the temporal consistency of vectors [67] computed at different time points (S13A and S13B Fig). In both areas, the defined during the transient interaction time window (0–100 ms after air-puff onset) showed low consistency with s derived from later time windows, indicating that early patterns uniquely encoded stimulus identity rather than the animals’ behavioral response that emerges closer to the choice period. We next computed similarity values based on a late (1–1.1 s), corresponding to the late delay period. These similarity values were significantly lower than those based on the early in both areas, and this effect was specific to the feedforward direction (S13C–S13F Fig; total FF; S1_cd versus S1_{late cd}: p = 9.4 × 10⁻⁵, z = 3.9, M2_cd versus M2_{late cd}: p = 1.0 × 10⁻⁷, z = 5.3, total FB; S1_cd versus S1_{late cd}: p = 0.16, z = 1.4, M2_cd versus M2_{late cd}: p = 0.15, z = 1.4, noise FF; S1_cd versus S1_{late cd}: p = 0.0004, z = 3.6, M2_cd versus M2_{late cd}: p = 1.1 × 10⁻⁷, z = 5.3, noise FB; S1_cd versus S1_{late cd}: p = 0.88, z = 0.1, M2_cd versus M2_{late cd}: p = 0.68, z = 0.4 Wilcoxon signed-rank tests, n = 50 sessions).

These results confirm that the early feedforward interactions between S1 and M2 contain transient sensory coding patterns relevant for task performance but do not persist into later, preparatory, or action-related periods.

A short delay of S1–M2 interactions suggests a potentially direct, feedforward coupling

To further probe the directionality of S1 and M2 interactions, we examined the temporal delay of interacting signals. Averaging CCV across sessions as a function of delay revealed a feedforward bias, with stronger coupling at positive delays (Fig 5A–5C, top). In many sessions, CCV peaked at a short and positive delay, indicating rapid transmission of sensory information from S1 to M2, compatible with a potentially direct feedforward pathway (Fig 5C, top). Correspondingly, the magnitude of cosine similarity in both regions also exhibited a feedforward bias in most of the sessions (Fig 5A–5C, middle and bottom). The feedforward bias of the CCV peak was especially prominent in high-selectivity sessions (Fig 5D, top, median peak delay: 0 ms versus +5 ms for low- and high-selectivity sessions, respectively, noise: p = 0.018, t = 2.4, total: p = 0.046, t = 2.0, df = 48, linear mixed-effects using each sessions selectivity as a fixed-effect and subjects’ identity as a random effect).

Considering total correlations, similarities showed a more dynamic temporal profile (S14A Fig); both S1 and M2 similarities showed significantly different peak time between early (0–50 ms) and late (50–100 ms) post-stimulus time windows (S14B and S14C Fig, S1: p = 0.018, z = 2.4, M2: p = 0.002, z = 3.1 CCV: p = 0.875, z = 0.16 Wilcoxon signed-rank test). The peak of M2 similarity prominently shifted from a feedforward delay in the early window (median = +5 ms) to a feedback delay in the late window (−15 ms). This temporal shift was invariant to the choice of window around stimulus onset and reproduced when a was defined from 100 to 200 ms after stimulus onset (early window: +7.5 ms; late window: −20 ms, p = 0.001, z = 3.3, Wilcoxon signed-rank test), indicating that the similarity profile emerges due to stimulus selectivity and not merely due to the choice of window. Overall, these results reveal rapid, stimulus-selective feedforward communication from S1 to M2 within 100 ms of the air puff onset. A weaker, delayed feedback communication reflected M2 stimulus selectivity patterns but likely arose through indirect pathways. Notably, a consistent feedforward bias was also observed during inter-trial intervals (S14D–S14G Fig) and in CCVs of second to fourth canonical dimensions, even though these dimensions did not resemble stimulus coding patterns (Figs 6A–6C and S15A–S15C). This suggests that, in parallel to stimulus identity, additional sensorimotor information is likely transferred from S1 to M2 and may contribute to shaping behavioral output.

We additionally tested whether the delay-period interactions potentially align with a choice subspace of long RT sessions. For this analysis, we included long RT trials in which licking occurred within an extended response window (1.2–5 s from air puff) and excluded trials with premature responses (0–1.2 s from air puff). A choice coding subspace was then computed in each area using neural activity during the delay period (0.25–1 s after air puff onset). We derived one choice coding vector for each consecutive 250 ms window (4 total analysis windows), analogous to the stimulus coding direction cd. We then repeated the CCA analyses and computed the similarity between the communication subspace and the choice coding subspace at each 250-ms time window and averaged them across all 4 windows and all sessions. Across all the first 4 canonical dimensions, similarities were flat, near-zero for all delays and for both total and noise correlations (S16 Fig), indicating that the delay-period communication subspace did not align with the choice-coding subspace.

To further consolidate these results, we repeated the CCA with 1-ms binning to more precisely localize the peak of CCVs (S17A–S17C Fig). Across sessions, peak CCV delay systematically decreased as the fraction of selective S1 neurons increased (S17D Fig, median = +3 ms). CCV peak magnitude dropped across canonical dimensions (dim) regardless of S1 selectivity (S1sel) (S17E, S17F, S17H and S17I Fig, top; p_dim = 7.1 × 10⁻¹⁵, t_dim = −8.1, p_S1sel = 0.29, t_S1sel = 1.07, df = 392, linear mixed-effects, 50 sessions from 18 mice), the absolute peak delay depended on both dimension and selectivity (S17G and S17J Fig, top, p_dim = 0.0013, t_dim = −3.2, p_S1sel = 0.0003, t_S1sel = −3.6, p_{dim × S1sel} = 0.043, t_{dim × S1sel} = 2.03, df = 392, linear mixed-effects). High-selectivity sessions had a median absolute peak delay of 4 ms in the first dimension, indicating a rapid feedforward [73,74] coupling from S1 to M2 (S17G Fig, top), whereas low-selectivity sessions showed longer delays overall (S17J Fig, top). To ensure that these delays did not arise from slow or indirect shared dynamics [75], we computed a jittered-CCV [43] and subtracted it to obtain a jitter-corrected CCV. This revealed that the first canonical dimension retained the short median absolute peak delay (4 ms) after correction (S17G Fig, bottom), with no significant interaction between correction and selectivity (p_correction = 0.50, t_correction = −0.68, p_{correction × S1sel} = 0.71, t_{correction × S1sel} = 0.37, df = 392, linear mixed effects).

These results provide convergent evidence that sensory information is transmitted from selective S1 populations to M2 via a fast, potentially direct feedforward pathway, and that this short-delay coupling depends on the degree of stimulus selectivity in S1.

We performed control and alternative analyses to test the reliability of peak delay estimations across a range of parameters and conditions. First, the average CCV maps in Fig 3 suggest that the positive delays are not uniformly produced across all conditions by the CCA analysis itself; in No-go trials (Fig 3C and 3I), the interaction peak appears at zero delay, despite being analyzed with the identical CCA procedures. We, therefore, considered the peak delays of Go and No-go trials separately (S18A–S18H Fig). This analysis clearly shows that no-go trials exhibited delay distributions centered at zero (median = 0 ms, p = 0.17, z = 1.37, Wilcoxon signed rank test, n = 50 sessions), whereas Go trials showed a significantly positive delay (median = +5 ms, p = 0.0006, z = 3.43, n = 50 sessions). Thus, the analysis framework is capable of distinguishing near-synchronous interactions from delayed interactions, arguing against the possibility that the observed positive delay in Go trials arises trivially from temporal autocorrelation, binning, or shared-input structure alone.

Furthermore, we simulated a simple inter-areal interaction to validate the CCA peak delay estimator. We modeled the spiking activity of two neuronal populations with a simple Poisson process (70 neurons each, 150 trials), where a delayed coupling between the two populations was increased 50-fold from a negligible baseline during a predefined 50-ms interaction period (green shading, delay = 10 ms, S18I and S18J Fig). We then applied the identical time-shifted CCA procedure used for the experimental data to these simulated population activities to test the reliability of our peak delay estimation (S18K Fig). Across all tested conditions, the CCV peak closely matched the imposed ground-truth delay, demonstrating that the estimator reliably recovers the direction and magnitude of inter-population temporal offsets rather than producing spurious short positive delays. These results support the interpretation that the positive delays observed in the experimental data reflect genuine delayed interactions captured by the analysis framework.

We additionally tested a range of parameters to check the robustness of the estimated delays per session, using CCA analysis. These include all possible combinations of; bin size (1, 5, or 10 ms), window length (10, 30, or 50 ms), step size (10, 30, or 50), and jitter window (10, 30, or 100 ms). We first applied CCA using all combinations of parameters on ground truth data (Poisson simulation for 300 trials) with +4 ms delay. Peak delay (circles in S19A Fig) showed perfect consistency across parameters (+4 ms) with a bin size of 1 ms. As expected, with a bin size of 5 ms, the peak moved to +5 ms and with a bin size of 10 ms, it was either at 0 or +10 ms (S19A Fig, middle and bottom). Next, we extended the analysis over parameters to an example session (0–100 ms from air puff and ± 50 ms delays) and found similar consistency of the peak delay across most of the parameter combinations (S19B Fig). A window size of 30 ms, together with finer steps, produced a stable delay estimation in the example session. Interestingly, jitter window did not affect the stability of the delay estimation, indicating that the relevant spike timing alignment between S1 and M2 underlying the communication subspace occurs at timescales substantially shorter than jitter windows tested (<10 ms). We applied time-shifted CCA across all sessions and parameter combinations and obtained uncorrected and jitter-corrected CCV peak delays (S20A and S20B Fig). This extensive analysis revealed the stability of the peak CCV timing across combinations of parameters. For a 1 ms bin-size, the distributions of peak delays across parameters remained significantly shifted toward positive values within ±10 ms range in both uncorrected and jitter-corrected analyses (S20C Fig, top row; uncorrected average across sessions = +1.7 ms, p = 0.004; jitter-corrected average = +1.1 ms, p = 0.021; n = 50 sessions, one-sided mouse-level sign-flip permutation tests; see Methods). With a larger 5-ms bin-size, the distributions similarly showed a positive bias, although the estimates became discretized at 0 or +5 ms due to the coarser temporal resolution (S20C Fig bottom row). Furthermore, this positive-delay bias in the distributions was more prominent in high-selectivity sessions compared to low-selectivity sessions (S20D Fig; high-selectivity: uncorrected average = +3.1 ms, p = 0; jitter-corrected average = +2.0 ms, p = 0.010; n = 25 sessions; low-selectivity: uncorrected average = +0.26 ms, p = 0.38; jitter-corrected average = +0.25 ms, p = 0.35; n = 25 sessions, one-sided mouse-level sign-flip permutation tests), consistent with the main findings. These results show that the central conclusion—a small positive coupling latency consistent with rapid feedforward communication from S1 to M2—is robust across a broad range of reasonable analysis parameters.

Furthermore, we applied two alternative directed measures of inter-areal interactions to consolidate these results: (1) pairwise cross-correlation analysis and (2) generalized linear modeling approach (GLM).

First, we investigated cross-correlations between the spiking activity of the two areas. We selected spikes within 100 ms from the air puff. Plotting the cross-correlogram (CCG) in an example session revealed that a large proportion of pairs had a small positive delay (S1 unit leading M2 unit) compared to their shuffled CCGs (where the CCG matrix elements were randomly resampled across delays and units) (S21A Fig). Jitter-corrected CCGs [75] (with 50 ms jittering window) showed a similar, high proportion of short-positive delays compared to the shuffled data (S21B Fig). Pooling all the pairs across sessions, peak delay histograms revealed a prominent number of pairs with a short positive peak delay in both raw and jitter-corrected CCGs across all sessions compared to the peak delay of the shuffled data (S21C Fig; raw average peak delay = +1.5 ms, p = 0.002; shuffled average peak delay = −0.07 ms, p = 0.8; S21D Fig; jitter-corrected average peak delay = +1.0 ms, p = 0.003; shuffled average peak delay = +0.05 ms, p = 0.1; n = 25 high-selectivity sessions, one-sided mouse-level sign-flip permutation tests on data within ±10 ms). The distribution of the peaks closely matches those of CCA analysis (S17D Fig) and reinforces the idea that the two areas are coupled in a predominantly feedforward fashion compatible with a direct pathway. To further illustrate the direction of coupling between strongly connected pairs we filtered all pairs across sessions based on the strength of correlation (peak CCG > 3.29 × standard deviation). The peak delay of most well-connected pairs was at a short positive delay (S21E Fig: average raw peak delay = +2.2 ms, p = 0.002; shuffled average peak delay = −0.14, p = 0.84; S21F Fig: average jitter-corrected peak delay = +1.8 ms, p = 0.002; shuffled = −0.009, p = 0.47; n = 25 sessions, one-sided mouse-level sign-flip permutation tests on data within ±10 ms).

GLM has been widely used as an alternative technique to analyze inter-areal coupling [40]. We therefore applied the same core principle to test interaction directionality in our dataset. We first tested the reliability of delay estimation using a ground truth simulation based on the previously described Poisson model. We ran the simulation for three predefined delays (−4, 0, and +4 ms). We then applied standard GLM to compute average weights of source neurons onto each target neuron across a range of delays (−20 to +20 ms). As shown in S22A Fig, the peak of average weights onto target neurons reliably recovered all three ground-truth coupling delays simultaneously.

Next, we applied this approach to our experimental data and revealed positive peak interaction delays across target (M2) neurons (S22B Fig). We averaged the regression weights across source neurons onto each target neuron for each delay as a proxy of population-level interactions at each delay. The median absolute peak delay across sessions was 4 ms (S22C Fig), matching the CCA-based delay estimator. Notably, the peak delay of individual target neurons was highly dependent on the session’s S1 selectivity (p = 8.3 × 10⁻⁹, t = 5.8, linear mixed-effects, n = 1,953 neurons with at least one non-zero weight across delays were included), with positive peak delays preferentially observed in high-selectivity sessions (S22D Fig), consistent with the CCA analysis (Fig 5D). Importantly, neurons with zero weights in the lasso model were not uniformly distributed across sessions. Instead, the fraction of zero-weight neurons per session was strongly negatively correlated with the session’s S1 selectivity (r = −0.688, p = 3.5 × 10⁻⁸, n = 50 sessions, Pearson’s correlation), and most zero-weight neurons (~86%) originated from low-selectivity sessions.

GLM is commonly applied to neuronal data using Poisson distributions to account for the statistical properties of the data. Using the Poisson-GLM approach, we found results consistent with those from the standard GLM described above. The median absolute peak delay across sessions was 4 ms, and the positive peak delay of each target neuron was significantly correlated with session’s S1 selectivity (p = 6.6 × 10⁻⁶, t = 4.5, linear mixed-effects, n = 2,917 neurons with non-zero weights were included; 111 neurons with zero weights were excluded, all from low-selectivity sessions), consistent with the results obtained using the standard linear model (LM) (S22E Fig). Therefore, our estimation of the coupling delay between S1 and M2 is reproducible via alternative analyses such as GLM-based coupling.

Finally, we considered the potential effect of unit-isolation quality on our delay estimations. Using ~80% of high-quality units with a low ISI violation rate (Methods, S23A Fig), we found consistent results regarding the short positive delay of the CCV across sessions (S23B–S23E Fig, median = +3). We consistently found significant effect of S1 selectivity (S1sel) and dimension (dim) on the absolute peak delay using either all the units (S17G and S17J Fig) or high-quality fraction of the units (S23F–S23H Fig, p_dim = 0.007, t_dim = −2.7, p_S1sel = 0.006, t_S1sel = −2.8, p_{dim × S1sel} = 0.38, t_{dim × S1sel} = 0.88, df = 392, linear mixed-effects). Consistently, we found a small absolute peak delay (median = 5 ms) using the first canonical dimension in the high-selectivity sessions even after jitter correction. These analyses demonstrate that the observed coupling structure and short positive delays are robust to stricter unit-quality criteria and are not attributable to multiunit contamination or spike-sorting artifacts.

These complementary analyses show that CCA-based estimation of the coupling delay between S1 and M2 is reproducible via alternative analyses such as GLM-based coupling and pairwise cross-correlations. Overall, we found delays consistent with the range (5–10 ms) of reported latencies between mouse sensory and motor cortices [73,74,76]. However, establishing these short delays as direct interactions by ruling out potential common input requires more targeted experiments.

Altogether, we revealed a transient, selective, and feedforward-biased communication between S1 and M2 that predicts behavioral performance. This rapid interaction was temporally locked to stimulus onset and contained task-relevant sensory coding information. These findings suggest that S1 sends selective task-relevant information to M2 shortly after stimulus onset to guide behavioral choice (go versus no-go). Thus, our findings uncover a robust feedforward functional connectivity motif linking somatosensory and premotor cortices that underlies successful sensorimotor transformation.

Ethics statement

Experimental procedures were all approved by the Animal Care and Use Committee (IACUC) at Nanyang Technological University (AUP No. A19062 and A24054) and adhered to Singapore’s National Advisory Committee for Laboratory Animal Research (NACLAR) guidelines.

Animals

Adult, male and female mice (Ai32, Jackson Lab stock 012569) were used for this study. Ai32 mice were crossed with PV-Cre mice (Jackson Laboratory stock 017320), and the offspring (PV-ChR2 mice) were used for optogenetic experiments. Animals were housed on a dark/light cycle of 12 h/12 h (light on from 7:00 to 19:00) at a temperature ~21 °C and humidity of 62%. Experiments were carried out in the light cycle between 9:00 and 18:00.

Head-bar implantation and behavioral training

Mice were initially anesthetized with 3% isoflurane in oxygen in a chamber and continued with 1.5%–2.0% on a stereotaxic frame (RWD Life Sciences). Animals received 0.05 mg/kg buprenorphine subcutaneously before surgery and 5 mg/kg meloxicam after surgery. A custom-made head-bar was implanted behind and parallel to the lambdoid sutures on the skull using veterinary tissue glue (3M Vetbond) and black dental cement (Contemporary Ortho-Jet, Black, Lang Dental). Mice were returned to the home cage with free access to food and water and monitored for a week to recover from the surgery. After the recovery period, mice were water-restricted (1 ml per day) for 2 days and were habituated for 2–4 days to head-fixation before training in a delayed go/no-go task controlled by the Presentation software (Neurobehavioral Systems). Each trial began with a 200-ms air puff delivered onto the distal half of the long whiskers on either side of the head. An air puff was followed by a 1-s delay, and white noise was played via a speaker for a maximum of 2 s to cue the response window. We built an IR beam-breaker circuit using separate emitting and detecting LEDs. The LEDs were carefully positioned to register any licking behavior. Beam interruption by tongue was registered via a National Instruments board (USB-6001) and fed to Presentation software to control task events. In go trials, mice were required to lick during the auditory cue to receive a 2.3 μL drop of water via a metal lick port. Miss trials, defined by the absence of a lick during the response window, resulted in neither reward nor punishment. A lick response in no-go trials resulted in a 6-s time-out, added to the inter-trial interval (5.1–7.0 s). Mice typically performed around 300 trials. In most recording sessions, a maximum of 300 trials was predetermined. As a result, most of the sessions (76%) contained 300 trials and 90% contained between 250 and 400 trials.

Craniotomy and electrophysiology

After reaching a stable daily performance, mice were anesthetized with 3% isoflurane in the chamber, followed by 1.5%–2.0% isoflurane on the stereotaxic frame. Two craniotomies (diameters up to ~1 mm) were drilled over the right hemisphere area corresponding to S1 (centered at AP: −1.7 to −2 from bregma and ML: 3.0 to 4.0 from midline) and M2 (centered at AP: 1.3 to 2.2 mm from bregma and 0.5 to 1 mm from midline). Craniotomy sites were sealed with a silicone elastomer (Kwik-Cast, World Precision Instruments). Electrophysiological recording began at least 4 hours after craniotomy to allow mice to recover from anesthesia. Two silicon probes (ASSY-37 H8b, Acute 32-channel H8b probe, 1 shank per recording area, 30 µm site spacing, 9 mm length, Cambridge NeuroTech) were coated on the back side with DiD stain (Invitrogen) and lowered into S1 and M2 using a pair of motorized micromanipulators (MM-500, RWD Life Science) independently at a depth range of 1,000–1,400 µm. Wideband neural signals were sampled at 30 kHz with RHX data acquisition software and Intan 128-channel acquisition system (M4200, Intan Technologies) during the task. Signals from task events and the IR beam-breaker circuit were recorded along with the electrophysiological data at sync. Kwik-Cast was reapplied after each daily recording session to seal the craniotomies. In subsequent recording sessions, silicon probes were inserted into a fresh site, at least 200 μm away from the previous sites. Recording sessions were performed in subsequent days. Each mouse was recorded for 1–4 sessions, with most animals undergoing 3 sessions of recording. After the final recording session, mice were deeply anesthetized and perfused with 0.1 M PBS, followed by 4% paraformaldehyde dissolved in PBS. Brain samples were removed and post-fixed overnight at 4°C, then transferred to 30% sucrose in PBS for 2–3 days at 4 °C. After embedding in OCT, brains were frozen and cut into 100 μm coronal sections using a Leica CM1950 cryostat. The slices were mounted on microscopic slides with Fluoroshield containing DAPI (Sigma-Aldrich) and were imaged with Axio Scan.Z1 microscope (ZEISS).

Behavioral data analysis

Behavioral data were analyzed using custom-written MATLAB codes. Four trial types were defined based on their stimulus time and outcome; (1) hit trials when left air puff triggered licking response, (2) miss trials when mice failed to respond during the response window, (3) false alarm (FA) when mice responded incorrectly during the white noise presentation following a right air puff and (4) correct rejection (CR) when mice withheld licking response during the white noise presentation following a right air puff. Hit and CR rates were defined as the fraction of correct go and no-go trials, respectively. Behavioral performance (discriminability of the stimuli) was computed as d-prime:

(1)

where is the inverse of the cumulative distribution function of the standard normal distribution.

Stimulus and lick event data were analyzed in MATLAB using custom-written code to identify lick and stimulus onset times. Lick events were detected by thresholding digital signals with debouncing at 0.8 V and then removing subsequent high voltage values up to 50 ms, which ensured sufficient debouncing while staying below the licking cycle across all animals. A first-lick response was defined by setting a 1-s refractory period on lick events and selecting the first lick in each bout. For response alignment, first licks occurring within 0.2–6 s after a left air puff were selected. To isolate the epochs of poor engagement, a local hit rate was computed by fitting a local regression function (smooth function in MATLAB with lowess method and a span of 10 trials) on zero and one hit response vectors. A threshold of 0.2 on the local hit rate effectively isolated epochs of poor engagement.

Electrophysiology data analysis

Neural data for each recording probe were separately analyzed for spike isolation using the HerdingSpikes algorithm [112] as part of SpikeInterface [113] in Python, and manually curated with Phy (https://github.com/cortex-lab/phy) as described previously [109]. Broadband data were filtered at 0.3–3 kHz and re-referenced before spike sorting. Spike clusters were manually curated to exclude noise and further isolate spikes based on features provided by Phy’s graphical user interface. These features included amplitude, ISI distributions, correlograms, firing rates, individual and average spike waveforms, and cluster localizations on the shank. Whenever required, clusters were separated from each other based on amplitude drift patterns and ISI distributions, suggesting multi-unit activity. Two or more clusters were merged if the spike localization, correlograms, and amplitude drift patterns suggested that they belonged to the same neuron. All valid spike clusters were included for further analyses regardless of single-unit isolation quality. To test the potential effect of single-unit isolation quality on the results, we computed an ISI violation rate defined as the fraction of ISIs smaller than 2 ms among all ISIs smaller than 200 ms. By thresholding the ISI violation rate at 0.01 in M2 and 0.04 in S1, we included ~80% of units with low ISI violation rates in each session. Average firing rates for the example sessions were calculated using 1-ms binning and 5-ms smoothing (moving average). Baseline firing rate (0–15 ms from air puff) was subtracted from each neuron’s average firing rate.

Optogenetic suppression of S1

PV-ChR2 mice, in which parvalbumin-positive inhibitory neurons express channelrhodopsin-2, were used for optogenetic suppression of the S1 using a fiber-optic-coupled 473 nm blue laser (Shanghai Laser & Optics Century Co.). The open end of the fiber optic was fitted with a zirconia ferrule, and an aluminum foil sleeve was custom fitted onto the ferrule to minimize exposure to scattered laser. During head-bar implantation surgery, the skull above whisker S1 on both hemispheres was thinned and a layer of UV-curing adhesive (Norland Optical) was applied to the area after a thin layer of 3M’s Vetbond to gain optical access into the cortex [12]. To minimize scattered laser light reaching the animals’ eyes, the exposed skull was surrounded by a high wall (4–6 mm high) made of dark acrylic cement. Mice were trained on the go/no-go task previously described. After they reached 70% correct rate for 3 consecutive days, a control experiment was performed by shining a spot of laser (1.5 mm in diameter and a power density of 6.4 mWmm⁻²) on a piece of aluminum foil that was covering the whole skull. The spot was positioned at the right S1 coordinates by attaching the fiber to a motorized micromanipulator. Laser was applied continuously from 100 ms before to 100 ms after the air puff onset in every one out of one-third of the go and one-third of the no/go trials. An additional 100 ms linear ramp-down to zero was applied to prevent the rebound activity in excitatory neurons. Trials were pseudo-randomly interleaved in a manner that only allowed 2 subsequent laser trials with at least a 5.3-s interval between subsequent laser onsets in such trials. Following the control experiment, the laser was delivered to the right S1 (centered at −1.7 mm from bregma and 3.5 mm from the midline) and then to the left S1 in subsequent sessions using the same mice.

Selectivity index

To compute the selectivity index (SI), peri-stimulus spikes of each neuron were binned using 100-ms intervals. A Receiver Operating Characteristic (ROC) analysis was performed for multiple 100-ms windows relative to the time of air puff using spike counts of go trials as “positive signal” and those of no-go trials as “negative signal” distributions. By calculating true positive rates and false positive rates at multiple thresholds, an ROC curve was constructed. Area under ROC curve (AUC) was computed with the trapezoid technique ( in MATLAB). AUC represents the difference in spiking activity between left and right stimuli, considering trial-to-trial variability in firing rate. AUC was scaled between −1 and 1 to represent SI:

(2)

To test statistical significance, we used a resampling technique. Two shuffled spike count vectors were constructed by randomly resampling from the pool of spike counts of all trials. SI was recalculated based on these resampled data and the process was repeated 250 times. The resulting SIs formed a null distribution. A p-value was computed by thresholding the null distribution at the value of the original SI. Sessions were divided into two groups of selectivity based on the significance of SI (p < 0.05, without correcting for multiple comparisons) in the first post-stimulus time window (0–100 ms). Sessions were classified as “high-selectivity” if they were above the median for the fraction of selective S1 neurons (median = 29%). This measure was highly correlated with the average selectivity of all S1 neurons in each session (r = 0.95, p = 2.06 × 10⁻²⁶, Pearson’s correlation, n = 50 sessions) indicating that the fraction of selective neurons is a reliable proxy for the overall sensory selectivity of the S1 population. Most downstream statistical analyses use population selectivity as a continuous variable rather than a thresholded categorical grouping (e.g., linear mixed modeling relating S1–M2 coupling to behavioral performance).

Canonical correlation analysis

Population activity of neurons in each area evolves in a multi-dimensional space where each neuron represents one dimension. Canonical correlation analysis (CCA) relates the spiking activity of two neuronal populations by finding latent dimensions (canonical pairs) representing shared patterns of activity between two regions [39]. The first canonical pair represents the activity along a one-dimensional subspace that captures maximum correlation between the two areas and subsequent canonical pairs capture weaker correlations along independent subspaces. The objective of CCA is to find axes and along which the activity is maximally correlated between the two regions:

(3)

where is and is matrices of activity where denotes samples. and denote the number of dimensions (e.g., neurons) in and , respectively. and are 1 vectors where is a minimum of and . These vectors can be thought of as dimensions or axes in the neuronal activity space denoted as , along which the two areas interact. Canonical correlation value (0 ≤ CCV ≤ 1) captures the magnitude of correlation between and . In practice, CCA returns matrices and where columns of each matrix are independent axes we term “canonical dimensions”. These dimensions are sorted by the magnitude of correlation along each axis. For all analyses, we used the first canonical dimension unless explicitly stated.

We used a time-shifted CCA [43], where for each time window in an equal but shifted time window of was used as input to CCA. In each CCA iteration, we computed up to 4 canonical dimensions. Data was prepared for CCA analysis by isolating spike times of each unit in each trial before and after the first lick or stimulus. Event-aligned spike trains were then binned in each trial using 5-ms intervals. For each iteration of CCA, we used a 50-ms time window consisting of 10 spike count bins, moving with 10-ms steps. Because the CCA algorithm requires the inversion of the covariance matrices [114], neurons with no spikes within the analysis window across trials should be excluded. After this filtering, most sessions (~70%) contained at least 50 neurons in each area (S3H Fig). The activity of each neuron was then z-scored for all CCA analyses. To isolate moment-to-moment co-fluctuations (i.e., “noise correlations”), trial-averaged activity of each neuron was subtracted from the z-scored activity of each trial. For each time window, a was chosen at a delay of , where is the first sample in the 50-ms chunks of the data. We moved in either direction, one sample (5 ms bin) at a time, to obtain shifted iterations of the data. We assigned S1 activity to and M2 activity to . Therefore, the sign of delay determines the direction of interactions; a feedforward interaction (S1 → M2) should occur in positive delays and a feedback interaction (M2 → S1) in negative delays. To evaluate temporal dynamics of interactions in each session, CCVs (Fig 2C–2H) were averaged across 5 delays This yielded a CCV vector showing peri-stimulus dynamics of CCV for each session. To test significance of the fluctuations relative to the air puff onset, a 95% confidence interval for the mean of baseline CCV across sessions (CCV CI) was computed by first taking CCV from −0.25 s to 0 relative to the air puff onset, yielding a matrix. The CI was computed using function in MATLAB (R2019b) with the following parameters; , , and as well as the input data matrix (. The output was a matrix for upper and lower confidence bounds for which the mean values were used to assess deviations from baseline. For response-aligned data, the baseline CCV was taken from −3 s to −2 s relative to the first lick onset. Time-averaged values of CI crossing points in each side of CI were compared against the time-averaged values of the baseline CCVs to verify significant deviations from the baseline using Wilcoxon signed-rank tests. A shuffled CCV was computed for each session using the same procedure as the original CCV, except using trial-shuffled S1 activity to remove the trial-to-trial correspondence between S1 and M2 activity while preserving the temporal structure of the activity in each area. To obtain a higher temporal resolution, spike trains were binned within 1-ms intervals, resulting in 50 samples for each iteration of CCA. This 50-ms window was moved without overlapping (with a step of 50 ms) and each delay step was 1 ms. Corrected CCV maps were computed by subtracting jittered-CCV to reveal potential direct interactions [43]. Jittered-CCV was computed by randomly resampling spikes within 25-ms windows across trials [115].

Geometric mean firing rate

The geometric mean firing rate (gmFR) of each area was computed by averaging across trials and neurons within each 50-ms window, in parallel with CCA analysis for the same period of activity in both regions:

(4)

where and are grand average spike counts of and , respectively in a 50-ms window. gmFR was also calculated for delayed iterations, yielding the same temporal map as with the CCV. In parallel with CCV, to evaluate the temporal dynamics of firing rate in each session, gmFR was averaged across 5 delays .

Coding direction

Coding direction () is an axis that maximally separates evoked activity of left and right air puffs and can be computed by fitting a LM (i.e., a linear classifier). Due to the high dimensionality of the population activity, a limited number of trials and often high trial-to-trial variability in response, linear models such as linear discriminant analysis or support vector machines can overfit the data and distort the actual coding direction. A simpler LM is advisable in this situation. Therefore, we computed coding direction () based on a procedure frequently used in the literature [67,68]. In each area, average spike counts within 100 ms from left air puff onset were subtracted from those from corresponding spike counts evoked by right air puffs:

(5)

Each element of the vector was the trial-average spike count for each neuron in the population. To evaluate the consistency of over time, we computed a for multiple 50-ms windows relative to the stimulus and then constructed a correlation matrix using Pearson’s correlation coefficients. High correlation between in each window with a subsequent window indicates that the same pattern of activity is sustained, and a low correlation indicates a dissimilar activity pattern.

Trial-level coupling

To quantify trial-level coupling, we applied CCA analysis on residual activity (for noise correlations) using 2-ms bins and a 100-ms sliding window (10-ms steps), which yielded 50 temporal samples for each trial. All trials within each session were used to obtain a reliable estimate of the first pair of canonical weights ( and ) and projected all trials along the first canonical dimension ( and ). For each window and delay combination, a Pearson’s correlation coefficient was calculated between the projected population activity of S1 and M2 in each trial. These correlation coefficients represent the within-trial coupling strength between the two areas. Within-trial coupling strength was estimated as the average of these correlation coefficients within 100 ms of the air puff in either the feedforward (FF: 2–20 ms) or feedback (FB: −20 to −2 ms) direction.

As a complementary analysis, we repeated the trial-level coupling using 10-fold cross-validation. For each fold, the first pair of canonical weights was estimated from 90% of the trials (training trials) and applied to the remaining 10% of held-out trials (test trials). Test-trial activity was projected onto the canonical dimensions obtained from the training trials ( and ) and within-trial coupling was calculated from the Pearson’s correlation coefficient between the projected S1 and M2 activity.

Simulation of interareal interactions

To simulate interareal interactions, we modeled spiking activity of two populations with a simple Poisson process. For neuron in trial , the spike probability was defined as

(4)

where is the sum of three components: (1) a baseline firing rate of 5 Hz, (2) a trial-to-trial fluctuation randomly sampled between 0 and 10 Hz for each trial, and (3) a coupling term. For population 1, the coupling term was set to 0, implementing negligible feedback coupling from population 2 to population 1. For population 2, the coupling term was defined as 0.2 times the summed spiking activity of population 1 at a given time point, shifted by the ground-truth delay, thereby implementing feedforward coupling from population 1 to population 2. During a predefined 50-ms interaction period (green shading S18I and S18J Fig), the coupling strength between the two populations was transiently increased 50-fold, where the coupling term for population 2 increased to 10 times the summed spiking activity of population 1, while the coupling term for population 1 remained 0. The ground-truth temporal offset between the populations was systematically varied from −10 ms (population 1 lagging population 2) to +10 ms (population 1 leading population 2). We repeated the simulation for 150 trials and for 70 neurons for each population.

For GLM-based peak delay estimation, we used a variation of the same Poisson model. In each of 300 simulated trials, a 50-ms coupling window was assigned and population activity in each 1-ms bin was defined as the sum of a 5 Hz baseline, 0–10 Hz trial-to-trial fluctuation (shared between populations) and a coupling term. For population 1, the coupling term was set to 0, whereas for population 2, it was set to 70 times the summed spiking activity of population 1 shifted by ground-truth delay. We ran the simulation for three predefined delays (−4, 0, and +4 ms, 20 neurons per population) and then concatenated neurons before applying the lasso regression across 21 delays (−20 to +20 ms).

Pairwise cross-correlation analysis

A cross-correlation analysis was conducted between pairs of higher-quality S1 and M2 neurons (~80% in each area with a lower ISI violation rate as described above). Spikes within 100 ms of the air puff were binned in 1 ms intervals and concatenated across trials. CCG for neuron pair was computed over trials as [75,116];

(4)

where is the number of time bins, is the time delay between the spike train and and are spike counts in a time bin of each trial for S1 and a delayed time bin for M2 neurons. Normalizing denominators consist of a triangular function; which corrects triangular decay over non-overlapping bins and geometric mean firing rate denoted by which corrects for CCG estimation errors caused by variable firing rates across pairs [57]. A jitter window of 50 ms was used to compute a jittered CCG for each pair [75,117]. CCGs were further corrected by subtracting this jittered CCG as follows;

(5)

Jitter correction accounts for slower firing co-fluctuations caused by stimulus- or behavior-locked activity and slow cortical oscillations which may not relate to direct communications. CCG was computed within ±50 ms delays and the peak delay was calculated after applying a 5 ms moving average and limiting CCG to ±40 ms.

GLM-based peak delay estimations

To apply GLM to our inter-areal, we used lasso function in MATLAB, we systematically quantified neural coupling between S1 and M2 after z-scoring the activity of each neuron. Ten candidate values of the LASSO regularization parameter (λ) were evaluated, and 3-fold cross-validation was used to assess model performance across these values, with the λ yielding the lowest cross-validated deviance selected for the final model. The maximum number of optimization iterations was set to 10⁵ for lasso. The resulting model was used to predict each M2 target neuron from time-shifted activity of S1 source neurons. Each M2 neuron’s activity (0–50 ms of air puff, 25 bins of 2 ms × number of trials in each session) was regressed by its own history (−20 to −2 ms, thus 10 self-history predictors) and time-shifted activity of all S1 neurons (−20 to +20 ms, thus 21 × N_S1 predictors), across all trials within each session. We averaged the regression weights across source (predictors) onto each target neuron for each delay as a proxy of population-level interactions. Neurons with zero weight across all delays were excluded from the subsequent statistical analysis testing the association between each target neuron’s peak delay and S1 selectivity. To ensure generalization of our results to the original spiking activity with its unique statistical properties, we repeated the analysis using lassoglm in MATLAB with a Poisson distribution (Poisson GLM, 10 λ with 3-fold cross-validation and 10³ as maximum iterations). Data parameters (e.g., temporal resolution was 2 ms) were identical to those chosen for the standard GLM described above.

Statistics

All statistical analyses were performed using MATLAB and are described in the text and figure legends. All statistical tests were two-sided unless otherwise stated.

To test if variability in population selectivity is due to variable quality of the recordings, a linear mixed-effects model was fit to neuronal data with subject and session identifiers as random effects and the fraction of selective neurons in each area (), as a fixed effects and baseline firing rate of each neuron as the response variable ().

(8)

Another linear mixed-effects model was constructed to test the behavioral relevance of CCV measures with behavioral d-prime as the response variable. Fixed effects included the , fraction of selective M2 neurons () and average CCV within 100 ms of air puff () with full interaction terms as:

(9)

To examine variables predicting delay of the CCV peak, the following linear mixed-effects model was used with (delay of the peak of CCV), (high or low selectivity), (with or without jittered-CCV subtraction) and as fixed effects:

(10)

To test the effect of trial-level coupling on trial outcome, we applied a binomial mixed-effects model with trial outcome (correct versus incorrect) as the response variable, trial-level coupling as a fixed effect, and sessions and subjects as random effects:

(11)

To examine the possible association between the motion index () and trial outcome, we applied a binomial mixed-effects model with trial outcome as the response variable, trial-level motion index a fixed effect, and sessions and subjects as random effects:

(12)

Potential interaction between trial-level motion index and coupling was examined in a binomial mixed-effects model as follows:

(13)

To examine if the motion index can predict coupling at the trial level, the following linear mixed-effects model was used:

(14)

The association between GLM-based peak delay of individual source neurons and was tested using the following linear mixed-effects model:

(15)

To assess the statistical significance of the positive peaks obtained across CCA parameter settings and across CCG pairs, we used a one-sided sign-flip permutation test at the mouse level. Session values were first averaged to obtain a single value for each session, and the grand average across all sessions was used as the observed statistic. For each of 1,000 permutations, an independent random sign (+1 or −1) was assigned to each mouse and applied to all session-wise average values from that mouse. The grand average of the sign-flipped session values across all sessions was then recalculated to generate a null distribution. The p-value was calculated as the proportion of permuted statistics that were greater than or equal to the observed statistic. This procedure tests the null hypothesis that the distribution of session-level values is symmetric around zero, while preserving the dependence among sessions from the same mouse.

S1 Fig. Behavioral performance of mice in the delayed go/no-go task.

(A–C) Hit rate (A), Correct rejection (CR, B) rate and d-prime (d′, C) of 18 mice in the task (Mean ± SEM). (D) Performance of mice in 50 recording sessions. Color range indicates d-prime (d’) for each session. (E and F) DAPI (blue) staining marks cell nuclei and DiD (pink) stain marks probe tracks in S1 (E) and M2 (F). Scale bars: 1 mm. Overlayed gray outlines depict area outlines (modified from Allen Mouse Brain Atlas, P56 coronal sections, E: image 70 and F: image 34) [118].

https://doi.org/10.1371/journal.pbio.3004095.s001

(TIF)

S2 Fig. Neuronal responses in two example sessions indicate reliable and selective responses in both regions.

(A) Responses of simultaneously recorded neurons in S1 within 100 ms of air puff, sorted by depth of each neuron (rows), averaged over hit (left) and correct rejection (CR, right) trials. (B) Average firing rate of S1 across neurons in hit vs. CR (left) and hit vs. miss (right) trials (n = 73 neurons). (C) Responses of simultaneously recorded neurons in M2 in hit (left) and correct rejection (CR, right) trials. (D) Average firing rate of M2 across neurons in hit vs. CR (left) and hit vs. miss (right) trials (n = 95 neurons). (E and F) Quantification of S1 firing rates for Example session 2 (n = 129 neurons). (G and H) Quantification of M2 firing rates for Example session 2 (n = 114 neurons). Average plots are accompanied by the standard error of the mean (±SEM).

https://doi.org/10.1371/journal.pbio.3004095.s002

(TIF)

S3 Fig. Firing rate and selectivity of neurons across 50 sessions.

(A) Fraction of S1 and M2 neurons with significant initial selectivity (0–0.1 s from air puff onset) in each session (gray lines). The arrow points to the median (arrow) of the data in S1. The data underlying this panel can be found in https://doi.org/10.5281/zenodo.23338822 as S3A.xlsx. (B) Cumulative distribution of baseline firing rate (500 ms before stimulus) of all S1 (left) and M2 (right) neurons grouped based on the median fraction of selective S1 neurons. Population selectivity in neither region predicted the baseline firing rates of individual neurons in the same region (S1: p = 0.84, t = −0.2, df = 4,971, M2: p = 0.056, t = −1.91, df = 5,876, linear mixed-effects). (C) Grand average firing rate (in 100-ms bins) of S1 (left) and M2 (right) neurons across all trials in high selectivity (blue) and low selectivity (red) divisions. Black and gray arrows point to air puff and auditory response cue times, respectively. (D) Firing rate of S1 neurons averaged across all trials for high (left) and low selectivity (right) sessions. (E) Average firing rate of M2 neurons. (F) Left vs. right selectivity of S1 neurons in high (left) and low (right) selectivity divisions. (G) Left vs. right selectivity of S1 neurons in high (left) and low (right) selectivity divisions. Note that in F and G, neurons with no significant selectivity from 0 to 3s after air puff were excluded and time bins with a non-significant selectivity index are set to 0. (H) Distribution of the number of recorded neurons (left) and neurons included in the CCA analysis (right) per session in S1 and M2.

https://doi.org/10.1371/journal.pbio.3004095.s003

(TIF)

S4 Fig. Temporal dynamics of interactions in trial types grouped by outcome and in response-aligned data.

(A–D) Total interactions for each trial type; hit (A, p_up = 1.9 × 10⁻⁵, z_up = 4.3 and p_down = 1.6 × 10⁻⁸, z_down = −5.7, upper and lower CI crossing points respectively, within 1 s from air puff vs. baseline, Wilcoxon signed-rank test), correct rejection (CR, B, p_up = 0.00013, z_up = 3.8, p_down = 4.3 × 10⁻⁸, z_down = −5.5), miss (C, p_up = 0.0002, z_up = 3.7) and false alarm (FA, D) respectively. Geometric mean firing rates (gmFR, mean ± SEM) are plotted on the right y-axis in red. (E–H) Noise interactions in hit (E, p_down = 1.9 × 10⁻⁸, z_down = −5.6), correct rejection (F, p_down = 4.3 × 10⁻⁸, z_down = −5.5), miss (G, p_down = 6.8 × 10⁻⁸, z_down = −5.4), and false alarm (H), respectively. (I, J) Response-aligned data show decreases in total interactions in all sessions (I, p_down = 1.9 × 10⁻⁹, z_down = −6.0) and sessions with long median reaction time (mRT > 1.2 s, I, p_down = 5.4 × 10⁻⁷, z_down = −5.0). Lavender dashed lines represent 95% CI of the baseline CCV and the gray dotted line represents the average trial-shuffled data (one shuffle per session). The shaded area represents the standard error of the mean (SEM). Wilcoxon signed-rank test was used to test the differences between average upper and lower CI crossing points within 1s from the air puff onset vs. baseline (n = 50 sessions).

https://doi.org/10.1371/journal.pbio.3004095.s004

(TIF)

S5 Fig. Distribution of trial-level coupling between S1 and M2 in correct and incorrect trials.

(A) Probability density function (pdf) of trial-level coupling in the feedforward direction (FF, S1 → M2) for high selectivity sessions. Correct trials exhibited significantly stronger coupling than incorrect trials (p = 0.0045, t = 2.84, df = 5,242, linear mixed-effects, 25 sessions from 12 mice). (B) Distribution of trial-level coupling in the feedback direction (FB, M2 → S1) for high selectivity sessions with no significant difference between correct and incorrect trials (p = 0.24, t = 1.2, df = 5,242, linear mixed-effects, 25 sessions from 12 mice). (C) Distribution of trial-level coupling in the feedforward direction for low selectivity sessions with no significant difference between correct and incorrect trials (p = 0.80, t = −0.26, df = 5,065, linear mixed-effects, 25 sessions from 13 mice). (D) Distribution of trial-level coupling in the feedback direction for low selectivity sessions with no significant difference between correct and incorrect trials (FB: p = 0.53, t = 0.62, df = 5,065, 25 sessions from 13 mice.

https://doi.org/10.1371/journal.pbio.3004095.s005

(TIF)

S6 Fig. Lack of association between motion index and the S1 → M2 coupling.

(A) Mean motion index (±SEM) averaged across 25 high-selectivity sessions. The right axis (red) shows the average lick rate. (B) Trial-level measurement of z-scored motion index and S1 → M2 coupling show no significant correlation (r = −0.02, p = 0.145, 5,244 trials from 25 high-selectivity sessions, Pearson’s correlation). (C) Grand average trial-level coupling (top; CCV_total, middle; CCV_noise) and motion index (bottom) across high-selectivity sessions (5,244 trials from 25 sessions).

https://doi.org/10.1371/journal.pbio.3004095.s006

(TIF)

S7 Fig. Communication subspaces of S1 and M2 align with their local stimulus coding subspace in total correlations.

(A–C) Average CCV_total across sessions in each time window and delay (A), showing significant interaction in both directions (B, S1 → M2; p_up = 4.1 × 10⁻⁸, z_up = 5.5 and p_down = 1.7 × 10⁻⁶, z_down = −4.9, C, M2 → S1; p_up = 0.0009, z_up = 3.3, p_down = 3.5 × 10⁻⁵, z_down = −4.1, Wilcoxon signed-rank tests). (D–F) Average magnitude of similarity between and _S1 in each time window across sessions showing feedforward (E, S1 → M2; p_up = 5.9 × 10⁻⁸, z_up = 5.4) and feedback (F, M2 → S1; p_up = 0.001, z_up = 3.3, p_down = 0.07, z_down = −1.8, Wilcoxon signed-rank tests) components. (G–I) Average magnitude of similarity between and _M2 with significant feedforward (H, S1 → M2; p_up = 3.6 × 10⁻⁸, z_up = 5.5, and p_down = 0.002, z_down = −3.1) and feedback components (I, M2 → S1; p_up = 0.0002, z_up = 3.7, p_down = 0.017, z_down = −2.4, Wilcoxon signed-rank tests). Average data are represented as mean ± SEM. White lines indicate delay intervals selected for feedforward (FF) and feedback (FB) components. Lavender dashed lines represent 95% confidence interval of the baseline and gray dotted lines represent the shuffled CCV_total. Wilcoxon signed-rank was used to test the differences between average upper and lower CI crossing points within 200 ms from the air puff onset vs. baseline (n = 50 sessions).

https://doi.org/10.1371/journal.pbio.3004095.s007

(TIF)

S8 Fig. Requirement of early S1 activity for behavioral discriminability and early M2 selectivity.

(A) Schematic illustration of optogenetic suppression of S1 (blue). (B) Top: Behavioral discriminability (d’) in trials without (OFF) and with (ON) laser stimulation during the foil control, right S1 and left S1 sessions. Each color represents an individual mouse (n = 4), and connected diamonds show the within-session laser effect of laser stimulation. Bottom: Within-session change in behavioral discriminability induced by laser stimulation (∆d′ = d′_ON − d′_OFF). Brief suppression of the right S1 produced a significantly larger reduction in d′ than either the foil control or left S1 suppression (repeated-measures ANOVA: p = 0.015, F = 9.2, df = 2, Fisher’s LSD: Control vs. Right S1: p = 0.019, Left S1 vs. Right S1: p = 0.042). The data underlying this figure can be found in https://doi.org/10.5281/zenodo.23338822 as S8B.xlsx. (C) Left: Average time course of M2 selectivity index during trials without (OFF, black) and with (ON, blue) right S1 suppression (mean ± SEM, n = 147 neurons). Right: Selectivity index of individual M2 neurons exhibiting significant selectivity within the first 100 ms after air puff during OFF and ON trials. Brief right S1 suppression significantly reduced early M2 selectivity (linear mixed-effects model: p = 0.0017, t = −3.17, df = 292; 9 sessions from the same 4 mice shown in B). Each gray line represents one neuron. Boxplots represent the median and interquartile range (IQR); whiskers extend to the most extreme data points within 1.5 × IQR of the first and the third quartiles, and notches represent the 95% confidence intervals of the median. The data underlying this figure can be found in https://doi.org/10.5281/zenodo.23338822 as S8C.xlsx

https://doi.org/10.1371/journal.pbio.3004095.s008

(TIF)

S9 Fig. No change in coupling direction and coding-communication alignment in S1 and M2 after controlling for single-unit isolation quality.

(A) Average autocorrelogram of well-isolated single units (blue; ~80% of recorded units) and excluded units with high inter-spike interval (ISI) violation rate (red) in S1 (left) and M2 (right). Shaded areas represent mean ± SEM across units. (B) Analysis of CCV_noise using only well-isolated single units. Top: Average CCV_noise, S1 similarity, and M2 similarity maps across sessions. Middle: Feedforward coupling (left, p_up = 6.8 × 10⁻⁵, z_up = 4.0), S1 similarity (middle, p_up = 2.2 × 10⁻⁵, z_up = 4.2), and M2 similarity (right, p_up = 0.002, z_up = 3.0) remained significant after excluding low-quality units. Bottom: Feedback coupling and similarity measures showed no significant upper confidence interval (CI) crossing points. (C) Same as (B) but CCV_total. Top: Average CCV_total, S1 similarity, and M2 similarity maps across sessions. Middle: Feedforward coupling (left, p_up = 4.5 × 10⁻⁸, z_up = 5.5), S1 similarity (middle, p_up = 2.8 × 10⁻⁸, z_up = 5.5), and M2 similarity (right, p_up = 1.4 × 10⁻⁶, z_up = 4.8) remained significant after excluding low-quality units. Bottom: Feedback coupling (left, p_up = 0.0014, z_up = 3.2), S1 similarity (middle, p_up = 3.0 × 10⁻⁴, z_up = 3.6), and M2 similarity (right, p_up = 1.2 × 10⁻⁴, z_up = 3.8) showed significance. Data are represented as mean ± SEM. White horizontal lines indicate the delay windows used to quantify feedforward (top) and feedback (bottom) components. Lavender dashed lines represent the 95% confidence interval of the baseline, and gray dotted lines represent shuffled CCV. Upper and lower CI crossing points within 200 ms after air-puff onset were compared with baseline using the Wilcoxon signed-rank test (n = 50 sessions).

https://doi.org/10.1371/journal.pbio.3004095.s009

(TIF)

S10 Fig. Split-half cross-validation confirms consistent feedforward coupling and coding-communication alignment in S1 and M2.

(A) Cross-validated average maps of CCV_noise, S1 similarity, and M2 similarity across sessions. Communication and coding subspaces were computed from independent halves of trials within each session. (B) Feedforward CCV_noise (left; p_up = 0.006, z_up = 2.7), cross-validated S1 similarity (middle; p_up = 0.0033, z_up = 2.9), and cross-validated M2 similarity (right; p_up = 0.033, z_up = 2.1), demonstrating consistent feedforward coupling and coding-communication alignment. (C) Feedback CCV_noise (left), cross-validated S1 similarity (middle), and cross-validated M2 similarity (right), showing no significant increase, consistent with the non-cross-validated results. (D) Cross-validated, average maps of CCV_total, S1 similarity, and M2 similarity across sessions, computed as in (A). (E) Feedforward CCV_total (left; p_up = 3.4 × 10⁻⁷, z_up = 5.1), cross-validated S1 similarity (middle, p_up = 5.8 × 10⁻⁶, z_up = 4.5), and cross-validated M2 similarity (right, p_up = 5.8 × 10⁻⁶, z_up = 4.5), confirming consistent feedforward coupling and coding-communication alignment. (F) Feedback CCV_total (left; p_up = 0.063, z_up = 1.9), cross-validated S1 similarity (middle; p_up = 0.03, z_up = 2.2) and cross-validated M2 similarity (right; p_up = 0.01, z_up = 2.5). Data are represented as mean ± SEM. White horizontal lines indicate the delay windows used to quantify feedforward (top) and feedback (bottom) components. Lavender dashed lines represent the 95% confidence interval of the baseline, and gray dotted lines represent shuffled CCV_total. Upper and lower CI crossing points within 200 ms after air-puff onset were compared with baseline using the Wilcoxon signed-rank test (n = 50 sessions).

https://doi.org/10.1371/journal.pbio.3004095.s010

(TIF)

S12 Fig. Consistency of cosine similarity sign between S1 and M2.

(A) Two example sessions showing CCV_noise (1st column), signed S1 similarity (2nd), signed M2 similarity (3rd), and the similarity product (4th) for each time bin. (B) The average similarity product across sessions revealed a strong, transient alignment between coding and communication subspaces selectively in the feedforward (FF) direction, consistent with the results obtained using absolute cosine similarities. Left: Average similarity product across 50 sessions for each time bin. Middle: Feedforward similarity product, showing a significant increase within the first 100 ms after stimulus onset (p_up = 5.9 × 10⁻⁷, z_up = 5.0, upper confidence interval (CI) crossing points compared with baseline, Wilcoxon signed-rank test). Right: Feedback (FB) similarity product, showing no significant change (p_down = 0.18, z_down = −1.3, lower CI crossing points compared with baseline, Wilcoxon signed-rank test).

https://doi.org/10.1371/journal.pbio.3004095.s012

(TIF)

S13 Fig. Low consistency of coding directions over trial time.

(A and B) Consistency (Pearson’s correlations) of calculated in each 100-ms time bin relative to the air puff onset, in S1 and M2, respectively. Arrows point to consistency values relative to the used to compute the similarity (0–0.1s). The late (1–1.1 s) is chosen for comparison in C–F. (C and D) Similarity calculation across sessions for feedforward (FF) components in total (C, S1_cd vs. S1_{late cd}: p = 9.4 × 10⁻⁵, z = 3.9, M2_cd vs. M2_{late cd}: p = 1.0 × 10⁻⁷, z = 5.3) and noise (D, S1_cd vs. S1_{late cd}: p = 0.0004, z = 3.6, M2_cd vs. M2_{late cd}: p = 1.1 × 10⁻⁷, z = 5.3 Wilcoxon signed-rank tests, n = 50 sessions) interactions. (E and F) Similarity calculation across sessions for feedback (FB) components in total (E, S1_cd vs. S1_{late cd}: p = 0.16, z = 1.4, M2_cd vs. M2_{late cd}: p = 0.15, z = 1.4) and noise (F, S1_cd vs. S1_{late cd}: p = 0.88, z = 0.1, M2_cd vs. M2_{late cd}: p = 0.68, z = 0.4, Wilcoxon signed-rank tests, n = 50 sessions) correlations. Dark green lines and circles represent sessions. Light green lines are the median of the data and blue boxes show interquartile ranges. The data underlying C–F can be found in https://doi.org/10.5281/zenodo.23338822 as S13CDEF.xlsx

https://doi.org/10.1371/journal.pbio.3004095.s013

(TIF)

S14 Fig. Total correlations predominantly show a feedforward bias, but a delayed feedback alignment with selectivity emerges in the late activity of both regions.

(A) Baseline-corrected CCV_total profile (top), S1 similarity (middle) and M2 similarity (bottom) averaged across sessions. (B) Early (left) and late (right) CCV_total (top), S1 similarity (middle) and M2 similarity (bottom) averaged (±SEM) across sessions within the first and the second 50 ms, respectively, for each delay. (C) Early (left) vs. late (right) delay of the peak CCV_total (top), S1 similarity (middle), and M2 similarity (bottom) for each session (CCV_total: p = 0.875, z = 0.16, S1: p = 0.018, z = 2.4, M2: p = 0.002, z = 3.1, Wilcoxon signed-rank tests). (D) Total (top) and noise (bottom) interaction maps during inter-trial intervals (up to 200 ms before air puff) reveal a feedforward bias. (E) Average (±SEM) of CCV_total (top) and CCV_noise (bottom) during ITI for each delay across sessions. (F) Delay of the peaks in CCV_total (top) and CCV_noise (bottom) showing feedforward but larger values compared to post-stimulus CCV (C, top). (G) Delay of the peaks as a function of population selectivity for CCV_total (top) and CCV_noise (bottom). Each circle represents a session (n = 50, 18 mice).

https://doi.org/10.1371/journal.pbio.3004095.s014

(TIF)

S15 Fig. Subsequent canonical dimensions display a feedforward bias in total correlations.

(A) Average (±SEM) of CCV_total (top row), S1 similarity (middle row), and M2 similarity (bottom row) for each delay across canonical dimensions (columns). (B) Peak magnitude of CCV_total (top), S1 similarity (middle), and M2 similarity (bottom) across 4 canonical dimensions. (C) Delay of the CCV_total (top), S1 (middle), and M2 (bottom) similarity across canonical dimensions. Each gray line represents one session. Boxplots represent the median (green) and interquartile range (IQR, blue) of the data in each dimension and the whiskers extend to the last data point not considered an outlier. Red plus signs represent outliers (1.5 times IQR below and above the first and third quartiles, respectively). Notches represent 95% confidence intervals of medians. The data underlying B and C can be found in https://doi.org/10.5281/zenodo.23338822 as S15BC.xlsx

https://doi.org/10.1371/journal.pbio.3004095.s015

(TIF)

S16 Fig. Consistently weak delay-period coupling across delay intervals, canonical dimensions and correlation types (total and noise), with no alignment between choice-coding and communication subspaces in either S1 or M2.

Top three rows: CCV_total (top), S1 similarity (middle), and M2 similarity (bottom). Columns 1–4: Average values (±SEM) for the 1st through 4th canonical dimensions. Column 5: Peak magnitude across 4 canonical dimensions. Column 6: Peak delay across canonical dimensions. Gray lines represent individual sessions. Boxplots indicate the median (green line) and interquartile range (IQR; blue box); whiskers extend to the most extreme data points within 1.5 × IQR of the first and third quartiles, red crosses denote outliers, and notches represent the 95% confidence interval of the median. Bottom three rows: Same analyses as above, shown for CCV_noise. The data underlying boxplots (right two columns) can be found in https://doi.org/10.5281/zenodo.23338822 as S16_BoxPlots.xlsx

https://doi.org/10.1371/journal.pbio.3004095.s016

(TIF)

S17 Fig. Higher temporal resolution (1 ms) and jitter correction revealed short-delay coupling in the feedforward direction.

(A) Total canonical correlation values (CCVs) were computed using 1 ms spike time binning in 50 ms windows with a 1 ms step size. (B) Average CCV within 50 ms of air puff shows a clear peak at a positive delay. (C) Distribution of the peak magnitudes of each session’s CCV across delays. (D) Delay of the peaks as a function of S1 population selectivity. (E) Average (±SEM) of raw (top row) and jitter-corrected (bottom row) CCVs as a function of canonical dimensions (columns) across 25 high-selectivity sessions within 50 ms of air puff. (F) Magnitude of peaks in raw (top), and jitter-corrected (bottom) CCVs across dimensions. (G) Absolute delay of the peaks in raw (top) and jitter-corrected (bottom) CCVs. (H) Average of raw (top row) and jitter-corrected (bottom row) CCVs as a function of canonical dimensions (columns) across 25 low-selectivity sessions within 50 ms of air puff. (I) Magnitude of peaks in raw (top), and jitter-corrected (bottom) CCVs in low-selectivity sessions. (J) Absolute delay of the peaks in raw (top) and jitter-corrected (bottom) CCVs in low-selectivity sessions. Boxplots represent the median (green) and interquartile range (IQR, blue) of the data in each dimension and the whiskers extend to the last data point not considered an outlier. Red plus signs represent outliers (1.5 times IQR below and above the first and the third quartiles, respectively). Notches represent 95% confidence intervals of medians. The data underlying F, G, I and J can be found in https://doi.org/10.5281/zenodo.23338822 as S17FGIJ.xlsx

https://doi.org/10.1371/journal.pbio.3004095.s017

(TIF)

S18 Fig. Coupling delay estimation using time-shifted CCA in go vs. no-go trials.

(A-D) Total canonical correlation values (CCV_total) computed using go trials only (5-ms spike bins, 50-ms analysis window, 10-ms step size). (A) Average CCV map across sessions, (B) Average CCV as a function of delay. (C) Peak CCV magnitude versus peak delay of CCV for individual sessions. (D) Peak delay versus the fraction of selective S1 neurons (Frac. Sel. S1), showing a median peak delay of +5 ms (p = 0.0006, z = 3.43, Wilcoxon signed-rank test, n = 50 sessions). (E-H) CCV_total computed using no-go trials only. (E) Average CCV map across sessions. (F) Average CCV as a function of delay (G) Peak CCV magnitude versus peak delay for individual sessions. (H) Peak delay versus the fraction of selective S1 neurons (Frac. Sel. S1), showing no significant delay bias (median peak delay = 0 ms; p = 0.17, z = 1.37, n = 50 sessions). (I, J) Raser plots (black dots) and firing rates (red) of two simulated populations (70 neurons each) in an example trial. Population 1 (I) drove population 2 (J) within a predefined 50-ms coupling window with a ground-truth delay of 10 ms. (K) CCV_noise computed for simulated populations with ground-truth delays ranging from −10 ms (population 1 lagging population 2) to +10 ms (population 1 leading population 2). Time-shifted CCA accurately recovered the imposed coupling delays across the full range of simulated conditions.

https://doi.org/10.1371/journal.pbio.3004095.s018

(TIF)

S19 Fig. Robustness of CCA-based coupling delay estimation across a range of parameters.

(A) Stability of CCA-based peak delay estimation in simulated populations with a ground-truth delay of 4 ms across a range of analysis parameters: spike bin size (b), analysis window size (w), step size (s), and jitter window size (jw), all expressed in milliseconds. The top, middle, and bottom panels show the estimated peak delay (green circles) for parameters combinations using 1-, 5-, and 10-ms spike bin sizes, respectively. (B) Stability of CCA-based peak delay estimation in an example recording session across the same range of analysis parameters. The format is the same as (A).

https://doi.org/10.1371/journal.pbio.3004095.s019

(TIF)

S20 Fig. Stability of CCA-based coupling delay estimation across the parameters for all recording sessions.

(A) Distribution of peak delays for uncorrected CCV_total across the 50 recording sessions over a range of analysis parameters (b, spike bin size; w, analysis window size; s, step size; jw, jitter window size; all in milliseconds) using a 1-ms spike bin size. Green dots represent individual sessions and demonstrate consistent delay estimates across parameter combinations. (B) Distribution of peak delays for jitter-corrected CCV_total across the same parameter combinations using a 1-ms spike bin size, showing that the estimated delays are largely insensitive to the choice of jitter window. (C) Top: Probability density functions (pdfs) of estimated peak delays within ±10 ms for uncorrected (left, average = +1.7 ms, significantly larger than zero, p = 0.004, n = 50) and jitter-corrected (right, average = +1.1 ms; p = 0.021, n = 50; one-sided mouse-level sign-flip permutation tests) CCVs using a 1-ms spike bin size across parameter combinations. Bottom: Corresponding distributions obtained with a 5-ms spike bin size, showing consistent but more discretized peak-delay estimates. (D) Distribution of estimated peak delays in high- and low-selectivity sessions using a 1-ms spike bin size across parameter combinations. Top: Uncorrected CCVs (high-selectivity: average = +3.1 ms, p = 0, n = 25; low-selectivity: average = +0.26 ms, p = 0.38, n = 25). Bottom: Jitter-corrected CCVs (high-selectivity: average = +2.0 ms, p = 0.010, n = 25; low-selectivity: average = +0.25 ms, p = 0.35, n = 25; one-sided mouse-level sign-flip permutation tests). Left and right panels correspond to high- and low-selectivity sessions, respectively.

https://doi.org/10.1371/journal.pbio.3004095.s020

(TIF)

S21 Fig. Pairwise cross-correlation analysis recapitulates the dominance of feedforward coupling between S1 and M2.

(A) Left: Raw (uncorrected), z-scored pairwise cross-correlograms (CCGs; example session) computed within the first 100 ms after air puff. Many unit pairs exhibited peaks at short positive delays (red box, S1 → M2), consistent with fast, feedforward coupling. Right: Time- and unit-shuffled CCGs from the same session. Unit pairs are sorted by their peak delay. (B) Same as (A) but for jitter-corrected CCGs. (C) Distribution of the peak delays of raw CCGs across all unit pairs from high-selectivity sessions showing a peak at short positive delay (average = +1.5 ms, p = 0.002, average shuffled peak delay = −0.07 ms, p = 0.8, n = 25, one-sided mouse-level sign-flip permutation test on data within ±10ms). (D) Same as (C) but for jitter-corrected CCGs (average = +1.0 ms, p = 0.003, average shuffled peak delay = +0.05 ms, p = 0.1, n = 25). (E, F) Raw (E) and jitter-corrected (F) CCGs of unit pairs with the highest peak magnitudes (top 5%). Most pairs exhibited short positive delays consistent with the fast feedforward S1 → M2 coupling revealed by CCA. E: average raw peak delay = +2.2 ms, p = 0.002, shuffled data; −0.14 ms, p = 0.84 and F: average jitter-corrected peak delay = +1.8 ms, p = 0.002, shuffled data; −0.009 ms, p = 0.47, n = 25 high-selectivity sessions, one-sided mouse-level sign-flip permutation tests within ±10 ms).

https://doi.org/10.1371/journal.pbio.3004095.s021

(TIF)

S22 Fig. Generalized linear modeling of S1–M2 interactions confirms the predominance of feedforward coupling identified by CCA.

(A) Linear model (LM) applied to simulated neuronal populations accurately recovered the ground-truth interaction delays of −4 ms (neurons 1–20), 0 ms (neurons 21–40), and 4 ms (neurons 41–60). (B) LM applied to an example recording session within the first 100 ms after air puff onset. Most M2 neurons exhibited peak coupling weights at positive delays (S1 leading M2), consistent with feedforward interactions. (C) Distribution of peak interaction delays across sessions as a function of S1 selectivity recapitulated the peak delay intervals observed with CCA (S17D Fig). (D) Distribution of peak interaction delays across target neurons, estimated by the LM in high- and low-selectivity sessions (p = 8.3 × 10⁻⁹, t = 5.8; linear mixed-effects model with session-level selectivity as a fixed effect; n = 1,953 neurons with at least one non-zero coupling weight across delays). (E) A Poisson generalized linear model (GLM) yielded similar results, with peak interaction delays remaining predominantly positive and varying systematically with S1 selectivity (p = 6.6 × 10⁻⁶, t = 4.5; linear mixed-effects model; n = 2,917 neurons with at least one non-zero coupling weight across delays).

https://doi.org/10.1371/journal.pbio.3004095.s022

(TIF)

S23 Fig. Unit isolation quality does not affect CCA-based delay estimation.

(A) Average autocorrelograms of well-isolated single units (blue; ~80% of recorded units) and excluded units with high inter-spike interval (ISI) violation rates (red) in S1 (top) and M2 (bottom). Shaded areas represent mean ± SEM across units. (B) Average CCV_total map across sessions computed using only well-isolated single units (1-ms spike bins, 50-ms analysis windows, 50-ms step size). (C) Average CCV_total as a function of delay within the first 50 ms after air puff onset, showing a clear peak at positive delays. (D) Distribution of peak CCV_total magnitudes as a function of delay for individual sessions. (E) Peak delay as a function of S1 population selectivity. (F–H) CCV_total results for high-selectivity sessions. (F) Average raw (top) and jitter-corrected (bottom) CCVs across canonical dimensions within the first 50 ms after air-puff onset (mean ± SEM). (G) Peak CCV magnitude across canonical dimensions. (H) Absolute peak delay across canonical dimensions. (I–K) Same analyses as in F–H for low-selectivity sessions. Boxplots represent the median (green line) and interquartile range (IQR; blue box). Whiskers extend to the most extreme data points within 1.5 × IQR of the first and third quartiles, red crosses denote outliers, and notches represent the 95% confidence interval of the median. The data underlying G, H, J and K can be found in https://doi.org/10.5281/zenodo.23338822 as S23GHJK.xlsx

https://doi.org/10.1371/journal.pbio.3004095.s023

(TIF)