Nested spatiotemporal theta–gamma waves organize hierarchical processing across the mouse visual cortex

Domain: physics
Published: 2026-01-30
Editorial outcome: major_revision

Abstract
The brain processes visual information across diverse spatial and temporal scales, but the dynamics that bridge these scales remain poorly understood. Here, we analyze Neuropixels recordings from mouse visual cortex during a visual change-detection task and uncover a nested spatiotemporal structure linking low-frequency θ, high-frequency γ, and neuronal spiking across cortical layers and areas. θ forms broad traveling waves that flexibly switch direction across a stimulus presentation. Following stimulus onset, top-down θ propagates from deep to superficial layers and from higher to lower-order areas, but after stimulus offset, bottom-up θ travels in the reversed direction. These nonstationary θ waves predict behavioral responses to changes in stimulus content. In contrast, γ activity appears as brief, spatially localized packets that sharpen with stimulus onset. Extending conventional cross-frequency coupling analyses, we show that both γ packets and neuronal spiking are phase-locked to the traveling θ waves. This cross-scale coupling exhibits systematic gradients across layers and areas aligned with the visual hierarchy. Together, these findings reveal a spatiotemporal θ–γ neural code that multiplexes bottom-up and top-down information, linking cross-scale dynamics to adaptive hierarchical visual processing.

Introduction
Neural activity unfolds across multiple spatial and temporal scales to support flexible brain function. From microscopic spiking packets1,2 to macroscale traveling waves3,4, diverse dynamical motifs emerge at various levels to enable efficient neural computation amidst high-dimensional, scale-free fluctuations5. Interactions across scales are thought to play a critical role in the flexible, context-dependent processing of complex, multimodal information6. Yet, computational demands shift along the hierarchy of brain regions7, with lower regions performing simple feature extraction and higher regions being involved in cognitive control. Processing across hierarchical levels is integrated by an intricate web of feedforward and feedback neuronal connections between different brain regions and cortical layers8,9. A key challenge, therefore, is to understand how cross-scale dynamics are jointly organized in space and time, and how this spatiotemporal organization supports flexible hierarchical processing in the brain.

Nested oscillations are a leading framework for understanding interactions between dynamics at multiple scales. For example, low-frequency θ oscillations have been associated with synchronizing neuronal activity, conveying top-down information such as expectation10–12, and supporting multimodal integration13. By contrast, high-frequency γ oscillations reflect local processing14,15 or targeted communication16,17, and propagate visual information through the visual hierarchy18 According to the classical θ–γ code, γ bursts occur at specific θ phases to encode different information and provide a temporal scaffold for neural processing19. Recent work has broadened the classical θ–γ framework, extending the dynamical mechanisms underlying working memory20 to describe the functional roles of cross-scale interactions in the cortex2,6,21. θ–γ interactions—and phase–amplitude coupling more generally—have been observed across brain regions, species, and tasks22,23, suggesting nested oscillations serve as a key mechanism for general neural computation. Despite their ubiquity, we still lack a comprehensive account of how nested oscillations are organized across cortical laminae and inter-areal networks. Understanding the role of cross-scale coupling in distributed, hierarchical neural processing requires moving beyond a local, univariate picture.

Recent studies highlight that cortical dynamics can be structured in space as well as time, with traveling waves playing an important role in coordinating distributed neuronal assemblies24–32. In the mouse visual cortex, waves that travel parallel to the cortical surface33,34 have been found in the θ band, propagating in feedback directions, as well as in the γ band, propagating in feedforward directions35. However, the functional significance of propagation across laminae remains unknown. Typical approaches have analyzed traveling properties using trial-averaged measurements of the local field potential (LFP) or current source density36, which can obscure trial-by-trial variability and convert discrete activation events into spurious propagations. Understanding how nested dynamics unfold across both space and time requires a new class of methods capable of capturing the instantaneous, nonstationary, and highly variable nature of neural population activity.

Here, we extend the local θ–γ code into a distributed, wave-based framework. In this framework, large-scale θ waves propagate across cortical laminae and along the visual hierarchy, modulating localized γ wave packets and spiking to shape visual behavior. Using high-density electrophysiological recordings from the mouse visual cortex37, we first characterize the spectral structure of LFPs, identifying hierarchically-organized low-frequency and high-frequency components superimposed on a 1/f background. Based on the dominant spectral peaks, we operationally define θ as 3–10Hz and γ as 30–100Hz. Using wave-based analytical tools, we then show that translaminar and inter-areal θ traveling waves flexibly switch between top-down and bottom-up modes depending on visual stimulation and task demands. At stimulus onset during active task engagement, θ travels from deep to superficial layers and backward along the visual hierarchy, from higher-order to lower-order areas. This onset-related θ propagation is strongest during engagement in a visual task and predicts behavioral responses to changes in a naturalistic stimulus. At stimulus offset, θ reverses direction, traveling from superficial to deep layers and forward along the hierarchy, even during passive viewing. These opposing modes reconcile debate over whether θ supports bottom-up or top-down processing10,38,39, and situates θ coherence as a snapshot of an underlying traveling wave40,41. Finally, we show that θ waves coordinate with task-relevant γ activity and neuronal spiking across layers and areas. γ activity occurs as brief, spatially localized bursts (wave packets) that align to different θ phases across trial periods. Spikes are nested within localized γ packets and large-scale θ waves: spike timing is modulated by γ amplitudes and θ phases, and firing rates covary with θ propagation and task performance. Together, these findings link cross-scale θ–γ wave dynamics to hierarchical visual processing, and suggest that a distributed spatiotemporal θ–γ code supports flexible computation in mouse visual cortex.

Results
We propose that distributed visual computation harnesses spatiotemporal patterns with three key components, illustrated in Fig. 1a: (i) low-frequency, large-scale θ waves that travel vertically (between cortical layers) and horizontally (between cortical regions); (ii) high-frequency, short, and localized γ packets (spatiotemporal bursts) that reflect focused processing or targeted projections; and (iii) cross-scale coupling between θ traveling waves, γ wave packets, and neuronal spiking. Here we investigated these rich dynamics using state-of-the-art high-density, multi-region laminar recordings from the mouse visual cortex, available as the Allen Neuropixels–Visual Behavior dataset37,42. These recordings, illustrated in Fig. 1b, feature fine spatial resolution, with LFP channels spaced vertically by 40 μm, and extensive spatial coverage, with probes in each of six visual cortical regions (VISp, VISl, VISrl, VISal, VISpm, and VISam, ordered by increasing anatomical hierarchy rank9). Temporally, these data comprise high-resolution 1250Hz LFP and spike data from hundreds of neurons over recording sessions lasting ~3 h, with a variety of stimuli, including: a static gray screen (spontaneous), bright or dark flashes (flash_250ms), and natural images paired with a change detection task (Natural_Images). During the task, depicted in Fig. 1b, natural images were flashed for 250ms, followed by 500ms of an intervening gray screen, and mice were trained to lick a reward spout when the image content changed. These stimuli—especially the change-detection task—and the large cohort of 81 mice let us study the mechanisms by which nested spatiotemporal dynamics support top-down and bottom-up processing over distributed brain areas.

We leveraged the dense vertical sampling of the Neuropixels probes, illustrated in Fig. 1b, and a suite of wave-based methods to quantify nested spatiotemporal dynamics in the LFP signal. After identifying the dominant timescales of visual cortical activity, we filtered the raw LFP (shown in Fig. 1c) into two frequency bands: θ (3–10Hz) and γ (30–100Hz). We characterized the θ component as a traveling wave by computing the instantaneous phase, depicted in Fig. 1d, which describes the progression of a time series through its peaks and troughs using an angular value from −π to π. This approach normalizes intrinsic amplitude differences between channels that may arise from variations in laminar electrical response properties43 (unlike methods that quantify laminar activity using current source densities) and is particularly suited to detecting traveling patterns. By examining the spatiotemporal organization of the instantaneous phase, we found that θ activity is organized as large-scale traveling waves that propagate between laminae and across cortical areas. Specifically, we measured the directionality of θ-wave propagation across cortical layers using the negative spatial gradient of θ phase (as illustrated in Fig. 1e). To emphasize the importance of correctly orienting the derivative under different phase conventions, and for conciseness, we refer to this negative phase gradient as the wavenumber [see Eq. (6)]. Then, to examine the spatiotemporal structure of the higher-frequency γ component, we used the instantaneous amplitude, depicted in Fig. 1f. The instantaneous amplitude tracks the envelope of the high-frequency oscillation, revealing the presence of non-stationary, spatially localized γ packets across cortical layers. Finally, we investigated how large-scale traveling θ waves coordinate with localized γ packets and neuronal spiking across cortical layers to form cross-scale spatiotemporal structures, as illustrated in Fig. 1g.

We propose that these nested dynamics form a spatiotemporal θ–γ code, which builds on the classical univariate θ–γ code by accounting for interactions between propagating θ waves and localized γ packets. Figure 1a illustrates how this principle implements two key features of flexible computation: multiplexing and modulation. Consider a hypothetical scenario involving paired γ packets, I and II, that occur near the onset of a visual stimulus. Packet I arises in the deep layers of a lower-order visual area, reflecting top-down information conveyed through feedback connections. Shortly thereafter, packet II emerges in the superficial layers, carrying a potentially independent bottom-up signal. Critically, packets I and II are separated in space and time—spatiotemporally multiplexed—to minimize interference. Packet I is locked to a θ wave that propagates up across the cortical laminae. When the θ wave reaches the superficial layers, its trough coincides with packet II, thereby enhancing—modulating—the bottom-up neural information it carries. In a contrasting, independent scenario, consider weaker activity in deep layers, corresponding to a lower-amplitude γ packet III and no θ propagation. Now when packet IV occurs in the superficial layers, it encounters a θ peak rather than a trough, diminishing its amplitude and leading to a weaker response. Thus, traveling θ may provide an additional degree of freedom for modulating the strength of the γ packet and its downstream effect. By implementing this principle, the nested spatiotemporal dynamics we uncover here might play an important mechanistic role in organizing broad translaminar modulation of neural activity, beyond the detailed circuit-level organization that stratifies the cortex.

Visual cortical LFPs comprise θ, γ, and 1/f
We first characterized the dominant temporal scales of neural activity, and their organization across laminae and visual cortical regions. As shown in Fig. 2a, we calculated the mean power spectral density of LFP signals across all channels for each cortical region during the flash_250ms stimulus. The resulting spectra exhibit oscillatory peaks above a 1/f background, which we fit using the FOOOF algorithm44 between 1–300Hz. To isolate the oscillatory peaks from the 1/f background, we calculated the ratio of the original spectrum and the 1/f component across cortical layers. This residual spectrum, shown by Fig. 2b for VISl during the flash_250ms stimulus, exhibits peaks at two distinct temporal scales: a slower, low-frequency component centered at ~7Hz, and faster, high-frequency activity between 30–100Hz. Notably, all regions exhibit this broadband γ component above the 1/f background, but only VISp shows an additional narrow-band γ peak18,45 at ~55Hz (both the mains electricity frequency and screen refresh rate were 60Hz). Based on these spectral peaks, we decomposed neural activity into two operational bands for subsequent analysis: a lower-frequency θ band (3–10Hz) and a higher-frequency γ band (30–100Hz). Our primary criterion was to capture most of the spectral power that exceeds the 1/f background across all visual cortical regions and stimuli (see Supplementary Fig. 5 and Supplementary Fig. 6); additional criteria, such as minimizing filter artifacts, are detailed in Methods.

Having identified two dominant frequency bands, denoted as θ and γ, we next mapped the distribution of power in these bands across cortical layers and visual areas. For each band, we calculated the residual power, given by Eq. (5), as the ratio of the band-limited power above the 1/f component and the total power of the signal. During the flash_250ms stimulus, residual θ power increases systematically from superficial to deep layers, as shown in Fig. 2c: across all areas, the median Kendall’s τ between residual θ power and cortical depth ranges from τ = 0.23 (p = 1 × 10−8, 95% CI = [0.15, 0.34], n = 20, two-sided Mann-Whitney U-test) for VISam, to τ = 0.44 (p = 6 × 10−25, 95% CI = [0.38, 0.48], n = 20, two-sided Mann-Whitney U-test) for VISpm. The prominence of θ power in deep layers may reflect a feedback origin, given the role of deep-layer soma in feedback pathways46,47. To assess this feedback hypothesis more directly, we compared the layer-wise residual θ power to the anatomical hierarchy score of each region using Kendall’s τ correlation. These hierarchy scores were calculated for visual cortical areas by Siegle et al.7 from the tracer-based anatomical connectome of Harris et al.9. As shown in Fig. 2e, residual θ exhibits a strong positive gradient along the anatomical hierarchy across all depths, peaking in layer 2/3. We also found similar patterns in θ organization across different stimuli, including the spontaneous (Supplementary Fig. 5) and Natural_Images (Supplementary Fig. 6) conditions.

We then calculated residual γ power and observed no strong trend across cortical depths. Instead, residual γ tends to peak in layers 1 and 2/3 (especially in VISp), as shown in Fig. 2d. The strong γ power in supragranular layers of VISp aligns with recent work emphasizing the role of γ-frequency spiking in the transmission of feedforward, thalamocortical sensory information18,48. Unlike θ, residual γ during the passive flash_250ms stimulus has a weak correlation to the visual hierarchy, as shown in Fig. 2e. However, the hierarchical organization of γ power is highly task-dependent; during the spontaneous stimulus, we found a significant negative hierarchical gradient in most layers (as shown in Supplementary Fig. 5), whereas in the active Natural_Images condition the gradient is positive in all laminae except layer 1 (Supplementary Fig. 6).

We also observed a hierarchical gradient in the 1/f component49,50, with the 1/f α exponent increasing toward deeper cortical layers in all regions (Fig. 2f): the median Kendall’s τ between the 1/f exponent and cortical depth ranges from τ = 0.20 (p = 5 × 10−9, 95% CI = [0.084, 0.25], n = 20 depths, two-sided Mann-Whitney U-test) for VISp, to τ = 0.55 (p = 1 × 10−35, 95% CI = [0.53, 0.60], n = 20 depths, two-sided Mann-Whitney U-test) for VISam. Similarly, the normalized 1/f intercept is enhanced in the superficial layers of VISp and the deep layers of higher-order regions (Fig. 2g). Both the 1/f exponent and intercept are significantly anti-correlated to the anatomical hierarchy in superficial layers, and the intercept is also positively correlated in deep layers (Fig. 2h). Finally, we found that while the 1/f exponent shows behavior-related variation between hit trials and miss trials in the Natural_Images task, the magnitude of this effect is small (see Supplementary Fig. 1).

In addition to the flash_250ms stimulus, we systematically examined spectra for the spontaneous and Natural_Images epochs. A full set of precomputed spectra for all stimuli, mice, and cortical regions (as well as the dorsal lateral geniculate nucleus, LGd) is provided in the accompanying data repositories51,52. This work maps gradients in LFP spectral structure across cortical layers in mouse visual cortex and reveals a tight correspondence between spectral features and the visual hierarchy. The presence of these spectral patterns across laminae, brain regions, and stimulus conditions suggests that θ and γ are general components of mouse visual cortical dynamics that may act as vehicles of information processing across different levels of sensory processing and cognition. However, this stationary, time-averaged view of θ and γ dynamics provided by the power spectrum is insufficient to fully address their functional roles, particularly amid the highly non-stationary dynamics evoked by trial-based stimuli. To resolve these transient dynamics, we next applied wave-based methods to characterize the moment-by-moment spatiotemporal organization of the θ and γ components.

θ waves propagate bidirectionally across laminae and predict task performance
Having identified two dominant timescales of visual cortical LFPs, we next sought to understand the spatiotemporal organization of non-stationary, task-dependent dynamics at each scale. We used wave-based methods to describe the moment-by-moment spatiotemporal structure of the fast and slow LFP components across cortical layers during individual trials. First, in this section, we focus on the slow 3–10Hz band, which we operationally refer to as θ. We used a band-pass filter to isolate the θ component from the raw LFP signal. As visualized for a single trial in Supplementary Movie 1, θ shows clear signs of spatial autocorrelation and propagation along the Neuropixels probes. For every time point, LFP channel, trial, and recording session, we then calculated a θ wavenumber k θ that describes the strength and directionality of translaminar θ propagation. The wavenumber, having units of mm−1, is inversely proportional to the speed of the wave, and takes positive values for propagation from superficial to deep layers and negative values for propagation in the reversed direction. Unlike conventional approaches to analyzing traveling LFP components, such as trial-averaged current source density analysis53,54, our wave-based approach provides instantaneous, single-trial measurements of LFP propagation and is agnostic to layer-wise variation in the magnitude of the LFP signal that may be confounded by differences in the passive properties of the cortical medium, rather than underlying neural activity43. After removing periods of negative frequency (which reverse the relationship between wavenumber and propagation direction, but occur for <2% of trial timepoints), we found that θ-scale LFP appears to travel robustly across cortical layers, and that the visual change-detection task elicits consistent patterns of propagation.

Figure 3a–c illustrate the key patterns of θ propagation by showing the average wavenumber across task trials in VISl. The first key pattern is the non-stationary, bidirectional nature of θ propagation: θ travels from deep to superficial layers following an onset of a change stimulus, but from superficial to deep layers following stimulus offset. Figure 3a shows the wavenumber averaged across hit trials (i.e., Natural_Images presentations in which the mouse correctly identified a change in image content), across cortical depths and time relative to the stimulus. Prior to the onset of the stimulus, between −0.25–0.0s, the average wavenumber is small. By contrast, around the offset of the stimulus, ~0.25s, the average θ wavenumber becomes highly positive, indicating propagation from superficial layers to deep layers. Second, the directionality of θ around the onset of a stimulus is tied to task engagement. During miss trials (for which the mouse failed to identify a change in image content), the average wavenumber is weaker, as shown in Fig. 3b. As a control comparison, Fig. 3c shows the average wavenumber for the flash_250ms stimulus, when mice were not engaged in a task. We expected weaker top-down modulation for this disengaged condition, and we found the average wavenumber in most layers is positive throughout the trial. Furthermore, the strong, bidirectional translaminar propagation during task engagement diminishes in higher-order regions such as VISpm and VISam; we plot the average instantaneous wavenumbers for all visual cortical regions in Supplementary Fig. 7. We also show average current source density maps in Supplementary Fig. 8, highlighting how this conventional analysis uncovers layer-specific patterns but does not detect translaminar propagation. It should be noted that we used a zero-phase bandpass filter to minimize phase distortion. However, as this filter is non-causal it can broaden sharp fluctuations forward and backward in time55. Accordingly, small changes in the wavenumber prior to stimulus onset do not necessarily reflect anticipatory neural activity. Nevertheless, the zero-phase filter preserves peak timings, and we replicated our analyses using a causal filter to confirm that the θ phase relationships between channels are not artifactual (see Supplementary Fig. 15). Together, these analyses suggest that task trials elicit bidirectional traveling θ waves that are closely tied to task engagement.

We next sought to map translaminar θ propagation in different regions of the visual hierarchy, and quantify its relationship to task performance. To do so, we summarized the coherence of θ propagation by an order parameter R θ, which computes, at each moment in time, the mean direction of the θ wave (given by the sign of the wavenumber) across cortical layers. Figure 3d shows the average time course of the order parameter during task trials, for each visual cortical area. Following the onset of the stimulus, the mean order parameter from 0.0–0.15s is most negative in VISl and VISal, reflecting deep-to-superficial propagation in these intermediate areas of the hierarchy. By contrast, after the offset of the stimulus the average order parameter peaks, indicating a more coherent propagating θ wave from superficial to deep layers. This peak is stronger in lower-order visual areas such as VISp; to quantify this hierarchical trend, we computed the individual-level and group-level correlation of the mean order parameter after stimulus offset with the anatomical hierarchy scores of the visual areas7,9. The trial-averaged, within-session post-offset (0.25–0.4s) order parameter is strongly negatively correlated to the anatomical hierarchy, with an individual-level Kendall’s τ = −0.60 ± 0.13 (median ± IQR over sessions, p = 4 × 10−10 against a probe-shuffled null, two-sided Wilcoxon signed-rank test) and a group-level median Kendall’s τ = −0.54 (p < 10−6, 95% CI = [−0.57, −0.50], permutation test, N shuffles = 106).

Does the order parameter—a simple, low-dimensional description of the propagating spatiotemporal structure of low-frequency LFP—indeed contain information relevant to the mouse’s perception and cognition? We trained a regularized linear discriminant analysis (LDA) classifier to differentiate hit trials and miss trials using the order parameter time courses. We found that the classifier, when trained only on the 500ms of data prior to the stimulus offset (and prior to the mouse’s reaction), can predict whether the mouse responded with a median balanced accuracy of 0.64 ± 0.11 (median ± IQR across sessions; 5 folds, 20 repeats), as shown in Fig. 3f. This performance is significantly greater than the null balanced accuracy for probe-shuffled pre-offset data (0.5 ± 0.04), with a median difference of 0.15 (p = 2 × 10−9, 95% CI = [0.11, 0.19], right-tailed Wilcoxon signed-rank test, n = 53). The pre-offset accuracy is also comparable to the performance of a classifier trained on the 500ms of data following the stimulus offset (which can include the mouse’s reaction, termed post-offset; median balanced accuracy of 0.68 ± 0.10). Moreover, the order parameter gives more accurate predictions than the channel-wise mean θ LFP (z-scored over time) in the same pre-offset time window, which has a median balanced accuracy of 0.59 ± 0.09: the median difference between pre-offset order parameter and θ LFP accuracies is 0.044 (p = 5 × 10−5, 95% CI = [0.019, 0.061], right-tailed Wilcoxon signed-rank test, n = 53). That the directionality and coherence of translaminar θ propagation after stimulus onset is more informative of the mouse’s reaction than the average θ LFP suggests that the spatiotemporal dynamics of low-frequency activity are a key signature of the mouse’s perception of visual change. Finally, to determine the relative importance of θ propagation in each visual cortical region, across time, we extracted the weights of an LDA classifier trained on the whole trial period, as shown in Fig. 3g. The pre-offset order parameter, particularly from 50–150ms after the stimulus onset, is most informative of hit/miss performance in three lower-order visual regions: VISp, VISl, and VISal. Here a negative average order parameter and a negative LDA weight indicate that θ propagation from deep layers to superficial layers is more coherent during hit trials. However, after the stimulus offset, VISp has a more positive weight than any other region, indicating that superficial-to-deep propagation in VISp is most informative of the mouse’s behavior after the stimulus is removed.

In summary, offset θ exhibits distinctive signatures of bottom-up processing: (i) it propagates from layer 2/3 and layer 4—keys sites of bottom-up input—to deeper layers (Fig. 3a); (ii) it is strongest and most informative in the lowest-order visual area, VISp (Fig. 3g); (iii) it robustly decreases in coherence along the cortical hierarchy (Fig. 3d); and (iv) it is present even when mice are disengaged from a task (Fig. 3e). Conversely, task-related onset θ appears to be more reflective of top-down processing. First, during task trials, after the stimulus onset, θ propagates from layer 6—which is thought to have a high concentration of cell bodies receiving feedback connections56,57—toward superficial layers. Second, this deep-to-superficial propagation is only present during active task engagement—not during passive flash stimuli—and carries information that can be used to predict the mouse’s response. Hence the deep-to-superficial onset θ may arise from top-down input from higher-order regions, reflecting the mouse’s engagement with the task or expectation of a change in the stimulus. By contrast, the superficial-to-deep offset propagation may reflect a bottom-up effect at the disappearance of the stimulus, temporally synchronizing feedforward information and creating brief quiescent windows to reduce interference with processing in higher-order or extra-visual areas58. Thus we propose that these nonstationary, bidirectional θ waves have dual roles, being associated with bottom-up or top-down processing at different periods during the visual task trials.

Inter-areal θ propagation switches hierarchical directionality
Evoked θ waves travel upward across cortical laminae after stimulus onset, and downward following stimulus offset. These bidirectional propagations also show strong signatures of top-down and bottom-up processing during visual tasks. To clarify the structural basis of this dual functional role of θ, we next studied the propagation of θ waves between regions of the visual cortical hierarchy using anatomical and functional metrics of feedforward and feedback neuronal connectivity7,9.

We first calculated anatomical and functional distance matrices using the anatomical scores (visualized in Fig. 4a) and functional directionality scores (Fig. 4c) published by Siegle et al.7. While the two hierarchies are correlated, they were measured from independent datasets: the anatomical scores reflect the strength of feedforward and feedback axonal connections, whereas the functional scores reflect the average correlation and temporal delay between spiking activity in different regions. To quantify θ propagation between visual areas, we computed the moment-by-moment phase differences between pairs of channels at equivalent cortical depths on different probes. To better estimate the coherence of propagation along the loose hierarchies, we considered all pairs of probes and normalized the phase differences by the corresponding anatomical or functional distance. As for the translaminar wavenumber, we took the negative of the phase differences to ensure that positive values indicate propagation in a feedforward direction, where θ in lower-order regions leads θ in higher-order regions. The mean of the resulting propagation directions gives the anatomical order parameter A θ and the functional order parameter F θ, which describe the coherence and directionality of θ propagation across the respective hierarchies. Crucially, the hierarchical order parameters do not trace the detailed structural projections between regions and layers; instead, they are better viewed as tracking the sequence of neural activity in equivalent layers across the hierarchy, which are thought to have similar functional roles.

We found that task-evoked θ propagation across the hierarchy, like translaminar propagation, switches directionality across a trial: from higher- to lower-order areas (feedback-oriented) following stimulus onset, but from lower- to higher-areas (feed-forward oriented) following offset. Figure 4b shows the mean anatomical order parameter during the task period, which is highly negative (approaching −1) in deeper layers at 0.1s following the onset of the task stimulus, and highly positive (approaching 1) in more superficial layers at 0.1s following the offset. Similarly, Fig. 4d shows the mean functional order parameter, which switches from a minimum value of approximately −0.5 at 0.1s after onset to ~0.5 at 0.1s after offset; in both cases, the peak positive and negative order parameters are significant against a probe-shuffled null. Notably, the switch from feedforward-oriented to feedback-oriented propagation of slow LFP components along both the anatomical and functional hierarchies matches the reversal of translaminar propagation direction between task onset and task offset. We also found that during the flash stimulus, when mice passively view bright or dark flashes, there is no significant inter-areal propagation along either the functional or anatomical hierarchy (see Supplementary Fig. 16).

Together, these complementary findings indicate that θ-scale activity plays a dual role in visual processing, flexibly switching between feedback-oriented (top-down) and feedforward-oriented (bottom-up) modes according to the visual stimulus and task demands. Figure 4e, f summarize how translaminar and inter-areal θ may switch between these two processing modes during task trials. First, as depicted in Fig. 4e, we suggest that stimulus onset triggers a rapid transmission of fine-grained information about the visual scene48 along a feedforward pathway. This feedforward pass appears to dominate γ LFP activity in layers 2/3 and 4 (see Fig. 2), with the earliest signatures of visual input appearing in VISp less than 50ms after stimulus onset7. A short time later, between 50–150ms, a feedback-oriented θ signal arrives in superficial layers of the lower-order visual areas, being transmitted from higher-order regions along a deeper pathway59,60. This feedback θ carries top-down information, likely tied to task engagement (see Supplementary Fig. 3), that travels between layers in each region to meet and modulate feedforward γ in the superficial layers. At stimulus offset, however, θ plays a different role, as shown in Fig. 4f. Offset-evoked θ oscillations may synchronize the feedforward γ packets that travel from lower-order regions to higher-order regions in a superficial pathway. These θ oscillations, which have a characteristic frequency of ~4Hz in VISp (as shown in Supplementary Fig. 2), resemble the evoked activity previously observed in lower visual cortical areas11,58,61,62. Here, we have shown how this offset-evoked θ oscillation propagates across the visual hierarchy and cortical layers, which may serve to periodically shield neural activity in higher regions from the feedforward signal, reducing the interference of bottom-up information with cognitive processing. These findings suggest that the spatiotemporal structure of θ affords it a flexible, context-dependent functional role that adapts to the computational demands of visual tasks.

Localized γ packets are nested within traveling θ waves
Next, we aimed to quantify the spatiotemporal relationship between θ waves and γ activity to understand the mechanisms by which broad θ waves may modulate and multiplex local neural processing. We hypothesized that fast activity—dominated by γ oscillations—would be structured as short-lived, localized packets of heightened oscillatory power. To confirm the presence of spatiotemporal γ packets, we used univariate peak-finding methods illustrated in Fig. 5a to identify local maxima in the γ amplitude r γ of each LFP channel. We found that γ packets are short in duration—typically (40.0 ± 5.5) ms (mean ± standard deviation across all trials and sessions) during the Natural_Images stimulus. These packets are also localized across space, having a typical spatial width of 457 ± 102 μm, although their size varies across the trial periods. As shown in Fig. 5b, which plots the trial- and session-average spatial width of γ packets across the trial period and visual cortical areas, bursts are narrower [357 ± 62 μm] in the 250ms following stimulus onset and wider [476 ± 102 μm] in the 250ms following offset. Moreover, the width of γ packets following stimulus offset decreases across the anatomical hierarchy at both the individual level, with Kendall’s τ = −0.47 ± 0.40 (median ± IQR across sessions, p = 6 × 10−6 against a probe-shuffled null, two-sided Wilcoxon signed-rank test), and the group level, with a median Kendall’s τ = −0.28 (95% CI [−0.33, −0.21], p < 10−6, permutation test, N shuffles = 106). The packet width is a biologically sensible index of γ activity, with γ bursts being tied to local stimulus luminance on the retinotopic map and circuit-level motifs across laminae48. However, to our knowledge no earlier work has analyzed the spatial size of γ bursts in response to visual stimuli. We propose that increased localization of γ packets following stimulus onset may reflect multiplexed processing of many local features of the visual stimulus.

Having found that γ oscillations are organized into localized packets, we then asked how the timing of γ packets is related to the low-frequency translaminar waves. We studied three aspects of the relationship between θ traveling waves and γ packets: (i) how cross-scale coupling varies with frequency, using standard phase–amplitude coupling methods; (ii) trends in PAC across layers and cortical areas, averaging over time; and (iii) the nonstationary properties of nested oscillations in response to a visual stimulus, by computing PAC jointly across space and time during the course of an average trial.

First, we analyzed global phase–amplitude coupling (PAC) to determine whether our chosen θ and γ bands exhibit cross-scale interactions. We computed the classical phase–amplitude modulation index63 across a range of phase and amplitude frequency bands, averaging over all channels and time points for each visual area and stimulus (see the full set of comodulograms in Supplementary Fig. 9). We found that all areas tend to exhibit strong PAC at one or both of two phase–amplitude frequency pairs. First, phase frequencies near 7Hz and amplitude frequencies from 50 to 100Hz tend to be coupled during the spontaneous conditions, and in higher visual areas during task and flash trials. Second, phase frequencies near 4Hz and broadband amplitude frequencies >30Hz are coupled in lower-order areas during flash and task trials. These frequency bands closely agree with the power spectral peaks in Fig. 2, suggesting a coupling between θ waves and γ packets and supporting our chosen bands of 3–10Hz (θ) and 30–100Hz (γ).

We next examined how phase–amplitude coupling—the strength of the relationship between θ waves and γ packets—changes across cortical layers and regions, establishing cross-frequency organization aligned with the visual hierarchy. Using the θ and γ bands we identified in Fig. 2, we computed the modulation index (across time) for each channel in each visual cortical region during the spontaneous stimulus. Shown in Fig. 5c, the layer-wise PAC is strongest in the supragranular layers of lower regions and in the infragranular layers of higher regions. We also found that γ packets are phase locked to the trough of the θ cycle in lower-order regions (VISp, VISl), but to the peak of θ in higher-order regions (VISrl, VISam), as summarized in Fig. 5d. That γ packets are differentially locked to θ phases between layers and cortical areas suggests that θ may represent a general principle for coordinating γ activity across the visual hierarchy. However, classical PAC metrics cannot capture detailed spatiotemporal interactions between θ and γ, and thus provide poor descriptions of the non-stationary activity that occurs during task trials.

Following our study of the spectral and spatial features of PAC, which shows how propagating θ waves carry localized γ bursts, we sought to determine if PAC fluctuates over the course of a trial. To resolve phase–amplitude coupling over space and time, we computed the modulation index for each channel at each time point, taking the distribution of phases across trials and subjects. The resulting spatiotemporal modulation index is high when γ packets in a given spatial location, at a particular point in time relative to the stimulus onset, occur at a consistent phase of the θ cycle. We show the spatiotemporal modulation index during the average task trial for a lower-order area (VISl) in Fig. 5e, for a higher-order area (VISam) in Fig. 5f, and for all other visual cortical areas in Supplementary Fig. 10. In both lower-order areas, VISp and VISl, the spatiotemporal modulation index peaks in layers 2/3 and 4 following stimulus offset (Fig. 5e), aligning with the emergence of broad γ packets (see Fig. 5b) and also θ waves propagating from superficial to deep layers. This heightened coupling between θ and γ in layers 2/3 and 4 after stimulus offset also agrees with the feedforward orientation of offset θ, given the strong feedforward cortico-cortical network that resides in these layers64. Also in these lower-order areas, the spatiotemporal modulation index peaks soon after stimulus onset, in layer 1. Layer 1, which houses apical tufts of pyramidal cells that receive feedback input65 and cell bodies of GABAergic inhibitory interneurons66, plays a critical role in the simultaneous integration of top-down and bottom-up information65–67. In higher-order areas, however, spatiotemporal θ–γ coupling is strongest in deep layers following stimulus onset (Fig. 5f), especially layer 6, which houses a strong cortico-cortical feedback network64. That layer 1 (in lower areas) and layer 6 (in higher areas) have a high θ–γ PAC following stimulus onset further indicates that onset θ reflects top-down propagation meeting bottom-up γ packets in the superficial layers. Overall, these results demonstrate how γ packets and θ waves are nested across spatial and temporal scales, with the flexibility to switch directionality and adapt dynamically to changing demands of visual processing.

Nested θ and γ reflect underlying patterns of neuronal spiking
Finally, we asked how the θ and γ oscillation are coupled to underlying spiking activity. For θ, we computed a spike–phase coupling statistic—the pairwise phase consistency (PPC)68—between spikes from each unit (putative individual neurons) and the θ LFP phase from the closest LFP channel. Figure 6a shows the the average PPC during the Natural_Images stimulus, showing strong spike–phase coupling in the supragranular layers of lower-order areas. A similar trend is reflected in spike–amplitude coupling between γ amplitudes and neuronal spiking, quantified as the mean normalized r γ across spike times. As plotted in Fig. 6b, the distribution of spike–amplitude coupling is strongest in the supragranular layers, as for spike–phase coupling. Together, these high values of spike–phase and spike–amplitude coupling corroborate our earlier findings of nested LFP oscillations, showing that the simultaneous cross-scale coordination of θ and γ extends to the level of precise neuronal spiking.

To elucidate the mechanistic relationship between spatiotemporal θ components and neuronal spiking, we calculated the average preferred θ phases across visual cortical areas and layers, plotted in Fig. 6c. We found that neurons in supragranular layers prefer troughs of the θ cycle, whereas neurons in infragranular layers prefer peaks; this trend is repeated across areas in the visual hierarchy, and broadly matches the preferred phases of γ packets (cf. Fig. 5). Do these preferred phases reflect a functional role of traveling θ components? We selected neurons that were sensitive to θ, with a pairwise-phase consistency of greater than 0.1 across the entire Natural_Images stimulus epoch, then asked how layer-wise differences in preferred θ phases and firing rate change between hit and miss trials. We found that units in supragranular layers of the four lowest-order areas have a significant tendency to prefer greater θ phases (closer to θ troughs) during hit trials, as shown in Fig. 6d. The greatest change in preferred phase occurs in the superficial and middle layers of intermediate areas, with the largest effect in layer 4 of VISl (rank biserial correlation of 0.85). By contrast, deeper layers show no significant change in preferred phases.

If neuronal spiking in supragranular layers becomes entrained to the trough of θ during trials where θ has the strongest deep-to-superficial propagation (see Fig. 3), is there a corresponding change in firing rate that would suggest the propagating θ trough is linked to enhanced neuronal activity? To address this question, we calculated the mean firing rate of θ-sensitive neurons, grouped into infragranular, granular and supragranular compartments, shown in Fig. 6e for hit trials and Fig. 6f for miss trials. We found that supragranular units have a significantly enhanced peak firing rate during hit trials for the four lowest-order areas (VISp–VISal), occurring at ~90ms after stimulus onset, with an effect size (rank biserial correlation) of between 0.37 (VISp; moderate effect) and 0.68 (VISrl; large effect). Collectively, these findings support the spatiotemporal θ–γ scheme summarized in Fig. 1a: during successful change detection, a θ trough carrying top-down information travels upward from deep layers and coincides with enhanced bottom-up γ packets and neuronal spiking in superficial layers. These task-related spike–phase relationships come together to form the picture of nested dynamics presented in Fig. 1a: that a traveling θ trough moves upward from deep layers during successful change-detection trials, where it coincides with enhanced neuronal responses in shallow layers.

Methods
LFP preprocessing
We obtained LFP time series from the Allen Neuropixels—Visual Behavior dataset for three stimuli: (i) spontaneous, a uniform gray screen with a duration of ~5 min; (ii) flash_250ms, 250ms bright or dark flashes spaced by 2s for ~5 min; and (iii) Natural_Images, corresponding to an active task in which natural images are flashed for 250ms, spaced by 500ms, for ~60 min (we did not analyze passive task replay). Individual trials were defined as time periods consisting of the 250ms prior to the stimulus onset, the 250ms during which the stimulus was presented, and the 500ms following stimulus offset. After filtering for channels located inside the cortex, we rectified their irregularly sampled LFP timestamps by rounding to a relative tolerance of 10−6, before mapping stimuli and spikes to the same timestamps. The resulting regular time series had a sampling period of Δt = (1/1250)s.

Power spectra
We calculated, for every LFP channel, the power spectral density P(f), where f is the frequency in Hz and P has units of V2/s, using a Welch periodogram set for a frequency spacing of 0.1Hz and a frequency range of [1 300] Hz (with a Hann window, 10,000 samples of zero-padding, and 50% window overlap). We identified spectral peaks with a local peak-finding algorithm, filtering by a relative prominence of 0.25; the relative prominence, measured as the height the local maximum above the local baseline divided by the absolute height of the peak, accounts for differences in base-level power across frequencies. We then fit the power spectral density for each channel, using the FOOOF algorithm44, to a model given by:

Pb(f) = exp(L(f) + Σn=1N Gn(f)), (1)

where Pb(f) has units of V2/s and n is the number of Gaussian peaks (at most N=8), each of the form:

Gn(f) = an exp(-(f - cn)2 / 2wn2), (2)

where an is the amplitude, cn is the center frequency, and wn is the bandwidth of the nth Gaussian component. L(f) is the 1/f component, given by:

L(f) = b - log(k + fα), (3)

where α is the exponent of the 1/f tail, b is an offset parameter, and k is a parameter that introduces a knee at low frequencies. Equation (3) was used to calculate the residual spectrum:

R(f) = P(f) / L(f), (4)

which is reported in decibels (see Fig. 2). Finally, we calculated the total residual θ and γ power as the proportion of power above the 1/f component in the θ and γ band:

RTot = ∫f1f2 [P(f) - L(f)] df / ∫f1f2 P(f) df, (5)

where (f1, f2) = (3Hz, 10Hz) for θ and (f1, f2) = (30Hz, 100Hz) for γ.

Statistics and reproducibility
When plotting power spectra, we calculated averages using the mean and quantified uncertainty using the standard deviation across sessions. For most other quantities—as specified in the main text—average values were calculated using the median and uncertainties were quantified by the 95% confidence interval across sessions. Confidence intervals were estimated by bias-corrected and accelerated bootstrapping with balanced resampling over 10,000 draws107. We applied the same bootstrapping method to estimate confidence intervals for circular averages.

To calculate the correlation between the anatomical hierarchy scores (retrieved from Siegle et al.7) and various features of neural dynamics, such as LFP spectral content or θ propagation coherence, we used Kendall’s τ, a tie-robust rank-based correlation coefficient, at an individual level and a group level, as in Harris et al.108. At the individual level, we computed Kendall’s τ between the anatomical hierarchy scores and neural features for each session separately, then calculated the median τ. To assess significance, we generated surrogate samples by shuffling the hierarchy ranks of each probe, ensuring that channels on the same probe had the same hierarchy rank even in the shuffled data, then computing τ for the shuffled ranks. To estimate p values, we performed a Wilcoxon signed-rank test between pairs of observed τ correlations and the surrogate samples, then adjusted for multiple testing using the Benjamini-Hochberg method. At the group level, we instead computed Kendall’s τ between the hierarchy scores and neural features pooled across sessions. We then computed p values with a permutation test, shuffling data in the same way as the individual-level analysis then calculating the proportion of shuffles that had a correlation greater than the true measured value.

Sample size and excluded data were determined by the session filtering criteria detailed in Session filtering; no statistical method was used to predetermine sample size or select data for exclusion. The experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment.

θ traveling waves
Based on our spectral analysis, we chose to operationally define the θ band as the frequency range between 3Hz and 10Hz. This choice of θ frequency was informed by three main factors. First, we aimed to cover the majority of low-frequency peaks above the 1/f component across all stimuli and visual areas (as shown in Fig. 2, Supplementary Fig. 5, and Supplementary Fig. 6) to fully resolve the task and region dependencies of the slow LFP component. Additionally, we show in Supplementary Fig. 2 how the dominant frequency of θ oscillations evoked by stimulus offset can vary between 4Hz and 6.5Hz across visual regions, suggesting a broader range of frequencies is needed to simultaneously capture the dominant low-frequency dynamics across visual areas. Second, since our methods are based on inferring propagation from the Hilbert transform of the spatiotemporal signal, the frequency range must be narrow enough to avoid significant negative-frequency components. Third, as we aimed to analyze task-related dynamics evoked by visual stimuli with a short (~250ms) duration, we chose a frequency band wide enough to reduce non-causal filter artifacts. This choice of θ also covers the majority of low-frequency bands that have been linked to visual processing in mice, including 3–7Hz61, 0–10Hz109, 3–5Hz58, 3–6Hz35, 2–8Hz34, 3–8Hz110, 4–8Hz111, and 3–10Hz112. We repeated our analyses for Fig. 3 and Fig. 4 using two narrower frequency bands: 3–5Hz (Supplementary Fig. 13) and 6–10Hz (Supplementary Fig. 14). We found qualitatively similar results but weaker behavioral correlations and decoding accuracy, with the 6–10Hz band showing no significant inter-areal propagation. We also found similar patterns of θ propagation for data filtered with a causal 8th-order Butterworth filter, with timings shifted by ~125ms due to the expected phase offset (see Supplementary Fig. 15).

Next, to measure θ propagation across cortical layers, we used wave-based methods that have previously been applied to horizontal cortical propagations measured by various recording modalities4. We represent a real spatiotemporal signal in a conventional complex form given by:

a = r exp[iϕ] = r exp[i(kx - ωt + b)], (6)

where r(x,t) is the instantaneous amplitude, ϕ(x,t) is the phase, k(x,t) is the wavenumber (or angular spatial frequency), ω(x,t) is the angular temporal frequency, and b(x) is a phase offset. Since ϕ = −kx + ωt − b, the angular temporal frequency ω(x,t) is the temporal gradient of the phase:

ω = ∂ϕ / ∂t (7)

Similarly, the instantaneous angular wavenumber is the negative spatial phase gradient:

k = - ∂ϕ / ∂x (8)

Together, the wavenumber and temporal frequency determine the direction of wave propagation. By Eq. (7), phase increases over time for signals with a positive instantaneous frequency, but by Eq. (8), phase decreases in the positive x direction for a positive wavenumber. Therefore, when both the instantaneous frequency and wavenumber are positive, a given wavefront (such as the contour of ϕ = 180∘, a trough) will begin at a smaller x position and move to larger x positions as time passes, and the wave as whole will travel in the positive x direction. When applied to the θ signal, the instantaneous wavenumber k quantifies the local, instantaneous directionality and coherence of θ propagation; positive k indicates the wave is traveling from superficial to deep layers, whereas negative k indicates deep-to-superficial propagation.

The transformation of a real signal to a complex wave [Eq. (6)] is not uniquely defined29,113,114: we used the Hilbert transform, H, which is a linear operator defined as a convolution with the Hilbert kernel, 1/πt. First, we extracted the low-frequency θ time series from each raw LFP time series using a fourth-order, zero-phase, non-causal Butterworth bandpass filter between 3Hz and 10Hz. For convenience, we notate this filtered time series as θ(x,t), where t is the time, in seconds, and x is the spatial position along the probe, in μm. The analytic signal is defined from the Hilbert transform as:

a(t;x) = y(t;x) + iH[y(t;x)] = r(t;x) exp[iϕ(t;x)], (9)

where y(t;x) denotes the univariate time series from a given spatial position x; the Hilbert transform is applied independently to each LFP channel. We applied Eq. (9) to the θ time series from each channel, giving the phasor a θ(x, t), and the instantaneous θ phase ϕ θ(x,t) = Arg[a θ(x,t)]. Finally, we used Eq. (8) to calculate the θ wavenumber k θ(x, t), taking central finite derivatives with circular differences.

We note that inferring the propagation direction from the wavenumber relies on the assumption that the instantaneous frequency is positive; the propagation direction, as measured by the phase velocity (or, the ratio of the instantaneous frequency and the wavenumber) has an opposite sign to the wavenumber when the instantaneous frequency is negative. For the analytic signal, based on the Hilbert transform, negative frequencies indicate periods where the signal contains more than one simultaneous spectral component. For the filtered θ signal, we found that periods of negative frequency were rare (occurring for <2% of trial timepoints) and that directionalities inferred from either the wavenumber or the phase velocity agreed closely. Since the wavenumber is more numerically stable than the phase velocity, and less sensitive to noise, we present results based on the wavenumber in the main text, and mask any negative frequency periods from our wave-based analyses.

Order parameter
To summarize the overall directionality of θ propagation along cortical layers, we introduce an order parameter, R θ(t), calculated as the average direction of the wavenumber:

R θ(t) = 〈Sgn[k θ(x,t)]〉 x, (10)

where 〈⋅〉 x is the mean over channels in a given probe (excluding wavenumbers that co-occur with negative instantaneous frequencies). The order parameter R θ is bounded between [−1, 1]; R θ=1 indicates that θ propagation is contiguously directed from superficial to deep layers, and R θ=−1 indicates contiguous propagation in the opposite direction. |R θ| < 1 indicated that the directionality of θ is mixed across cortical layers.

To quantify the relationship between the θ propagation and the mouse’s response to the visual change detection task, we used a linear discriminant analysis (LDA) classifier to predict hit trials (correct response to a change in the stimulus) and miss trials (no response to a change in the stimulus) based on the order parameter across time and visual cortical regions.

For each session, we flattened the order parameter sampled across time, channels, and visual cortical areas for each trial. Collecting the flattened data yielded an nf × nt, feature × trial input matrix X, where each feature corresponds to the order parameter at a given time point, for a given channel, in a given structure. To reduce the dimensionality of the input data, we coarse-grained over time by taking the mean R θ across groups of 8 time points. As a final preprocessing step, we standardized the input data to have zero mean and unit variance. For these normalized data, LDA finds the hyperplane that optimally discriminates the two classes (hit and miss) under assumptions of equal variance and within-class normality. LDA uses the scatter matrix of each class, defined for hit trials as:

Sh = (Xh - μh)(Xh - μh)T, (11)

where Xh is the nf × nh feature × trial matrix representing hit trials. Miss trials are defined similarly for miss trials. These single-class scatter matrices capture information about the covariance of each class, summarized by the within-class scatter matrix:

Sw = Sh + Sm + ϵλ1 I, (12)

where ϵ is a regularization parameter, λ1 is the largest eigenvalue of Sh + Sm, and I is the identity matrix. LDA seeks to find a projection that maximizes the separation between classes relative to the between-class scatter matrix defined from the class centers as:

Sb = nh(μh - μ)(μh - μ)T + nm(μm - μ)(μm - μ)T, (13)

where μ is the n column vector of feature means over all trials. The optimal projection is found by solving the generalized eigenvalue problem given by:

Sb w = λ Sw w, (14)

where w is an eigenvector with eigenvalue λ. Predictions are made by projecting data onto the leading eigenvector w1, then assigning positive values to the hit class and negative values to the miss class. Thus the elements of w1 are weights that indicate the importance of R θ at each channel and time point t to the overall classification, with negative weights indicating R θ is more negative for the positive hit trials. For each condition, we evaluated the accuracy of the classifier using 5-fold cross validation. We estimated a suitable regularization parameter ϵ via particle swarm optimization of the balanced accuracy for pre-offset data, giving ϵ = 0.5. We used the same regularization parameter and the same dimensionality nf (determined by the duration of the data) for all classification conditions to ensure an equivalent comparison. We also removed any hit trials with a reaction time of less than 0.25s, ensuring that all data points prior to stimulus offset also precede the mouse’s behavioral response (although the proportion of trials with reaction times less than 0.25s was low, typically less than 2%).

Hierarchical θ propagation
To determine how translaminar θ relates to low-frequency propagation across the cortical hierarchy, we compared phase differences between channels in different visual cortical regions against their hierarchical positions. Recognizing that many functional and anatomical properties of cortical layers are conserved across the cortex, we chose to compare phases at equivalent cortical depths across different visual areas. For each session, we selected the set of 19 channels of each probe closest to each element of a unified set of cortical depths: 5%, 10%, …, 95%. We denote the analytic phase (at a given time point) of the set of six channels closest to a given depth, one from each region, as ϕ s(i), where i corresponds to one of the six probes (or regions). To compared differences in θ phase to relative hierarchical levels, we retrieved the anatomical hierarchy scores and functional hierarchical distances published by Siegle et al.7 (where the anatomical scores had been recomputed from published by Harris et al.9). We computed the pairwise differences of anatomical hierarchy scores to generate an anatomical hierarchical distance matrix matching the functional distance matrix. Using the anatomical distance matrix, annotated as a(si) for a given region i and rescaled so that the maximum hierarchical distance was 1, we calculated an anatomical hierarchical order parameter as:

A θ = 〈Sgn[a(ϕj) - a(ϕi)]〉 i,j, (15)

where i, j ∈ {1, 2, …, 6}, i ≠ j (corresponding to six visual cortical regions), we use the circular statistics to evaluate the difference in phases, and the negative sign accounts for phase increasing over time, as in Eq. (8). The magnitude of A θ quantifies the consistency of θ propagation across the visual hierarchy, with larger values indicating more consistent propagation; a small A θ occurs when θ phases are oriented inconsistently across the hierarchy. The sign of A θ indicates the direction of propagation: positive values indicate that θ propagates from higher to lower hierarchy levels, whereas negative values indicate propagation in the opposite direction. Similarly, we applied Eq. (15) to the functional distance matrix to compute the functional hierarchical order parameter F θ. Note that we masked and discarded periods of negative frequency, as in as for analysis of the θ wavenumber.

γ packets
To identify periods of high γ power, we filtered each LFP time series in the band 30Hz to 100Hz using a high-frequency bandpass filter, giving the time series γ(x, t). Then, as for the low-frequency θ oscillation, we calculated the analytic signal of γ, a γ(x, t), using Eq. (9). However, rather than extract the instantaneous phase of γ, we calculated the analytic amplitude, r γ(x, t) = |a γ(x, t)|. To detect γ bursts, we first z-score normalized the analytic amplitude, giving rγs, then took the mean as:

Nh(rγ) = rγ - min(rγ) / σ(rγ) * (1 - π/2), (16)

where σ(rγ) is the standard deviation of rγ and min(rγ) is the minimum value. We identified temporal bursts as local maxima of Nh(rγ) using a univariate peak-finding algorithm over time, with a window width of 1 and a minimum prominence of 2. For the resulting channel-wise temporal bursts, channel-wise durations were calculated as the full width at half maximum, taking the height of the smallest adjacent minimum as the zero level. Overall channels, we masked time points that fell within a window as wide as the duration of each given channel-wise temporal burst, generating a mask of spatiotemporal γ packets. Finally, we identified connected components of the packet mask to infer total packet durations and spatial widths from their bounding boxes.

Phase–amplitude coupling
To quantify the relationship between slow-frequency θ oscillations and the high-frequency γ bursts, we used the phase–amplitude coupling (PAC) modulation index proposed by Tort et al.63. θ-γ phase–amplitude coupling reflects the degree of coincidence between the timing of high-frequency γ bursts and select phases of the low-frequency oscillation, or the entropy of the distribution of θ phases conditioned on the amplitude of γ bursts. To calculate the modulation index given a set of θ phases ϕ θ and corresponding γ amplitudes r γ, we partitioned γ amplitudes into nb equally spaced bins according to their corresponding θ phase values:

C = {ck : ck = (2k-1)/nb * π, k=1, 2, ..., nb} (17)

We then calculated the mean γ amplitude for a given phase bin centered at ck as:

<r(c)> = 〈rγ(i)〉 i, rγ(i) ∈ [c-π/nb, c+π/nb] (18)

After normalizing the mean γ amplitudes to have unit sum, giving br = <r> / Σ<r>, the PAC modulation index is the normalized Kullback-Liebler divergence between br and a uniform distribution:

PAC = ln(1/nb) * Σ c∈C br c ln(nb * br c), (19)

where PAC is bounded by [0, 1]; large values indicate a strong relationship between θ phase and γ amplitude. From the distribution br, we also calculated the preferred phases of θ as argmaxc br c, the phase bin center corresponding to the largest mean γ amplitude.

We calculated the PAC modulation index in three ways: over θ and γ central frequencies, across time for each LFP channel, and across trials for every channel and time point. To compute PAC over frequency bands, and construct phase–amplitude comodulograms as plotted in Supplementary Fig. 9, we filtered the LFP signals y in phase-frequency bands centered at fp = 3, 3.25, … 12 with widths of 2Hz, and amplitude frequency bands at fa = 25, 1, … 125 with widths of 20Hz. We found the instantaneous phase for each phase-frequency band, and the then the instantaneous amplitude for each amplitude-frequency band using Eq. (9), before calculating the PAC between each (fp, fa) pair. We then calculated the PAC over time—for each channel and session independently—during the spontaneous stimulus, using our original θ phases ϕ θ (corresponding to 3–10Hz) and γ amplitudes r γ (corresponding to 30–100Hz). We present this time-grouped PAC in Fig. 5c. Finally, to understand non-stationarity in θ–γ coupling over the course of a task trial, we calculated the trial-grouped PAC between ϕ θ and r γ for each channel, trial-relative time point, and session during the hit and miss trials of the Natural_Images stimulus (as plotted in Fig. 5e).

Spike-LFP coupling
As the final component of our proposed nested dynamical motif, we quantified the relationship between neuronal spiking and the θ and γ LFP. The Neuropixels data record spike times for hundreds of neurons, identified via spike sorting that assigns each detected spike to a given unit—each unit is likely to be an individual neuron. To calculate θ spike-phase coupling, we first extracted the θ phase from the LFP channel closest to each unit, at the closest time point to each spike, giving the set of spike-timed θ phases ϕ s. We then computed the spike-phase coupling using the PPC68:

SPC = 1 / (ns(ns-1)) * Σi=1ns Σj=i+1ns cos(ϕ s(i) - ϕ s(j)), (20)

where ns is the number of recorded phases. Negative SPC values occur when the mean size of differences in phase is greater than π/2. Since the lower bound of the average difference in phases decreases as sample size increases, negative PPC values are indicative of low sample sizes (i.e., sparse-firing neurons). We replaced these negative values with NaN to avoid bias by low sample sizes. To calculate the preferred θ phase for each unit, we used the angle of the resultant vector after averaging unit vectors with phases of ϕ s(i), as:

ϕ s = Arg [Σi=1ns exp(iϕ s(i))] (21)

For γ spike-amplitude coupling (SAC), we extracted γ amplitudes at each spike time from the closest LFP channel to each unit and applied a half z-score normalization [Eq. (16)], giving brs, then took the mean as:

SAC = 〈brs(i)〉 i (22)

Finally, we calculated the time courses of firing rates during hit and miss trials for individual units that were sensitive to θ (pairwise-phase consistency of >0.1), as presented in Fig. 6e, f. We calculated the peristimulus time histogram (PSTH) for each unit by aligning spike times to the onset of the stimulus change in each trial and binning the spikes in windows of 15ms. To normalize the firing rates for each unit, we z-scored it relative to the unit’s mean firing rate and standard deviation, which were calculated across the entire trial period for all trials. We then averaged z-scored PSTHs across hit or miss trials for each unit, before averaging across units into the superficial (L1, L2/3), middle (L4), or deep (L5, L6) compartments.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, Harris and Gong analyzed high-density electrophysiological recordings from the mouse visual cortex. They proposed that, during stimulation, a cross-scale oscillatory dynamic emerges, involving nested theta-gamma waves that are hierarchically organized and run alongside spiking activity. After disentangling the Local Field Potential (LFP) signal into theta, gamma, and the 1/f components, they first show that there are two distinct theta oscillation components: one that is offset-related, propagating from layer 2/3 to deeper layers, and another that is onset-related, reflecting feedback signaling. The offset-related theta has a high amplitude when the animal is disengaged from the task, whereas the onset-related theta appears only during active tasks. This task-evoked theta exhibits hierarchical organization and produces spatiotemporal dynamics compatible with previously described traveling waves. Gamma oscillations are organized in packets nested within theta, structured across both cortical areas and layers. This nested activity correlates with spiking activity during stimulation, suggesting a flexible, spatiotemporal mechanism for information multiplexing and modulation.\n\nThe manuscript analyzes a high-quality dataset and presents intriguing results and ideas. It provides a detailed description of theta dynamics across the mouse visual cortex, which aligns well with recent views on theta and gamma oscillation dynamics. However, it lacks a mechanistic characterization of how these nested oscillations contribute to information processing beyond their alignment with spiking activity. This limitation weakens the manuscript's central message. The following suggestions aim to strengthen the manuscript’s impact:\n\nThe primary concern is that, while the authors demonstrate significant LFP-spike coupling at both frequencies of interest, it remains unclear whether this coupling influences perceptual information processing or transmission. For example, the authors correctly perform a PPC analysis on neighboring LFP channels relative to the spike cluster of interest. However, the functional relevance of this coupling remains uncertain. Are spike-field locked clusters functionally meaningful? Can perceptual information be accurately decoded from spike-LFP coupling? Does this coupling facilitate information transfer across different cortical areas? To enhance the manuscript’s claims, the authors could quantify a potential causal relationship between LFP oscillations and spiking activity. This analysis could complement their results in Figure 3F on hit/miss decoding of theta waves and offer a more detailed characterization of the putative causality underlying their findings.\n\nA second concern is that, while the manuscript focuses on nested theta/gamma oscillations, certain additional analyses feel somewhat disconnected from the main narrative. For example, the decomposition of the LFP signal into its oscillatory and 1/f components is highly interesting, yet it is unclear whether the 1/f component plays any functional role beyond its description. The authors could strengthen their hypothesis by contrasting the functional relevance of theta/gamma nested activity against 1/f activity. If functional differences between these components can be identified, this would further support the manuscript’s key claims.\n\nA third concern relates to the origin and dynamics of theta activity. The supplementary video suggests that the traveling offset-related theta wave may have a subcortical (hippocampal?) origin. If so, some of the observed theta activity may not be intrinsic to the visual cortex but instead volume-conducted from subcortical regions. The authors should clarify whether this theta activity truly propagates within the cortex or is merely an effect of volume conduction. Although completely ruling out volume conduction is difficult, the authors could quantify phase differences across cortical areas to check for zero-phase relationships, which would support the hypothesis of volume conduction.\n\nFinally, while the discussion section is well-informed, it could reach a broader audience by providing a more comprehensive conceptual framework for the manuscript’s findings. For example, the authors correctly cite Lisman’s work on nested theta/gamma oscillations, but they should elaborate on the current state of this theoretical framework. What recent findings support or challenge this perspective? Additionally, the discussion should address well-known limitations of cross-frequency coupling analyses and clarify how this study mitigates these limitations.", "reviewer_2": "The paper analyzes electrophysiological data from mouse visual cortex during a visual task, with or without engaging attention. The authors show theta-band traveling waves propagating at different scales between regions and cortical layers. They also relate the wave's phase to the gamma-band activity and spikes.\n\nOverall, the paper is methodologically sound, and the results are interesting. I'm not entirely sure about the novelty regarding the previous literature about gamma-theta coupling, but the specificity may apply to the specific animal model used.\n\nBelow are some points that may need further clarification for a broad audience.\n\nHow does the intra laminar theta relate to between areas theta? More precisely, I don't understand the assumption that is used to calculate the order parameter, i.e., that the average translaminar theta propagation in one region informs about the theta propagation between areas.\n\nIn figure 1A, the relationship between packets I and II, and between packets III and IV is unclear. Clearly, they are all modulated by the slow theta wave, but it's unclear how they are related to each other. Doed packet I (III) cause packet II (IV) ? Is the fact that II and IV are enhanced and inhibited unrelated to the presence of packets I and III?\n\nThe residual power in Figure 2F seems similar across cortical regions (except for VISam). Isn't this the opposite of what one would expect considering the proposed feedback/top-down role?\n\nIt is unclear how the anatomical hierarchy used to obtain the correlations is defined, and how the results can be interpreted (figure 2E and 2H) other than showing a (well-known) dynamical cortical gradient.\n\nFrom an interpretational point of view, what does it mean that theta propagates from deeper to superficial layers? Shouldn't one expect the opposite, considering that the (theta activity in the) deeper layer encodes expectations rather than sensory information? Would this reflect mostly top-down modulation rather than sensory processing? Could the author decode some stimulus information from the theta activity and its propagation?\n\nI don't think there is evidence to support the authors' proposal (line 262) that theta provides feedforward processing: shouldn't one observe a forward pass also at the onset of the stimulus, then?\n\nFrom the power spectra, they analyze broadband gamma (30-100 in equation 5). Would the authors imply this is a fairly large broadband effect (possibly present in high-beta as well?), or would it be possible to localize a narrow band gamma frequency?\n\nFrom the method, it is unclear why the instantaneous angular wave number (i.e. the spatial frequency) is the negative spatial phase gradient (a mathematical or even an intuitive explanation would very much help). It's also unclear why this value would quantify the directionality of the wave.\n\nIt is also unclear why the 'periods of negative frequency -defined as the derivative of the phase- were rare'. What does this imply? That the directionality was mostly consistent across trials/electrodes? Can the author quantify and report how many trials were masked (as specified in lines 568-569)?\n\nOverall, it's unclear the relationship between translaminar theta and between regions. The results are extremely similar, as pointed out by the authors (lines 291 onward). Since the laminar direction is roughly orthogonal to the cortical sheet, one may wonder how and why the two results relate so strikingly.\n\nHow do the authors interpret the findings that spikes are aligned to the trough and peaks of theta at different cortical layers? This speaks against the oscillatory voltage modulation (i.e. neurons are more likely to fire at the peak of an oscillation), and it's also difficult (for me) to interpret if we consider theta as a translaminar traveling wave: depending on the wave's speed, spikes won't be synchronized by the traveling wave.", "reviewer_3": "Harris and Gong have performed an interesting and novel characterisation of oscillatory activity in primary and higher visual areas of the mouse cortex by analysing data collected by the Allen Institute. The authors provide novel insights into integration of visual information along cortical depths and across visual areas. It is proposed that theta and gamma oscillation carry feedback and feedforward related information at different temporal epochs. Specifically, they identify short gamma bursts that occur based on feedforward processing, which in turn multiplex with travelling theta waves in a feedback direction after the onset of the stimulus in a change detection task.\n\nOverall, the analyses are very interesting and bring new perspectives to the field. However, there are many concerns with the interpretation of the results in terms of feedback and feedforward signals both across depth and across areas. One additional aspect that is really interesting, but doesn’t seem to be the focus of the study is the difference between task related activity and non-task related flashes. While we are enthusiastic regarding the potential of this work, there remain many concerns to be resolved.\n\nMajor Concerns\n\n1. Theta is defined to be between 3-10Hz. This is a relatively wide range of low frequencies which happen to include both 3-5Hz rhythms and theta (6-9Hz). These 3-5Hz rhythms, which have been likened to human alpha (see Senzai et al. 2019 and Nestvogel and McCormick 2022), are a dominant feature of mouse visual cortex and should be considered when performing these analyses. We suggest the authors reanalyse the data as they have done throughout figures 1-3, but only considering a narrower theta rhythm (6-10Hz, for instance), which also matches their findings from the Phase-Amplitude Coupling presented in figures 5 and S5. In addition, it would be most relevant, also based on the results from Figure S5, for the authors to perform a similar analysis on these lower frequency 3-5Hz rhythms.\n\n2. The key result of the manuscript is that activity is feedback dominated right after stimulus onset, and feedforward right after stimulus offset. This result seems counter-intuitive and potentially in conflict with previous studies that characterise this (Samedo et al, Nat. Comm., 2022 and others). Below are a few thoughts regarding this.\n\na. Why counter-intuitive: When a visual stimulus turns on, information regarding it needs to move forward through the hierarchy, so different parts of the system get the information. Therefore, the finding that information is feedback after stimulus onset is counter-intuitive.\n\nb. Is the result above only related to theta waves, and instead spiking activity is feedforward as previously observed (e.g. Samedo et al, 2022)? It would be useful to understand this as it would increase the potential impact of the study.\n\nc. If it is all related to the task and different in the task unrelated viewing condition, it would be good to highlight these not just at the start, but across the various results of the manuscript.\n\n3. The general assumption across the results is that deep to superficial direction is feedback and vice versa. However, L4 / L6 are where thalamic inputs are coming into visual areas (particularly VISp and VISl) and feedforward is presumably seeded from here. How do the results of the manuscript compare with this anatomical connectivity? Perhaps a more general and true statement would be to call it deep to superficial, rather than calling it feedback in the results. Of course, this can be discussed as potentially feedback in the discussion part of the manuscript.\n\n4. Is this a very specific task related feature that is being discovered?\n\na. Because it is a change detection task, it would be useful to look at stimulus events, repetition events, and those that the change occurred and show the results separately.\n\nb. Is the feedback response after stimulus onset actually related to information from the previous stimulus – given it is a change detection task? Some of the Allen change detection tasks included omission trials. It would be really interesting to characterise those trials. Presumably, the feedback signals will persist at the time of expected stimulus, and the feedforward would be absent. This could be a very interesting and significant result.\n\nc. Related to the above, it has been suggested (Jordan and Keller, Neuron 2020) that feedback signals are actually arising in L2/3. Again this is contrary to the terminology loosely used in the manuscript.\n\n5. Is the offset result particular to the timing of the stimulus? Is it a delayed response to stimulus onset? For example, if a longer stimulus is tested (presumably there are different stimulus lengths in the dataset), is the feedback dominant period locked to a fixed time after stimulus onset or locked to stimulus offset? This would be important to check.\n\n6. For analysing the directionality between areas, the analyses are performed across the same cortical depth between areas. Such an analysis would be appropriate if the connectivity between areas is ordered in the similar way, where L4 is communicating to L4, L5 to L5 and so on. However, the actual biological connections are often interlaminar. Is the analysis then valid?\n\na. Why would top down/feedback information be processed in L6? If anything, L2-3 would be more likely to integrate feedback information from higher areas?\n\n7. Would be good to identify layers based on CSD. All the inferences are based on deep to superficial being L5 to L2/3 and it would be important to ensure L4 and L5 are not mislabelled. While the new methods do appear to be more reliable, it would give confidence to readers more familiar with CSD to have those methods to confirm the identification of layers.\n\n8. In terms of explaining the results a few descriptive illustrations would really help. The first part of figure 1 was very useful in explaining the concepts for the reader to follow the first few figures. Similar to this it would be beneficial to have:\n\na. a figure describing the inter-area analysis.\n\nb. a summary figure to explain all the results and their temporal order.\n\nMinor Revisions\n\nLine 87: We find that θ activity is organized as large-scale traveling waves propagating between laminae and across cortical areas.\nFor the reader that is not familiar with such characterisations, it would be useful to clarify the evidence for such statements.\n\nFigure 2, has 5 different color schemes. Panel A, C-G have colors in areas; B color on depth, E on exponent or intercept, and H color on freq band + an additional bands along the right edge. It is all too much and would be better to either split the figures or reduce the different color schema\n\nExample panels of the theta-gamma interactions are based on results from a higher visual area (VISL). I think the main text should be revisited to reinforce this, or the main figures should be changed to include a more general view of the primary visual cortex as well.\n\nLines 144-147, the higher residual power of theta rhythms in the deeper layers of Visual Cortex may also reflect passive volume conduction rather than an active process of top-down information. Can a hippocampal travelling wave be discounted?\n\nThe results of the theta wave propagation are interesting (Fig 3), but the title is misleading. This seems to be true only in L, not in other areas. Additionally I wonder if more superficial waves are being dominated by delta/alpha rhythms, while the deeper layers are more influenced by theta rhythms? Are these patterns present also in the spontaneous stimulus?\n\nLines 214-216: Can this be a result of VISp having higher amplitude responses?\n\nLines 246-247, is this wrong, it seems like its strongest in the deep to superficial direction?\n\nThe significance of the Aθ parameter presented in figure 4B is unclear. The methods state ‘A high Aθ indicates large, consistently oriented phase delays along the visual hierarchy, whereas a small Aθ indicates that phase delays are either small (meaning θ is more synchronous over the regions) or inconsistently oriented’ (Lines 644-646). So based on the figure, dark red means highly synchronous theta, while blue means consistently delayed?\n\nThe PAC analyses of all the areas and layers are really interesting. But the figures in the main text are for the spontaneous stimulus, and not the task-one that was presented in previous figures. In fact the task PAC shows little coupling between theta and gamma in VISL.\n\nHow can this difference be explained?\n\nIn the discussion (lines 470-472), the authors mention that theta was likened to alpha rhythms in primates. This is not true. In both cited papers, it was the 3-5Hz rhythm in mice that was suggested to be comparable to alpha, not theta. Please revise your main text accordingly.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "In this manuscript, the authors have analyzed the publicly available data (Allen Institute) to examine the spatiotemporal dynamics of LFP across cortical layers and visual cortical areas, while head-fixed mice were presented with flash, gray screen, and natural images (change-detection tasks). They have investigated the “traveling wave” property of theta (3-10 Hz) oscillations across cortical layers in different cortical areas, and some properties are correlated with the hierarchy scores of visual cortical areas. They have also detected gamma (30-100Hz) oscillation packets and how these are coordinated with 3-10 Hz oscillations in different cortical layers in different visual cortical areas. Overall, these analyses are great additions to the community's knowledge.\n\nIn the meantime, I have a few concerns related to the interpretation of the data.\n\n(1) The authors interpret the theta wave as signaling the feedback inputs, citing stronger theta power in deeper layers that are relevant to the feedback connectivity. However, this analysis is mostly based on the non-task-engaging experiments (flash screen and grey screen experiments, Fig2 and FigS1), not so much on the task-engaging experiments (natural image, FigS2). Thus, “feedback” interpretation during non-task-engaging experiments could not be a straightforward interpretation. Also, although the theta power is stronger in deep layers, time scale hierarchy correlations are stronger in superficial layers (Fig2H, FigS1H, FigS2H). How would this align with the authors’ interpretation of theta as “feedback”?\n\n(2) The authors interpret the gamma packet as signaling the feedforward inputs, citing stronger gamma power in the superficial layers is relevant to the feedforward connectivity. However, these analyses are again mostly based on the non-task-engaging experiments. Indeed, in the natural image task experiments, the gamma-hierachy score correlations are high in deep layers. Thus, further discussion on how gamma could correspond to the feedforward processing would be preferable.\n\nI also have a few minor technical comments/questions.\n\n(3) The authors have used “wavenumber” to evaluate the property of theta as traveling waves, claiming that it can capture single-trial measurement of LFP propagation. However, major analyses are done using the average of it, not so much different from the conventional CSD analysis.\n\n(4) Related to (3), are they really the translaminar traveling wave? Perhaps just interpreting the LFP fluctuations made by current sinks and sources as “traveling”. You can always assign “theta phase” after band pass filtering of the original LFP signals in the theta band, and thus can compute the wavenumber out of it."}, "author_response": "Authors’ Response to Reviews of\nNested spatiotemporal theta–gamma waves organize hierarchical\nvisual processing\nBrendanHarris,PulinGong\nNatureCommunications\nR:Reviewers’Comment, A:Authors’Comment, (cid:3)ManuscriptText\nWethankthereviewersfortheirinterestinourworkandthoughtfulevaluationofourmanuscript. Weprovidea\npoint-by-pointresponsetoeachcommentraisedbythereviewersbelow. Throughouttheresponse,wehave\nquotednewsectionsfromtherevisedmanuscriptwithchangedcontentindicatedingreen. Thelinenumbers\nprovidedinthisresponseletterrefertotheversionofthemanuscriptwithtrackedchanges.\nAsrequestedbytheeditor,wehaveupdatedthe‘DataAvailability’and‘CodeAvailability’sectionsofthemain\nmanuscripttoincludelinkstopubliclyhosteddataandcoderepositories. Alldisplayitemspresentedinthe\nmainmanuscriptandsupplementaryinformationcanbereproducedfromrawdataandcodethatarehostedon\npublicrepositories[1,2]andviatheCodeOceanplatform.\nDataavailability\nReduceddataunderlyingFigs.1to6andFigs.1to3andS5toS12,includingpowerspectra,wave-based\nstatistics,andcross-scalecouplingmeasures,arepubliclyavailableonFigshare[1].Thesefigurescan\nalsobereproducedonlinewiththepubliccodecapsuleat[xxxcodeoceanxxx].Thefullsourcedata\nunderlyingallanalyseshavebeendepositedontheScienceDataBank,at\nhttps://www.scidb.cn/en/s/2Unyyy[3],andourcompleteanalysisfromtheraw‘Allen\nNeuropixels–VisualBehavior’dataset[4]canbereproducedwiththeaccompanyingcode\nrepository[2].\nLines1033-1043\nCodeavailability\nTheanalysesandfigurespresentedinthisworkarefullyandfreelyreproducibleusingtheaccompanying\ncodepackage,SpatiotemporalMotifs.jl[2].\nAllanalyseswereperformedwithv1.10.10oftheJuliaprogramminglanguage[5].We retrievedthe\n‘AllenNeuropixels–VisualBehavior’dataset(atmanifestversionv0.5.0)withacustomJuliawrapperfor\ntheallensdkPythonpackage,AllenSDK.jl.AllenNeuropixels.jl[6]wasusedtoformat\nthedata,alignLFPtostimulustimes,infercorticaldepths,andproducevisualizationsofthedataset(such\nasVideoS1).WeusedtheTimeseriesTools.jlpackage[7]forpowerspectraandwave-based\nanalyses, andNormalization.jl[8]fordatanormalization.Finally,weproducedfigureswith\nMakie.jl[9],estimatedconfidenceintervalswithBootstrap.jl[10],andstandardizedour\nsoftwarewithDrWatson.jl[11].\nLines1044-1056\n1\n1. Reviewer#1\nR1.1: Inthismanuscript,HarrisandGonganalyzedhigh-densityelectrophysiologicalrecordingsfromthemouse\nvisualcortex. Theyproposedthat,duringstimulation,across-scaleoscillatorydynamicemerges,involving\nnested theta-gamma waves that are hierarchically organized and run alongside spiking activity. After\ndisentanglingtheLocalFieldPotential(LFP)signalintotheta,gamma,andthe1/fcomponents,theyfirst\nshowthattherearetwodistinctthetaoscillationcomponents: onethatisoffset-related,propagatingfrom\nlayer2/3todeeperlayers,andanotherthatisonset-related,reflectingfeedbacksignaling. Theoffset-related\nthetahasahighamplitudewhentheanimalisdisengagedfromthetask,whereastheonset-relatedtheta\nappearsonlyduringactivetasks. Thistask-evokedthetaexhibitshierarchicalorganizationandproduces\nspatiotemporal dynamics compatible with previously described traveling waves. Gamma oscillations are\norganizedinpacketsnestedwithintheta,structuredacrossbothcorticalareasandlayers. Thisnestedactivity\ncorrelates with spiking activity during stimulation, suggesting a flexible, spatiotemporal mechanism for\ninformationmultiplexingandmodulation. Themanuscriptanalyzesahigh-qualitydatasetandpresents\nintriguingresultsandideas. Itprovidesadetaileddescriptionofthetadynamicsacrossthemousevisual\ncortex,whichalignswellwithrecentviewsonthetaandgammaoscillationdynamics. However,itlacksa\nmechanisticcharacterizationofhowthesenestedoscillationscontributetoinformationprocessingbeyond\ntheiralignmentwithspikingactivity. Thislimitationweakensthemanuscript’scentralmessage.\nA1.1: We appreciate the reviewer’s positive assessment of our findings and constructive comments regarding the\nmechanisticinterpretationofourfindings. Assuggestedbythereviewer,wehavenowupdatedourmanuscriptto\nincludeamoredetailedinterpretationoftheroleofθwavesinvisualprocessing,backedbynewresultsshowing\ntask-relatedmodulationofneuronalspikingrelatedtopropagatingθ. Weintroduceeachofthesechangesin\ndetailinthefollowingresponsepoints.\nR1.2: Thefollowingsuggestionsaimtostrengthenthemanuscript’simpact: Theprimaryconcernisthat,while\ntheauthorsdemonstratesignificantLFP-spikecouplingatbothfrequenciesofinterest,itremainsunclear\nwhether this coupling influences perceptual information processing or transmission. For example, the\nauthors correctly perform a PPC analysis on neighboring LFP channels relative to the spike cluster of\ninterest. However,thefunctionalrelevanceofthiscouplingremainsuncertain. Arespike-fieldlockedclusters\nfunctionallymeaningful? Canperceptualinformationbeaccuratelydecodedfromspike-LFPcoupling? Does\nthiscouplingfacilitateinformationtransferacrossdifferentcorticalareas? Toenhancethemanuscript’s\nclaims, the authors could quantify a potential causal relationship between LFP oscillations and spiking\nactivity. ThisanalysiscouldcomplementtheirresultsinFigure3Fonhit/missdecodingofthetawavesand\nofferamoredetailedcharacterizationoftheputativecausalityunderlyingtheirfindings.\nA1.2: Assuggestedbythereviewer,wehavenowperformedadditionalanalysistoshowthattheinterplaybetweenθ\nLFPandneuronalspikingisrelevanttotaskperformance. Wefoundthatbothpreferredθphasesandneuronal\nfiringratescandistinguishhitfrommisstrials,withpreferredphasesmovingtowardθtroughsduring‘hit’trials.\nThesefindingsareconsistentwithacausalrelationshipwherebytop-downθenhancessuperficialspikingto\nraisesensitivitytostimuluschange. WenowpresentthesenewresultsinFigs.6dto6g,referredtointhemain\ntextquotedbelow. Wehavealsoaddedadiscussionofthepotentialmechanisticroleofθ–spikecouplingat\nlines615-636andlines637-659oftheDiscussion.\nNestedθandγreflectunderlyingpatternsofneuronalspiking\nFinally,weaskedhowtheθandγoscillationarecoupledtounderlyingspikingactivity.Forθ,we\ncomputedaspike–phasecouplingstatistic—thepairwisephaseconsistency[12]—betweenspikesfrom\neachunit(putativeindividualneurons)andtheθLFPphasefromtheclosestLFPchannel.Figure6a\nshowsthetheaveragepairwisephaseconsistency(PPC)duringtheNatural_Imagesstimulus,showing\nstrongspike–phasecouplingin thesupragranularlayersoflower-orderareas.Asimilartrendis\nreflectedinspike–amplitudecouplingbetweenγamplitudesandneuronalspiking,quantifiedasthe\n2\nmeannormalizedr acrossspiketimes.AsplottedinFig.6b,thedistributionofspike–amplitude\nγ\ncouplingisstrongestinthesupragranularlayers,asforspike–phasecoupling.Together,thesehigh\nvaluesofspike–phaseandspike–amplitudecouplingcorroborateourearlierfindingsofnestedLFP\noscillations,showingthatthesimultaneouscross-scalecoordinationofθandγextendstothelevelof\npreciseneuronalspiking.\nToelucidatethemechanisticrelationshipbetweenspatiotemporalθcomponentsandneuronalspiking,\nwecalculatedtheaveragepreferredθphasesacrossvisualcorticalareasandlayers,plottedinFig.6c.We\nfoundthatneuronsinsupragranularlayersprefertroughsoftheθcycle,whereasneuronsininfragranular\nlayerspreferpeaks;thistrendisrepeatedacrossareasinthevisualhierarchy,andbroadlymatchesthe\npreferredphasesofγpackets(cf.Fig.5).Dothesepreferredphasesreflectafunctionalroleoftraveling\nθcomponents?Weselectedneuronsthatweresensitivetoθ,withapairwise-phase-consistencyof\ngreaterthan0.1acrosstheentireNatural_Imagesstimulusepoch,thenaskedhowlayer-wise\ndifferencesinpreferredθphasesandfiringratechangebetween‘hit’and‘miss’trials.Wefoundthat\nunitsinsupragranularlayersofthefourlowest-orderareasshowedasignificanttendencytopreferθ\ntroughsduring‘hit’trials,asshowninFig.6d.Bycontrast,deeperlayersshowednosignificant\nchangeinpreferredphases.\nIfneuronalspikinginsupragranularlayersbecomesentrainedtothetroughofθduringtrials\nwhereθhasthestrongestdeep-to-superficialpropagation(seeFig.3),isthereacorrespondingchange\ninfiringratethatwouldsuggestthepropagatingθtroughislinkedtoenhancedneuronalactivity?To\naddressthisquestion,wecalculatedthemeanfiringrateofθ-sensitiveneurons,groupedinto\ninfragranular,granularandsupragranularcompartments,showninFig.6efor‘hit’trialsandFig.6f\nfor‘miss’trials.Wefoundthatsupragranularunitshaveasignificantlyenhancedpeakfiringrate\nduring‘hit’trialsforthefourlowest-orderregions(VISp–VISal),occurringat∼80msafterstimulus\nonset.Collectively,thesefindingssupportthespatiotemporalθ–γschemesummarizedinFig.1a:\nduringsuccessfulchangedetection,aθtroughcarryingtop-downinformationtravelsupwardfrom\ndeeplayersandcoincideswithenhancedbottom-upγpacketsandneuronalspikinginsuperficial\nlayers.Thesetask-relatedspike–phaserelationshipscometogethertoformthepictureofnested\ndynamicspresentedinFig.1a:thatatravelingθtroughmovesupwardfromdeeplayersduring\nsuccessfulchange-detectiontrials,whereitcoincideswithenhancedneuronalresponsesinshallow\nlayers.\nLines494-539\nR1.3: Asecondconcernisthat,whilethemanuscriptfocusesonnestedtheta/gammaoscillations,certainadditional\nanalysesfeelsomewhatdisconnectedfromthemainnarrative. Forexample,thedecompositionoftheLFP\nsignalintoitsoscillatoryand1/fcomponentsishighlyinteresting,yetitisunclearwhetherthe1/fcomponent\nplaysanyfunctionalrolebeyonditsdescription. Theauthorscouldstrengthentheirhypothesisbycontrasting\nthefunctionalrelevanceoftheta/gammanestedactivityagainst1/factivity. Iffunctionaldifferencesbetween\nthesecomponentscanbeidentified,thiswouldfurthersupportthemanuscript’skeyclaims.\nA1.3: Toaddressthefunctionalrelevanceofthe1/fcomponent,wehavenowaddedananalysisoftrial-by-trial1/f.\nThis analysis, presented in Fig. 1 and quoted below, indicates that there are significant changes in the 1/f\nexponentbetween‘hit’and‘miss’trials,butthemagnitudeofthischangeismuchsmallerthanthevariationin\n1/f exponentacrosscorticallayers. Thisnewresultsuggeststhatdynamicchangesinthe1/f componentare\nrelevant,butlessfunctionallymeaningfulthanevokedθorγ components;wehavenowaddedthispointtothe\nResults,atlines180-196.\n3\n-0.05 0 0.05 0.1\n1/f exponent (hit - miss)\n)%\n(htped\nlacitroC\n1/f contrast\n0\n50\nycarucca\ndecnalaB\n1/f classification\n1\n0.75\n0.5\n-0.2 0 0.2\nWeight\n)%\n(htped\nlacitroC\na b cClassifier weights\n0.0\nArea\nVISp\nVISl\n0.5 VISrl\nVISal\nVISpm\nVISam\n1.0\nFigure1:Smallvariationsin1/fexponentreflecthit/missperformance.\naDifferenceinthemedian1/fexponentbetween‘hit’and‘miss’trials,plottedacrosscorticaldepth\nforeachvisualarea.1/f exponentswereestimatedbetween10-300HzfortheLFPbetween\n-0.25-0.75srelativetotrialonsets.Bandsindicatethe95%confidenceintervals(CI),andfilledcircles\ndenoteasignificantdifference(p<0.01,Mann-WhitneyUtestagainsttrial-shufflednullswith\nBenjamini-Hochbergcorrection).Theoveralldifferencein1/f exponentbetween‘hit’and‘miss’is\nsmall,beinglargestinVISpatapproximately0.05.\nbBalancedaccuracyforalineardiscriminantanalysisclassifiertrainedtodistinguish‘hit’from\n‘miss’trialsusingeitherregionally-averaged1/fexponents(‘Regional’,left,yellow)oralllayer-wise\nexponents(‘Layerwise’,right,green).Coloredboxesshowtrueperformanceacrosssessions,while\ngrayboxesshowperformancefortrial-shufflednulls.Boxplotsshowtheinterquartilerange(IQR;\nboxes),median(line),outliers(dots),andmostextremepointswithin1.5×IQRofeachquartile\n(whiskers).\ncMediannormalizedweightsfromthelayer-wiseclassifieracrosssessions.\nR1.4: Athirdconcernrelatestotheoriginanddynamicsofthetaactivity. Thesupplementaryvideosuggeststhatthe\ntravelingoffset-relatedthetawavemayhaveasubcortical(hippocampal?) origin. Ifso,someoftheobserved\nthetaactivitymaynotbeintrinsictothevisualcortexbutinsteadvolume-conductedfromsubcorticalregions.\nTheauthorsshouldclarifywhetherthisthetaactivitytrulypropagateswithinthecortexorismerelyaneffect\nofvolumeconduction. Althoughcompletelyrulingoutvolumeconductionisdifficult, theauthorscould\nquantifyphasedifferencesacrosscorticalareastocheckforzero-phaserelationships,whichwouldsupport\nthehypothesisofvolumeconduction.\nA1.4: Assuggestedbythereview, wehavenowexpandedourdiscussionofvolumeconduction. Specifically, we\nhighlightthreekeysignsthatθpropagationisacorticalphenomenonratherthananeffectofvolumeconduction\nfromthehippocampus. First,ourfindingthatdeep-to-superficialθisstrongestinlower-orderregionsofthe\nhierarchy,diminishinginhigher-orderregions. Ifθwerevolumeconducted,emanatingbroadlyfromthedeeper\nhippocampus,wewouldexpectaroughlyequivalentstrengthacrossallregions(Fig.3). Second,ouranalysisof\nθphasedelaysbetweencorticalregionsshowsastronghierarchicaldirectionality,whichisinconsistentwiththe\nzero-phaserelationshipexpectedfromvolumeconduction(Fig.4). Finally,bothtranslaminarandhierarchical\npropagationreversedirectionalityatstimulusoffset,furtherindicatingthatθpropagationisindependentofa\nvolume-conductedsignalfromthehippocampus. Nevertheless,weagreethatvolumeconductionisdifficultto\ncompletelyruleoutasasourceofdeep-to-superficialθpropagation. Wehaveincorporatedthesepointsintothe\nDiscussionatlines679-695.\nR1.5: Finally,whilethediscussionsectioniswell-informed,itcouldreachabroaderaudiencebyprovidingamore\ncomprehensive conceptual framework for the manuscript’s findings. For example, the authors correctly\n4\nciteLisman’sworkonnestedtheta/gammaoscillations,buttheyshouldelaborateonthecurrentstateof\nthistheoreticalframework. Whatrecentfindingssupportorchallengethisperspective? Additionally,the\ndiscussionshouldaddresswell-knownlimitationsofcross-frequencycouplinganalysesandclarifyhowthis\nstudy mitigates these limitations. Remarks on code availability: The code contains a README file that\nsatisfactorilyexplainstheinstallationandreproductionoftheanalysisofthismanuscript\nA1.5: Assuggestedbythereviewer,wehavenowextendedourdiscussionoftheθ–γ codetohighlightthelimitations\nofexistingmethodsforquantifyingcross-frequencycoupling,suchassensitivitytofilterartifactsandnonlinear-\nity[13,14]. IntheDiscussion,wenowhighlightthatweusedbroadfrequencybands,single-trialwave-based\nmethods,andmultiplestimulusconditionstoaccountfortheselimitationsandseparatetask-relateddynamics\nfrombackgroundactivity. Wealsoemphasizethekeyadvanceofourθ–γ framework,whichisthatitaccounts\nfordistributedcross-scaledynamicsratherthanlocalcouplingorpairedlong-rangeinteractions.\nOurworkalsoputsforwardmethodologicaladvancesforstudyingtheθ–γcodeinspatiotemporal\nsettings.Conventionaltoolsforanalyzingcross-frequencycoupling[15,16]typicallyrequire\nband-passfilteringofbroad-bandneuralsignals;filteringcanintroducenon-causalartifactsand\ngenerateambiguitybetweennestedlinearoscillationsandnonlinearwaveforms,suchassharp\ntransients[14].Wemitigatedtheselimitationsbyselectingrelativelybroadfrequencybandsforθ\nandγ,therebyreducingpotentialfilterartifactsandcontrollingfornegativeinstantaneous\nfrequencies[13].Wealsoadoptedanon-stationary,wave-basedapproachthatquantifiedsingle-trial\ndynamicsratherthantrial-averagedphenomena,andstudiedvariousactive,passive,and\nspontaneousstimulusconditionstodisentangletask-relatedeffectsfromongoingbackground\ndynamics.Importantly,theslowerθcomponentwehaveidentifiedshouldnotbeinterpretedasa\npersistent,narrowbandoscillation.Wefoundbriefboutsofoscillatoryθfollowingstimulusoffset\n(seeFig.2),butduringtheshort250msstimuluspresentation,θisbetterviewedastheslow\ncomponentoftheevokedresponse.Nevertheless,θduringthestimulusperiodstillexhibitsreliable\nphasecouplingwithγpacketsandneuronalspiking,demonstratingthatnon-stationarycross-scale\ninteractions—ratherthanstrictlystationaryornarrow-bandoscillations—arecriticalcomponentsof\nvisualprocessingandvisualbehavior.\nLines662-678\n5\n2. Reviewer#2\nR2.1: Thepaperanalyzeselectrophysiologicaldatafrommousevisualcortexduringavisualtask,withorwithout\nengagingattention. Theauthorsshowtheta-bandtravelingwavespropagatingatdifferentscalesbetween\nregionsandcorticallayers. Theyalsorelatethewave’sphasetothegamma-bandactivityandspikes. Overall,\nthepaperismethodologicallysound,andtheresultsareinteresting. I’mnotentirelysureaboutthenovelty\nregardingthepreviousliteratureaboutgamma-thetacoupling,butthespecificitymayapplytothespecific\nanimalmodelused. Belowaresomepointsthatmayneedfurtherclarificationforabroadaudience.\nA2.1: Wethankthereviewerfortheirconstructivesuggestions. Weappreciatetherecognitionofthestudy’smethod-\nologicalrigorandtheinterestinourfindings. Weprovideapoint-by-pointresponsetoeachcommentraisedby\nthereviewerbelow.\nR2.2: How does the intra laminar theta relate to between areas theta? More precisely, I don’t understand the\nassumptionthatisusedtocalculatetheorderparameter,i.e.,thattheaveragetranslaminarthetapropagation\ninoneregioninformsaboutthethetapropagationbetweenareas.\nA2.2: Intherevisedmanuscript,wehavenowclarifiedthepotentialrelationshipbetweentranslaminarθandinter-areal\nθintheResults,asquotedbelow. Specifically,weindependentlyanalyzedpropagationbetweencorticallayers\nandcorticalregions,andouranalysisdidnotinitiallyassumearelationshipbetweentranslaminarandinter-areal\nθ. Our key finding was that translaminar θ propagation matches inter-areal θ propagation in both timing\n(onset,offset)anddirectionality(superficial–deepordeep–superficial,andfeedbackorfeedforward). Toclearly\nillustratethesepoints,wehaveaddedanewschematictoFigs.4eand4f. Finally,wehavealsoclarifiedour\ninterpretationofinter-arealθpropagation,asmeasuredbythehierarchicalorderparameter,atlines363-379\noftheResults.\nTogether,thesecomplementaryfindingsindicatethatθ-scaleactivityplaysadualroleinvisual\nprocessing,flexiblyswitchingbetweenfeedback-oriented(top-down)andfeedforward-oriented\n(bottom-up)modesaccordingtothevisualstimulusandtaskdemands.Figures4eand4fsummarizehow\ntranslaminarandinter-arealθmayswitchbetweenthesetwoprocessingmodesduringtasktrials.\nFirst,asdepictedinFig.4e,wesuggestthatstimulusonsettriggersarapidtransmissionof\nfine-grainedinformationaboutthevisualscene[17]alongafeedforwardpathway.Thisfeedforward\npassappearstodominateγLFPactivityinlayers2/3and4(seeFig.2),withtheearliestsignaturesof\nvisualinputappearinginVISplessthan50msafterstimulusonset[18].Ashorttimelater,\nbetween\n50–150ms,afeedback-orientedθsignalarrivesinsuperficiallayersofthelower-ordervisualareas,\nbeingtransmittedfromhigher-orderregionsalongadeeperpathway[19,20].Thisfeedbackθ\ncarriestop-downinformation,likelytiedtotaskengagement(seeFig.S3),thattravelsbetweenlayers\nineachregiontomeetandmodulatefeedforwardγinthesuperficiallayers.Atstimulusoffset,\nhowever,θplaysadifferentrole,asshowninFig.4f.Offset-evokedθoscillationsmaysynchronizethe\nfeedforwardγpacketsthattravelfromlower-orderregionstohigher-orderregionsinasuperficial\npathway.Theseθoscillations,whichhaveacharacteristicfrequencyof∼4HzinVISp(asshown\ninFig.2),resembletheevokedactivitypreviouslyobservedinlowervisualcortical\nareas[21,22,23,24].Here,wehaveshownhowthisoffset-evokedθoscillationpropagatesacrossthe\nvisualhierarchyandcorticallayers,whichmayservetoperiodicallyshieldneuralactivityinhigher\nregionsfromthefeedforwardsignal,reducingtheinterferenceofbottom-upinformationwith\ncognitiveprocessing.Thesefindingssuggestthatthespatiotemporalstructureofθaffordsitaflexible,\ncontext-dependentfunctionalrolethatadaptstothecomputationaldemandsofvisualtasks.\nLines396-417\nR2.3: Infigure1A,therelationshipbetweenpacketsIandII,andbetweenpacketsIIIandIVisunclear. Clearly,\ntheyareallmodulatedbytheslowthetawave,butit’sunclearhowtheyarerelatedtoeachother. Doedpacket\nI(III)causepacketII(IV)?IsthefactthatIIandIVareenhancedandinhibitedunrelatedtothepresence\nofpacketsIandIII?\nA2.3: We have now improved our description of spatiotemporal θ–γ code illustration presented in Fig. 1a. By\nperformingnewanalysesontherelationshipbetweenθphasesandspikingactivity,wefoundthatθtroughs\ncorrespondwithincreasedneuronalexcitabilityinsuperficiallayers(seeFig.6);thisisnowreflectedinthe\nschematicfigure. Wehavealsoclarified(atlines142-161)thatpacketsIandIIIarenotnecessarilyrelatedin\nthisschematic,asitdepictstwocontrastingscenarios(thepresenceandabsenceofθpropagation).\nR2.4: TheresidualpowerinFigure2Fseemssimilaracrosscorticalregions(exceptforVISam). Isn’tthisthe\noppositeofwhatonewouldexpectconsideringtheproposedfeedback/top-downrole?\nA2.4: Wethankthereviewerforraisingthispointofpotentialconfusion,whichhighlightsthereasoningbehindour\nlaternon-stationarywave-basedanalyses.Intherevisedmanuscript,wehavenowclarifiedthatthedistributionof\nresidualpowerforθisgenerallyconsistentwithatop-downrole,butthatinferringthefunctionalroleofθfrom\nstationaryspectralanalyseshassignificantlimitations. AsshowninFig.2h,andreproducedforotherstimuliin\nFig.S5handFig.S6h,theresidualθpowerincreasesrobustlyalongthecorticalhierarchy(withhighercortical\nregionshavingagreaterresidualθpower). Whilethisresultisconsistentwiththeinterpretationthatθhasa\nfeedbackrole(assuggestedbymanyearlierstudies),itiscrucialtonotethatthisspectralanalysisaggregates\novertheentirestimulusepoch.Inothercases,especiallyfollowingtheoffsetofimagepresentations,θcanappear\ntodecreaseinstrengthalongthehierarchy(seeFig.S2b). Henceweadoptawave-basednon-stationaryanalysis\nin the remainder of the paper, which reveals θ has a dual functional role and follows opposing hierarchical\ngradientsatdifferentperiodsaroundthevisualstimulus. WehavenowincludedthesepointsintheResults,at\nlines213-242.\nR2.5: Itisunclearhowtheanatomicalhierarchyusedtoobtainthecorrelationsisdefined,andhowtheresultscan\nbeinterpreted(figure2Eand2H)otherthanshowinga(well-known)dynamicalcorticalgradient.\nA2.5: We have now clarified how the anatomical hierarchy scores used in Fig. 2e and Fig. 2h are defined, at\nlines 197-212 of the Results. In summary, they were calculated by Harris et al. [25] and Siegle et al.\n[18]basedontheanatomicalorfunctionalconnectivitybetweenvisualareas,providingarobusthierarchical\norderingbasedonfeedforwardandfeedbackneuronalprojections. Toourknowledge,nopriorstudieshave\nestablishedthetherelationshipbetweenθ,γ,or1/f componentsoftheLFPsignalacrossthevisualcorticalhi-\nerarchyinmice. Inparticular,wefindstronghierarchicalgradientsandnon-trivialtask-relatedeffects,including\ntask-dependentgradientsofγ and1/f exponentbuttask-independentgradientsofθand1/f intercept. These\nfindingsgobeyondrecentworksthathavefocusedontheprimaryvisualarea[26]oraggregatedovermany\ncorticalregions[27]—ineithercase,missingdetailedhierarchicalvariations. Whilethesestationaryspectral\nfeaturesarenotthemainfocusofourwork,theyareacrucialmotivationforourdevelopmentofwave-based\nmethodstoanalyzethehighlynon-stationary,task-dependentpropertiesoftheθandγ components. Wehave\nnowstrengthenedthismotivationatlines213-242intheResults.\nR2.6: Fromaninterpretationalpointofview,whatdoesitmeanthatthetapropagatesfromdeepertosuperficial\nlayers? Shouldn’toneexpecttheopposite,consideringthatthe(thetaactivityinthe)deeperlayerencodes\nexpectationsratherthansensoryinformation? Wouldthisreflectmostlytop-downmodulationratherthan\nsensory processing? Could the author decode some stimulus information from the theta activity and its\npropagation?\nA2.6: Toaddressthepotentialpointofconfusionraisedbythereviewer,wehavenowaddedanewschematicfigure\n7\ntoFig.4. Withthisfigure,wesummarizeourfindingsinrelationtoθ propagationandclarifytheirpossible\nfunctionalroleinrelationshiptotop-downandbottom-upmodulationofγ activity. Particularly,wefoundthatθ\npropagatingfromdeeptosuperficiallayersinthe50–150mswindowfollowingstimulusonsetappearstocarry\ntop-downinformationrelatedtoexpectationorengagement,andthatthissignalisstrongestwhenthemouse\ncorrectlyidentifieschangesinstimuluscontent(seeFig.3). Indeed,aspointedoutbythereviewer,wecould\ndecodehitvs. misstrialsusingtheorderparameter(asummarystatisticofθpropagation;seeFig.3g). Wealso\nfoundstrongdifferencesinθpropagationfordifferentstimuli(suchaspassiveflashstimuli),assummarizedin\ntheDiscussionsectionquotedbelow. Tostrengthenourfindingthatonsetθcarriestop-downinformation,we\nhavealsoaddedacomparisonofθpropagationbetweenactiveandpassiveepochswithnaturalimagestimuli\n(seeFig.S3),showingthatθpropagationissignificantlystrongerinactiveepochs,eventhoughthestimulus\nsetisidentical. Together,ourfindingsondeep-to-superficialtop-downθ(whichoccursafterstimulusonset)\narecontrastedbyasecondsetoffindingsonsuperficial-to-deepbottom-upθ,whichweobserveafterstimulus\noffset(seeA2.7).\nSeveralkeyobservationsindicatethatfeedback-orientedθbetween50–150msfollowingstimulus\nonsetrepresentsanengagement-relatedtop-down‘priming’offine-grainedγandspiking\nactivity[28,29].First,translaminaronsetθwavestravelfromdeeplayerstosuperficiallayers\n(seeFig.3a);deeplayersareoftenthoughttobefeedbackdominated[28].Thistranslaminar\npropagationisstrongestinlower-orderareas(Fig.3d),butθalsotravelsacrossthehierarchy(Fig.4b)\nfromhigher-ordertolower-ordervisualareas;higher-orderregionsareknowntomodulatestimulus\nresponsivenessandguidecomplexaspectsofvisualcognitionsuchasdecisionmaking[30].Moreover,\nFig.5fandFig.S10showthatnestedθ–γinteractionsatdifferenttimesarestrengthenedindifferent\nlayers,whichalignwithknowncortico-corticalpathways.Inparticular,onsetθwavesarecoupled\nmoststronglytoγinlayer1oflower-orderregions,whichisakeysitefortheintegrationoftop-down\nandbottom-upsignals[31,32].Bycontrast,onsetθislockedtoγinlayer6ofhigher-orderareas,\nwhichisassociatedwithastrongcortico-corticalfeedbackpathway[33].Thisalignmentwith\nfeedbackanatomicalfeaturesisreflectedinthetop-downfunctionalroleofθ;deep-to-superficialθ\nwavesaremostcoherentwhenmicesuccessfullyrespondtoastimuluschange,andareabsentattheonset\nofpassiveflashstimuli(withoutanytaskdemands),suggestingthattheyarespecificallyassociatedwith\nactiveengagement(Fig.3e).Finally,neuronalspikepatternssuggestthatdeep-to-superficialθ\npropagationcorrespondstoincreasedexcitabilityofsuperficialneuronsafterstimulusonset;firing\nratesareenhancednearθtroughsinsupragranularlayersduring‘hit’trials(Fig.6).\nLines615-636\nR2.7: Idon’tthinkthereisevidencetosupporttheauthors’proposal(line262)thatthetaprovidesfeedforward\nprocessing: shouldn’toneobserveaforwardpassalsoattheonsetofthestimulus,then?\nA2.7: Wehavenowexpandedourdiscussionofthedifferentrolesofθandγ intheResults(lines396-417),andin\ntheDiscussion(asquotedbelow). Insummary,thefeedforwardpassatstimulusonsetislikelycomposedof\nrapidneuronalspikingandγ events[17]thatreachthefirstvisualareawithinthefirst∼50msafterstimulus\nonset. We find that θ, as the slow component of the LFP response, is most functionally relevant ∼ 100ms\nafterstimulusonsetinsuperficiallayers,andhasmanysignaturesofatop-downsignal(seeA2.6). Incontrast,\npropagating θ that arises after stimulus offset appears to have a role in bottom-up processing. Rather than\ncarryingfine-grainedinformationaboutthestimulus(theroleofγpacketsandneuronalspikes),weproposethat\nθorientedinafeedforwarddirectionactstosynchronizeγ andspikingactivityacrossthecorticalhierarchy.\nThisbottom-upinterpretationofoffsetθissupportedbyourfindingsthatoffsetθpropagationpersistswhenthe\nmouseisdisengagedfromthetask(seeFig.S3)orwhenthestimulusisasimpleflash(seeFig.3e),unlikeonset\nθ. Furthermore,theamplitudeofevokedθcomponentsisstrongerinlower-orderareas(seeFig.2),andthat\noffsetθ–γ lockingisstrongestinlayers2/3and4(inlower-orderareas).\nand,atthesametime,fromsuperficialtodeeplayers.Thesesuperficial-to-deepθwavesareevoked\nunderbothactiveandpassiveconditions(Fig.3c),aswellasphase-lockedtoγpacketsinlayers2/3and4\n(Fig.5f),whicharedominatedbyfeedforwardinput[28,33].Unliketheθcomponentpresentduring\ntheshortperiodinwhichthestimulusispresented,offsetθformsatransientoscillationthatis\nparticularlystronginlower-orderareas,withameanfrequencyof∼4Hz(seeFig.2).Suchevokedθ\noscillations[24]havepreviouslybeenassociatedwithaccumulatedexposuretoastimulus[34,35]or\nreward[21,36].Wesuggestthattheseθwavesmayfacilitatemultiplexedprocessinginhigher-order\nregionsbytransientlyinhibitingthebottom-uppassageofvisualinformation,therebyopeningperiodsfor\ninputfromcorticalregionsoutsidethevisualsystem[37].Thisinterpretationisconsistentwiththerecent\nsuggestionthat3–5Hzrhythmsinmicemaybefunctionallycognatewiththeαrhythm(approximately\n7–13Hz)inprimates[38,39].Inprimates,αhasbeenshowntomodulatehigher-ordercognitive\nfunctions[40,41]includingattention[42,43]andgoaldirection[44].Furthermore,whilepriorworkhas\ndescribednestedwavesthattravelhorizontallyacrossthemousevisualcortex[45],withafeedbackθand\nfeedforwardγ,ourworkshowsthatthesewavesarenotsimplystable,unidirectionalphenomena.Instead,\nbyconnectingtranslaminarandinter-arealθpropagationacrossthehierarchy,testingtheirdirectfunctional\nimplicationsforvisualchangedetection,andlinkingtravelingθtotask-relatedspiking,wefoundthatθ\nflexiblyswitchesdirectionalitytoimplementbothtop-downprimingandbottom-upsynchrony.\nLines637-659\nR2.8: Fromthepowerspectra,theyanalyzebroadbandgamma(30-100inequation5). Wouldtheauthorsimplythis\nisafairlylargebroadbandeffect(possiblypresentinhigh-betaaswell?),orwoulditbepossibletolocalizea\nnarrowbandgammafrequency?\nA2.8: Wethankthereviewerforraisingthispoint, andacknowledgethattheγ bandstudiedinourworkisbroad\ncomparedtopapersthatstudynarrow-bandγ oscillations(suchasShinetal.[46]). Inourworkspecifically,\nwechoseabroadγ bandtocapturethefullrangeofhigh-frequencycomponentsthatdominateallsixvisual\ncorticalareaswestudy,aswefindthatonlyVISphasastrongnarrow-bandγpeak(seeFigs.2,S5andS6). We\nnotethatthisbroadergammabandisalsocommoninliteratureonthemousevisualcortex[47,48,49,50]. We\nhavenowstrengthenedthismotivationintheResults,atlines163-179.\nR2.9: Fromthemethod,itisunclearwhytheinstantaneousangularwavenumber(i.e. thespatialfrequency)isthe\nnegativespatialphasegradient(amathematicalorevenanintuitiveexplanationwouldverymuchhelp). It’s\nalsounclearwhythisvaluewouldquantifythedirectionalityofthewave.\nA2.9: Assuggestedbythereviewer,wehavenowaddedanintuitiveexplanationofwhythewavenumberisthenegative\ngradientoftheanalyticphasetotheθtravelingwavessectionoftheMethods. Insummary,thewavenumber\nisdefinedasthenegativespatialphasegradientbecauseoftheconventionthatphaseincreasesovertimefor\nsignalswithpositiveinstantaneousfrequencies. Thispointisdetailedbelow:\nNext,tomeasureθpropagationacrosscorticallayers,weusedwave-basedmethodsthathavepreviously\nbeenappliedtohorizontalcorticalpropagationsmeasuredbyvariousrecordingmodalities[51].We\nrepresentarealspatiotemporalsignalinaconventionalcomplexformgivenby:\na=rexp[iφ]=rexp[−i(kx−ωt+b)], (1)\nwherer(x,t)istheinstantaneousamplitude,φ(x,t)isthephase,k(x,t)isthewavenumber(orangular\nspatialfrequency),ω(x,t)istheangulartemporalfrequency,andb(x)isaphaseoffset.Since\n9\nφ=−kx+ωt−b,theangulartemporalfrequencyω(x,t)isthetemporalgradientofthephase:\n∂φ\nω= . (2)\n∂t\nSimilarly,theinstantaneousangularwavenumberisthenegativespatialphasegradient:\n∂φ\nk=− . (3)\n∂x\nTogether,thewavenumberandtemporalfrequencydeterminethedirectionofwavepropagation.By\nEq.(7),phaseincreasesovertimeforsignalswithapositiveinstantaneousfrequency,butbyEq.(8),\nphasedecreasesinthepositivexdirectionforapositivewavenumber.Thereforewhenboththe\ninstantaneousfrequencyandwavenumberarepositive,agivenwavefront(suchasthecontourof\nφ=180◦,atrough)willbeginatasmallerxpositionandmovetolargerxpositionsastimepasses,\nandthewaveaswholewilltravelinthepositivexdirection.Whenappliedtotheθsignal,the\ninstantaneouswavenumberkquantifiesthelocal,instantaneousdirectionalityandcoherenceofθ\npropagation;positivekindicatesthewaveistravelingfromsuperficialtodeeplayers,whereasnegativek\nindicatesdeep-to-superficialpropagation.\nLines777-804\nR2.10: Itisalsounclearwhythe’periodsofnegativefrequency-definedasthederivativeofthephase-wererare’.\nWhatdoesthisimply? Thatthedirectionalitywasmostlyconsistentacrosstrials/electrodes? Cantheauthor\nquantifyandreporthowmanytrialsweremasked(asspecifiedinlines568-569)?\nA2.10: Atthesuggestionofthereviewer,wehavenowreportedtheproportionofnegativefrequencytimepointsfor\ntrialperiodsintheMethods,andclarifiedtheinterpretationofnegativefrequenciesinthecontextoftheanalytic\nsignal:\nWenotethatinferringthepropagationdirectionfromthewavenumberreliesontheassumption\nthattheinstantaneousfrequencyispositive;thepropagationdirection,asmeasuredbythephase\nvelocity(or,theratiooftheinstantaneousfrequencyandthewavenumber)hasanoppositesigntothe\nwavenumberwhentheinstantaneousfrequencyisnegative.Fortheanalyticsignal,basedonthe\nHilberttransform,negativefrequenciesindicateperiodswherethesignalcontainsmorethanone\nsimultaneousspectralcomponent.Forthefilteredθsignal,wefoundthatperiodsofnegative\nfrequencywererare(occurringfor<2%oftrialtimepoints)andthatdirectionalitiesinferredfrom\neitherthewavenumberorthephasevelocityagreedclosely.Sincethewavenumberismore\nnumericallystablethanthephasevelocity,andlesssensitivetonoise,wepresentresultsbasedonthe\nwavenumberinthemaintext,andmaskanynegativefrequencyperiodsfromourwave-based\nanalyses.\nLines818-828\nR2.11: Overall, it’s unclear the relationship between translaminar theta and between regions. The results are\nextremelysimilar,aspointedoutbytheauthors(lines291onward). Sincethelaminardirectionisroughly\northogonaltothecorticalsheet,onemaywonderhowandwhythetworesultsrelatesostrikingly.\nA2.11: We have now clarified the interpretation of our measurement of inter-areal θ propagation in the Results at\nlines363-379. Insummary,thehierarchicalorderparametersdonottracethedetailedstructuralprojections\nbetween regions and layers; instead, they are better viewed as tracking the sequence of neural activity in\nequivalentlayersacrossthehierarchy,whicharethoughttohavesimilarfunctionalroles. HavingaddedFigs.4e\nand4ftosummarizetherelationshipbetweentranslaminarandinter-arealθatstimulusonsetandoffset(see\nA2.2),wehavealsoincorporatedadeeperdiscussionofthepotentialmechanismsconnectingtranslaminarand\ninter-arealθtoγ andspikingactivityintheResults,atlines396-417.\nR2.12: Howdotheauthorsinterpretthefindingsthatspikesarealignedtothetroughandpeaksofthetaatdifferent\ncorticallayers? Thisspeaksagainsttheoscillatoryvoltagemodulation(i.e. neuronsaremorelikelytofireat\nthepeakofanoscillation),andit’salsodifficult(forme)tointerpretifweconsiderthetaasatranslaminar\ntravelingwave: dependingonthewave’sspeed,spikeswon’tbesynchronizedbythetravelingwave.\nA2.12: Wethankthereviewforidentifyingambiguityinoursummaryoftherelationshipbetweenθandspikingactivity.\nWehavenowperformednewanalysescomparingfiringratesandpreferredθphasesbetween‘hit’and‘miss’\ntrials,aspresentedinFigs.6cand6eto6g. Briefly,wefoundthatpopulationsofneuronsinsuperficiallayers\nthataresensitivetotheθ signalshowincreasesfiringrateduring‘hit’trialsandspikingactivitythataligns\ntowardθ. Thevariationofpreferredphasesbetweenhitandmisstrialsindicatesanadditionallinkbetweenthe\nstrengthofθpropagationandneuronalspikingeventhoughthiseffectmaynotbefullyseparablefrompotential\ndipole-relatedpolarityreversals[52]. WehavenowupdatedourschematicinFig.1tobetterreflectthesenew\nresults.\n11\n3. Reviewer#3\nR3.1: HarrisandGonghaveperformedaninterestingandnovelcharacterisationofoscillatoryactivityinprimary\nandhighervisualareasofthemousecortexbyanalysingdatacollectedbytheAllenInstitute. Theauthors\nprovidenovelinsightsintointegrationofvisualinformationalongcorticaldepthsandacrossvisualareas. It\nisproposedthatthetaandgammaoscillationcarryfeedbackandfeedforwardrelatedinformationatdifferent\ntemporalepochs. Specifically,theyidentifyshortgammaburststhatoccurbasedonfeedforwardprocessing,\nwhichinturnmultiplexwithtravellingthetawavesinafeedbackdirectionaftertheonsetofthestimulusin\nachangedetectiontask. Overall,theanalysesareveryinterestingandbringnewperspectivestothefield.\nHowever,therearemanyconcernswiththeinterpretationoftheresultsintermsoffeedbackandfeedforward\nsignalsbothacrossdepthandacrossareas. Oneadditionalaspectthatisreallyinteresting,butdoesn’tseem\ntobethefocusofthestudyisthedifferencebetweentask relatedactivityandnon-taskrelatedflashes. While\nweareenthusiasticregardingthepotentialofthiswork,thereremainmanyconcernstoberesolved.\nA3.1: Wethankthereviewerfortheirconstructivefeedbackandpositiveassessmentofofourwork. Wehavemade\nsubstantialadditionstothemanuscriptinresponsetothereviewer’scomments,whichweaddresspoint-by-point\nbelow. We believe these additions support our interpretation of feedforward and feedback signaling across\ncorticaldepthsandareas,andelaborateonthedifferencesbetweentask-relatedandpassivestimulus-evoked\nactivitytobetterdescribethetop-downandbottom-upfunctionalrolesofspatiotemporalθ–γ patterns.\nR3.2: MajorConcerns: 1. Thetaisdefinedtobebetween3-10Hz. Thisisarelativelywiderangeoflowfrequencies\nwhichhappentoincludeboth3-5Hzrhythmsandtheta(6-9Hz). These3-5Hzrhythms,whichhavebeen\nlikenedtohumanalpha(seeSenzaietal. 2019andNestvogelandMcCormick2022),areadominantfeature\nofmousevisualcortexandshouldbeconsideredwhenperformingtheseanalyses. Wesuggesttheauthors\nreanalysethedataastheyhavedonethroughoutfigures1-3,butonlyconsideringanarrowerthetarhythm\n(6-10Hz,forinstance),whichalsomatchestheirfindingsfromthePhase-AmplitudeCouplingpresentedin\nfigures5andS5. Inaddition,itwouldbemostrelevant,alsobasedontheresultsfromFigureS5,forthe\nauthorstoperformasimilaranalysisontheselower frequency3-5Hzrhythms.\nA3.2: Assuggestedbythereviewer,wehaverepeatedouranalysisofθpropagationfortwonarrowerbands: 3–5Hz\nand6–10Hz. Wefoundqualitativelysimilarresultsofbothtranslaminarandinter-arealθpropagationforboth\nbands. Fortranslaminarandhierarchicalpropagationafterstimulusonset,ourfindingsforbothnarrowerbands\naresimilartothosefor3–10Hzband,withtheexceptionofincreasednon-causalartifactseithersideofthe\nstimulusforthe3–5Hzband. Forpropagationafterstimulusoffset,the3–5Hzbandshowsstrongersuperficial-\nto-deeptranslaminarpropagationandfeedforwardhierarchicalpropagation. Wehavenowincludedtheseresults\nintheSupplementarymaterial,asFigs.S13andS14. Promptedbythisanalysis,whichsuggestsoffsetθismore\nprominentinthe3–5Hzband,weaimedtoexplicitlycharacterizetheoffset-evokedθoscillations,inparticular\ntheirdominantfrequency. Wehavenowaddedanewfigure(Fig.2)totheSupplementalMaterialthatshows\nthefrequency,amplitude,anderrorfordampedsinusoidsfittotheoffset-evokedθ. Thisanalysisindicatesthat\nthedominantfrequencyofevokedθoscillationssitsat∼4HzforVISp,butincreasesto∼6.5HzforVISrl,\nhighlightingthatanarrow3–5Hzbanddoesnotcapturethedominantlow-frequencydynamicsforallvisual\nareas:\n12\n4 5 6\nAverage frequency (Hz)\n)%\n(htped\nlacitroC\na Frequency b Amplitude c Error\nVISp\nVISl\nVISal\nVISpm\nVISam\n100\n1 1.25 1.5 0.75 0.8 0.85 0.9\nAverage normalized amplitude Average RMSE\nFigure2:Transientoffset-evokedθoscillationsarestrongerandslowerinlower-ordervisualareas.\nProny’smethodwasusedtofitdampedcomplexexponentialstothe750msofLFPfollowing\nstimulusoffset.TheLFPwasbandpassfilteredbetween1–30Hzandz-scoredacrosstimeforeach\nchannel.Thedominantcomponentforeachchannelwasselectedastheonewiththesmallest\nroot-mean-squareerrorrelativetotheoriginalsignal.\naMedianfrequencyofoffset-evokedθoscillationsacrosscorticalareasanddepths,averagedacross\nsessions.Bandsshow95%CI.\nbMedianamplitudeofthedominantoscillatorycomponent.\ncMedianroot-mean-squareerrorofthedominantcomponentrelativetotheoriginalsignal.\nBasedonthisnewanalysisofoffset-evokedθoscillations,wehavealsoaddedadditionalmotivationforour\nchoiceofthe3–10HzbandtotheMethods:\nBasedonourspectralanalysis,wechosetooperationallydefinetheθbandasthefrequencyrange\nbetween3Hzand10Hz.Thischoiceofθfrequencywasinformedbythreemainfactors.First,we\naimedtocoverthemajorityoflow-frequencypeaksabovethe1/f componentacrossallstimuliand\nvisualareas,asshowninFig.2andFigs.S5andS6,tofullyresolvethetaskandregiondependencies\noftheslowLFPcomponent.Additionally,weshowinFig.2howthedominantfrequencyofθ\noscillationsevokedbystimulusoffsetcanvarybetween4Hzand6.5Hzacrossvisualregions,\nsuggestingabroaderrangeoffrequenciesisneededtosimultaneouslycapturethedominant\nlow-frequencydynamicsacrossvisualareas.Second,sinceourmethodsbasedoninferring\npropagationfromtheHilberttransformofthespatiotemporalsignal,thefrequencyrangemustbe\nnarrowenoughtoavoidsignificantnegative-frequencycomponents.Third,asweaimedtoanalyze\ntask-relateddynamicsevokedbyvisualstimuliwithashort(∼250ms)duration,wechosea\nfrequencybandwideenoughtoreducenon-causalfilterartifacts.Thischoiceofθalsocoversthe\nmajorityoflow-frequencybandsthathavebeenlinkedtovisualprocessinginmice,including\n3–7Hz[23],0–10Hz[53],3–5Hz[22],3–6Hz[45],2–8Hz[54],3–8Hz[55],4–8Hz[56],and\n3–10Hz[57].Werepeatedouranalysisfornarrowerfrequencybands(3–5Hzand6–10Hz),finding\nqualitativelysimilarresultsbutweakerbehavioralcorrelationsanddecodingaccuracy(seeFigs.S13\nandS14).\nLines760-776\nR3.3: 2. Thekeyresultofthemanuscriptisthatactivityisfeedbackdominatedrightafterstimulusonset, and\nfeedforwardrightafterstimulusoffset. Thisresultseemscounter-intuitiveandpotentiallyinconflictwith\npreviousstudiesthatcharacterizethis(Samedoetal,Nat. Comm.,2022andothers). Belowareafewthoughts\nregardingthis. a. Whycounter-intuitive: Whenavisualstimulusturnson,informationregardingitneeds\ntomoveforwardthroughthehierarchy,sodifferentpartsofthesystemgettheinformation. Therefore,the\nfinding that information is feedback after stimulus onset is counter-intuitive. b. Is the result above only\nrelatedtothetawaves,andinsteadspikingactivityisfeedforwardaspreviouslyobserved(e.g. Samedoetal,\n2022)? Itwouldbeusefultounderstandthisasitwouldincreasethepotentialimpactofthestudy. c. Ifit\nisallrelatedtothetaskanddifferentinthetaskunrelatedviewingcondition,itwouldbegoodtohighlight\nthesenotjustatthestart,butacrossthevariousresultsofthemanuscript.\nA3.3: Assuggestedbythereviewer,wehaveclarifiedourinterpretationofthevariousrolesofθandγbyaddinganew\nschematicfiguretoFig.4. Wenowdiscussthisschematic,andourinterpretationofonsetandoffsetθandγ,in\ntheDiscussionatlines615-636andlines637-659. Specifically,afterstimulusonset,weobservetop-down\nθ (showingdeep-to-superficialtranslaminarpropagation,andhigher-to-lowerhierarchicalpropagation)that\nwe interpret as top-down priming. The presence of top-down θ does not preclude feedforward information\nflow[58]—fastγ packetsandneuronalspikingcanstillcarryrapidbottom-upsensoryinputsatstimulusonset,\nassummarizedinFig.4e. Moreover,transient,spatiallylocalisedγ packetscanmultiplexbothfeedforward\nandfeedbackinfluencesacrosslayersandpathways. Finally,aswehavefoundthatθafterstimulusonsethasa\ntop-downrole,itisnaturallytask-dependent,asindicatedbytheweakerorderparametersafterstimulusonset\npresentedinFig.S3(passivenaturalimages). Wehavediscussedthisresultinthemaintexttohighlightthat\nonsetθistask-dependent,whereasoffsetθislargelytask-independent.\nR3.4: 3. Thegeneralassumptionacrosstheresultsisthatdeeptosuperficialdirectionisfeedbackandviceversa.\nHowever, L4 / L6 are where thalamic inputs are coming into visual areas (particularly VISp and VISl)\nandfeedforwardispresumablyseededfromhere. Howdotheresultsofthemanuscriptcomparewiththis\nanatomicalconnectivity? Perhapsamoregeneralandtruestatementwouldbetocallitdeeptosuperficial,\nratherthancallingitfeedbackintheresults. Ofcourse,thiscanbediscussedaspotentiallyfeedbackinthe\ndiscussionpartofthemanuscript.\nA3.4: Wethankthereviewforhighlightingapotentialsourceofconfusioninourmanuscript. Insectionθ waves\npropagatebidirectionallyacrosslaminaeandpredicttaskperformanceoftherevisedmanuscript,wenowrefer\toonset(offset)θwiththephrasing‘deep-to-superficial’(or‘superficial-to-deep’;basedonourtranslaminar\nwavenumber analysis) and, after introducing our task-related findings, ‘top-down’ (or ‘bottom-up’). In the\nfollowing section, Inter-areal θ propagation switches hierarchical directionality, we have made explicit at\nlines355-362thatanyfeedback/feedforwardorientationisinferredfrominter-arealhierarchicalphasedelays\nrather than from the translaminar results presented in the earlier section. After presenting our findings on\nhierarchicalθpropagation,wethenusethephrasingof‘feedback-oriented’(‘feedforwardoriented’)torefertoθ\ndirectedfromlower-to-higher(higher-to-lower)areasintheanatomicalandfunctionalhierarchies. Notably,\nthespatialdistributionof‘feedforward’and‘feedback’θhierarchicalpropagationbroadlymatchesthelaminar\ndistribution of the dominant feedforward (superficial; layer 2/3 and 4) and feedback (deep, layer 5 and 6)\ncortico-corticalpathways[33]. Finally,ouranalyseshavebeendesignedtostudycortico-corticalθpropagation;\nthey do not specifically address thalamocortical inputs to layer 4, or thalamocortical interactions that may\ninfluenceapparentθpropagation.\nR3.5: 4. Isthisaveryspecifictaskrelatedfeaturethatisbeingdiscovered? a. Becauseitisachangedetection\ntask,itwouldbeusefultolookatstimulusevents,repetitionevents,andthosethatthechangeoccurredand\nshowtheresultsseparately. b. Isthefeedbackresponseafterstimulusonsetactuallyrelatedtoinformation\nfromthepreviousstimulus–givenitisachangedetectiontask? SomeoftheAllenchangedetectiontasks\nincludedomissiontrials. Itwouldbereallyinterestingtocharacterizethosetrials. Presumably,thefeedback\nsignalswillpersistatthetimeofexpectedstimulus, andthefeedforwardwouldbeabsent. Thiscouldbe\naveryinterestingandsignificantresult. c. Relatedtotheabove,ithasbeensuggested(JordanandKeller,\nNeuron2020)thatfeedbacksignalsareactuallyarisinginL2/3. Againthisiscontrarytotheterminology\nlooselyusedinthemanuscript.\nA3.5: Wethankthereviewerforraisingtheseimportantpointsabouttaskspecificity. First,regarding,pointa),wehave\nnowanalyzedadditionalstimulusepochs,addingasupplementalfigureshowingthestrengthofθpropagation\nduringthepassivenaturalimagesstimulusset,atFig.S3. Inthiscondition,micewereshownasimilarsetof\nnaturalimagesasduringtheactivechangedetectiontask,butwerenotrewardedforanybehavioralresponse.\nHowever,the‘AllenNeuropixels–VisualBehavior’datasetstudiedinourworkdoesnotincludeomissionsfor\nchangetrials(onlynon-changetrials).\nSecond,forpointb),wehavenowexpandedourDiscussionatlines679-695,suggestingthatfutureworkcan\nincorporatevariedstimulusconditionssuchaschange-omissiontrialstodeterminewhetheronsetθisrelatedto\ntheanimals’backgroundlevelofengagement,lockedtotheexpectedtimingofthestimulusonset,orispurely\ninitiatedbystimulusonset.\nThird,forpointc),wehaverefinedourphrasingthroughoutthemanuscripttoclarifythatthebottom-upand\ntop-downpatternsweobservearethesignaturesthatdominatethelocalfieldpotentialattheθandγ scales;\ntheseLFP-dominantsignaturesdonotnecessarilycaptureallcomponentsoftheinformationprocessingthat\noccursinagivenlayer. Indeed,JordanandKeller[59]usedintracellularrecordingsduringavisuomotortask\ntofindthatcoexistingsub-populationsofneuronsinlayer2/3respondinoppositedirectionstothedifference\nbetween‘top-down’motorinputsandbottom-upvisualinputs.\nR3.6: 5. Istheoffsetresultparticulartothetimingofthestimulus? Isitadelayedresponsetostimulusonset? For\nexample,ifalongerstimulusistested(presumablytherearedifferentstimuluslengthsinthedataset),isthe\nfeedbackdominantperiodlockedtoafixedtimeafterstimulusonsetorlockedtostimulusoffset? Thiswould\nbeimportanttocheck.\nA3.6: Wethankthereviewerforraisingthispoint. The‘AllenNeuropixels–Visualbehavior’datasetstudiedinthis\nworkdoesnotincludedstimuliofvaryingdurations. Weagreethatthedecisivetestistovarystimulusdurations\nandincludeomissiontrials;wehavenowincorporatedthislimitationandopportunityforfutureworkintothe\nDiscussion,atlines679-695.\nR3.7: 6. Foranalysingthedirectionalitybetweenareas,theanalysesareperformedacrossthesamecorticaldepth\nbetweenareas. Suchananalysiswouldbeappropriateiftheconnectivitybetweenareasisorderedinthe\nsimilarway,whereL4iscommunicatingtoL4,L5toL5andsoon.However,theactualbiologicalconnections\nareofteninterlaminar. Istheanalysisthenvalid?\na. Whywouldtopdown/feedbackinformationbeprocessed\ninL6? Ifanything,L2-3wouldbemorelikelytointegratefeedbackinformationfromhigherareas?\nA3.7: Wehavenowclarifiedthemotivationforourapproachtoanalyzinginter-arealθinInter-arealθpropagation\nswitcheshierarchicaldirectionality,atlines363-379. Specifically,ourprimarymotivationforcomparing\nphase delays between equivalent depths across cortical areas was to compare the functional activation of\nequivalentlayersacrosstheanatomicalhierarchy,basedontheassumptionthatthesamelayersinadjacentareas\ninthehierarchyhavesimilarfunctionalroles(especiallyasreflectedbytheθsignal). Withthisapproach,wecan\naddressthehypothesisofactivitypropagationbetweenareas,evenifwedonotcapturethedetailedinput-output\narchitectureoflaminaeortheprecisepatternofdelaysbetweenlayers.\nHowever,wehavefoundthatfeedforward-orientedinter-arealθpropagationismostcoherentinmoresuperficial\nlayers,whereasfeedback-orientedθpropagationismorecoherentindeeplayers;whilethisfindingdoesnot\nprecludethepresenceofotherinterlaminarpathways,itisbroadlyconsistentwiththeanatomyofaprominent\nfeedforwardpathway(residinginL2/3and4)andaprominentfeedbackpathway(residinginL6)[33]. We\nalsofoundlayer-specificθ–γ couplingthatalignswithotherfeedforwardandfeedbackanatomy. First,atonset,\nstrongerθ–γ phase–amplitudecouplingappearsinL1oflower-orderareas(acanonicalfeedbacktarget)andL6\nofhigher-orderareas. Second,atoffset,layer2/3and4oflower-orderareasshowthestrongestθ–γ coupling,\nconsistentwiththepresenceofasuperficialfeedforwardpathway. Wehavenowemphasizedthisreconciliation\nandadoptterminologyintheDiscussion(‘feedback-oriented’)thatisconsistentwithknownpathwaybiases.\nprojections;futureworkwillanalyzecross-layer,cross-areaphaserelationsandthalamocorticalcontributions\nexplicitly.\nR3.8: 7. WouldbegoodtoidentifylayersbasedonCSD.Alltheinferencesarebasedondeeptosuperficialbeing\nL5toL2/3anditwouldbeimportanttoensureL4andL5arenotmislabelled. Whilethenewmethodsdo\nappeartobemorereliable,itwouldgiveconfidencetoreadersmorefamiliarwithCSDtohavethosemethods\ntoconfirmtheidentificationoflayers.\nA3.8: Assuggestedbythereviewer,wehaveaddedananalysisofcorticallayersinferredfromcurrent-sourcedensity\nprofiles. AspresentedinFig.3, weusedtheacanonicalprincipletoidentifythecenteroflayer4: thatthe\nstrongest,earliestcurrentsinkmarkstheinitialthalamicinputtolayer4. WefoundthatCSD-basedestimatesof\nthelayer4centerlargelyagreedwithestimatesbasedonmappingchannelCCFv3coordinatestotheanatomical\natlas,butthatlayeridentificationwaspoorinhigher-orderregions,presumablyduetoweakerstimulus-related\nthalamicinputandmorecomplexarrangementsofsourcesandsinks. Basedonthesefindings,inadditionto\nrecentworkshowingthatCSD-basedlayeridentificationbesensitivetosubjectivejudgementswithoutadditional\ninformationsuchasprofiletemplates[60]ordynamicalfeatures[61],wechosetousethemappingtoCCFv3\ncoordinatesasourmainsourceofinformationregardingcorticallayerboundaries.\nCSD L4 Anatomical L4\n)mμ(\nhtpeD\nVISp\n0\nρ = 0.73 L1\nL4\n)mμ(\nhtpeD\nVISl\n0\nρ = 0.14 L1\nL4\n)mμ(\nhtpeD\nVISrl\n0\nρ = 0.5 L1\nL5\n)mμ(\nhtpeD\nVISal\n0\nρ = 0.45 L1\nL4\n)mμ(\nhtpeD\nVISpm\n0\nρ = 0.23 L1\nL5\n)mμ(\nhtpeD\na b\nc d\ne f VISam\n0\nρ = 0.55 L1\nL4\n800\nindependentlyusingeitherthelocationofthefirst,largestsinkinCSD,ortheCCFv3layer\nannotationsofeachchannel.\naThedistributionofestimatedL4centersinVISpacrosssessions,usingtheCSD(left)ortheatlas\nannotations(right).TheSpearmancorrelationbetweenCSDestimatesandatlasestimatesis\nannotated(topleft).\nb–fSameasa,butforothervisualcorticalareas.Theatlaslayershaveagreaterdegreeof\nconsistencythantheCSD-basedlayers,especiallyinhigher-orderareaswheretherelativestrengthof\nthalamicinputtolayer4isweaker.\nR3.9: 8. Intermsofexplainingtheresultsafewdescriptiveillustrationswouldreallyhelp. Thefirstpartoffigure\n1wasveryusefulinexplainingtheconceptsforthereadertofollowthefirstfewfigures. Similartothisit\nwouldbebeneficialtohave:\na. afiguredescribingtheinter-areaanalysis.\nb. asummaryfiguretoexplainalltheresultsandtheirtemporalorder.\nA3.9: Followingthereviewer’ssuggestion,wehaveaddedschematicfigures(Figs.4eand4f)thatsummarizeour\ninter-arealanalysis,therelationshipbetweentranslaminarandcross-areaθpropagation,theinteractionsbetween\nθandγ components,andthetemporalsequenceofeffects(onsetperiodtooffsetperiod).\nR3.10: MinorRevisions. Line87: Wefindthat activityisorganizedaslarge-scaletravelingwavespropagating\nbetweenlaminaeandacrosscorticalareas.\nForthereaderthatisnotfamiliarwithsuchcharacterisations,it\nwouldbeusefultoclarifytheevidenceforsuchstatements.\nA3.10: Wehavenowclarifiedthatbothtranslaminarandinter-arealθpropagationaremeasuredusingtheinstantaneous\nphasegradientatlines118-141oftheResults.\nR3.11: Figure2,has5differentcolorschemes. PanelA,C-Ghavecolorsinareas;Bcolorondepth,Eonexponent\norintercept,andHcoloronfreqband+anadditionalbandsalongtherightedge. Itisalltoomuchand\nwouldbebettertoeithersplitthefiguresorreducethedifferentcolorschema\nA3.11: WehavenowsimplifiedthecolorschemeforFig.2,reducingthenumberofpalettestotwo: oneforthesix\nvisualareas,andanotherforthefivecorticallayers. Thespectralmeasureshavebeenre-coloredforconsistency\nwiththevisual-areapalette,maintainingcontrastwiththeadjacentverticalbarsusedtoindicatelayerboundaries.\nR3.12: Examplepanelsofthetheta-gammainteractionsarebasedonresultsfromahighervisualarea(VISL).I\nthinkthemaintextshouldberevisitedtoreinforcethis,orthemainfiguresshouldbechangedtoincludea\nmoregeneralviewoftheprimaryvisualcortexaswell.\nA3.12: WehavenowclarifiedhowthepatternsofspatiotemporalPACchangeacrossvisualareasinthemaintext,at\nlines470-493. WechoosetopresentVISlinthemainfigureasitisamorerepresentativesampleofareasin\nthevisualhierarchythanVISp,givenVISpdoesnotreceivefeedforwardcorticalinput[62]. Nevertheless,we\nfoundanalmostidenticalpatternforspatiotemporalθ–γ PACinVISp,asshowninFig.S10.\nR3.13: Lines144-147,thehigherresidualpowerofthetarhythmsinthedeeperlayersofVisualCortexmayalso\nreflectpassivevolumeconductionratherthananactiveprocessoftop-downinformation. Canahippocampal\ntravellingwavebediscounted?\nhippocampus. Wenowprovideamorethoroughtreatmentofthevolumeconductionhypothesisandtheevidence\nagainstitintheDiscussion,atlines679-695.Whilewecannotexcludeacontributionfromvolumeconduction,\nthree observations make a hippocampal volume-conduction account less likely. First, deep-to-superficial θ\npropagationismostprominentinlower-ordercorticalareasandprogressivelyweakerinhigher-orderregions\n(see Fig. 3), but residual θ power is strongest in higher-order areas (see Fig. 2). If the signal were broadly\nvolumeconductedfromthehippocampus,itspropagationcoherenceandamplitudeshouldberelativelyuniform\nacrossthepatchofcortexoccupiedbythevisualareas. Second,inter-arealθshowsordered,non-zerophaselags\nwithclearhierarchicaldirectionality(seeFig.4),whichcontrastswiththenear-zerophasedifferencespredicted\nbyvolumeconduction. Third,bothtranslaminarandinter-arealphase-gradientdirectionsreverseatstimulus\noffsetandshowtasksensitivity—behaviornotexpectedfromastationary,externallygeneratedfield. Taken\ntogether, thesepointsargueagainstapurelyvolume-conductedhippocampalorigin, althoughwehavenow\nstatedexplicitlyintheDiscussionthatwecannotrule-outapartialcontribution.\nR3.14: Theresultsofthethetawavepropagationareinteresting(Fig3),butthetitleismisleading. Thisseemsto\nbetrueonlyinL,notinotherareas. AdditionallyIwonderifmoresuperficialwavesarebeingdominated\nbydelta/alpharhythms,whilethedeeperlayersaremoreinfluencedbythetarhythms? Arethesepatterns\npresentalsointhespontaneousstimulus?\nA3.14: WeobservequalitativelysimilarθpropagationinregionsVISp,VISl,VISrl,andVISal,asdemonstrated\ninFig.S7;topreventpotentialconfusion,wehavenowupdatedthetitleofFig.3toclarifythatthebidirectional\npropagationisstrongestinlower-orderareas.\nWehavealsoaddedasupplementaryfigure,Fig.S3,tohighlighthowθpropagationispresentacrossallstimuli,\nbutisnotablystrongerfornaturalimagesandflashes.\nR3.15: Lines214-216: CanthisbearesultofVISphavinghigheramplituderesponses?\nA3.15: Ourindexofpropagation—theorderparameter—isnotdirectlysensitivetotheLFPamplitude, onlytothe\ncoherenceofθ-phasespatialgradients. Thispropertyhasadvantagesfortrackingsignalpropagationacross\nlayers and regions that may have different intrinsic response amplitudes. Nevertheless, it is expected that\ningeneral,higher-amplituderesponsesarecorrelatedwithhigherpropagationcoherence. Assuch,wehave\nincludedcomparisonstothemagnitudeoftheLFPresponseforourclassificationresultsinFig.3f,findingthat\nthemeanLFPamplituderesponseislessinformative(especiallypriortothestimulusoffset)thantheorder\nparameterandcompartment-specificLFPresponses.\nR3.16: Lines246-247,isthiswrong,itseemslikeitsstrongestinthedeeptosuperficialdirection?\nA3.16: Wehavenowupdatedourphrasinginthissection,atlines332-352,toemphasizethatthesuperficial-to-deep\nθpropagationfollowingstimulusoffsethasbottom-upcharacteristics(incontrasttoθpropagationfollowing\nstimulusonset).\nR3.17: ThesignificanceoftheAparameterpresentedinfigure4Bisunclear. Themethodsstate‘AhighAindicates\nlarge,consistentlyorientedphasedelaysalongthevisualhierarchy,whereasasmallAindicatesthatphase\ndelaysareeithersmall(meaning ismoresynchronousovertheregions)orinconsistentlyoriented’(Lines\n644-646). Sobasedonthefigure,darkredmeanshighlysynchronoustheta,whilebluemeansconsistently\ndelayed?\nA3.17: We thank the reviewer for flagging this imprecise description of A . We have adjusted the relevant text at\nθ\n18\nlines881-899toreflectthatthemagnitudeofA determinesthecoherenceofphasedelaysacrossregions,\nθ\nandthesignreflectsthedirectionofapparentθpropagationalongthevisualhierarchy.\nR3.18: ThePACanalysesofalltheareasandlayersarereallyinteresting. Butthefiguresinthemaintextarefor\nthespontaneousstimulus,andnotthetask-onethatwaspresentedinpreviousfigures. InfactthetaskPAC\nshowslittlecouplingbetweenthetaandgammainVISL.Howcanthisdifferencebeexplained?\nA3.18: Theaveragedcomodulogramsshowninthesupplementarymaterialarecomputedacrossanentirestimulus\nepoch,which(forthetaskcondition)mixestogetherthevariousθ andγ dynamicsweobservedatstimulus\nonsetandoffset. Inparticular,thecomodulogramsforthetaskandflashconditionsinlower-orderareasappear\ntobedominatedbya4Hzθ component,whichislikelyrelatedtothe250msstimulusperiod. Assuch,we\nhaveextendedthisclassical‘stationary’θ–γ analysistothe‘spatiotemporal’θ–γ couplingcalculatedovertrials,\nwhichwepresentinthemaintextandbelievebettercapturesθ–γ couplinginthepresenceofstimulus-related\nnonstationarity.\nR3.19: Inthediscussion(lines470-472),theauthorsmentionthatthetawaslikenedtoalpharhythmsinprimates.\nThisisnottrue. Inbothcitedpapers,itwasthe3-5Hzrhythminmicethatwassuggestedtobecomparableto\nalpha,nottheta. Pleasereviseyourmaintextaccordingly.\nA3.19: Wethankthereviewerforhighlightingthispotentialsourceofconfusion. Wehavenowadjustedourphrasingin\ntheDiscussion(atlines637-659)toclarifythatthe3–5Hzrhythmismostsimilartoprimatealpharhythms,\nratherthanthebroaderθbandconsideredinourwork.\n19\n4. Reviewer#4\nR4.1: Ico-reviewedthismanuscriptwithoneofthereviewerswhoprovidedthelistedreports. Thisispartofthe\nNatureCommunicationsinitiativetofacilitatetraininginpeerreviewandtoprovideappropriaterecognition\nforEarlyCareerResearcherswhoco-reviewmanuscripts.\nA4.1: WethankthereviewerfortheirtimeandeffortinevaluatingourmanuscriptaspartoftheNatureCommunications\nECRinitiative. Thispointrequiresnochangestothemanuscript,butwearepleasedtorecordourgratitudehere.\n5. Reviewer#5\nR5.1: Inthismanuscript,theauthorshaveanalyzedthepubliclyavailabledata(AllenInstitute)toexaminethe\nspatiotemporaldynamicsofLFPacrosscorticallayersandvisualcorticalareas,whilehead-fixedmicewere\npresentedwithflash,grayscreen,andnaturalimages(change-detectiontasks). Theyhaveinvestigatedthe\n“travelingwave”propertyoftheta(3-10Hz)oscillationsacrosscorticallayersindifferentcorticalareas,and\nsomepropertiesarecorrelatedwiththehierarchyscoresofvisualcorticalareas. Theyhavealsodetected\ngamma(30-100Hz)oscillationpacketsandhowthesearecoordinatedwith3-10Hzoscillationsindifferent\ncorticallayersindifferentvisualcorticalareas. Overall,theseanalysesaregreatadditionstothecommunity’s\nknowledge. Inthemeantime,Ihaveafewconcernsrelatedtotheinterpretationofthedata. (1)Theauthors\ninterpretthethetawaveassignalingthefeedbackinputs,citingstrongerthetapowerindeeperlayersthat\narerelevanttothefeedbackconnectivity. However,thisanalysisismostlybasedonthenon-task-engaging\nexperiments(flashscreenandgreyscreenexperiments,Fig2andFigS1),notsomuchonthetask-engaging\nexperiments(naturalimage,FigS2). Thus,“feedback”interpretationduringnon-task-engagingexperiments\ncouldnotbeastraightforwardinterpretation. Also,althoughthethetapowerisstrongerindeeplayers,time\nscalehierarchycorrelationsarestrongerinsuperficiallayers(Fig2H,FigS1H,FigS2H).Howwouldthis\nalignwiththeauthors’interpretationofthetaas“feedback”?\nA5.1: Weappreciatethereviewer’spositiveassessmentandconstructivecomments. WehavenowclarifiedinVisual\ncorticalLFPscompriseθ,γ,and1/f,asquotedbelow,thatheightenedθpowerindeeplayers,togetherwitha\npositivegradientalongthecorticalhierarchy,ismostconsistentwithafeedbackorigin.Wehavealsoemphasized\nintheResults,atlines213-242,thatthisspectralanalysisisnotsufficienttofullyaddressthefunctionalrole\nofθorγ,especiallyinthecontextoftaskstimuli. Oursubsequentanalysisoftheθcomponentweretargeted\nat addressing the top-down nature of onset θ. First, we found that near stimulus onset, the coherence of θ\npropagationissignificantlystrongerduring‘hit’trialsthan‘miss’trials(seeFigs.3aand3b),weakerduring\npassivenon-changetrials(seeFig.S3),andabsentduringpassiveflashtrials(seeFig.3c);thispointstothe\ntop-downfunctionalroleofθ. Second,θhasaconsistentdirectionalityfromhigher-ordertolower-orderregions\nfollowingstimulusonset,especiallyindeeperlayers,suggestingthetop-downfunctionofθisimplementedby\nfeedbackpropagationfromhighercorticalareas(seeFig.4). Third,wefoundthatθ–γ couplingfollowingthe\nstimulusonsetalignswithfeedback-consistenttargets(beingstrongestinL1oflower-orderareasandL6of\nhigher-orderareas;seeFig.5). Fourth,weobservedthatsupragranularfiringforθ-lockedneuronswasenhanced\nandlockedmorecloselytoθtroughson‘hit’trialscomparedto‘miss’trials(seeFig.6). Finally,wenotethat\nakeyresultofourworkisthatθhasadualfunctionalrole;weusedasimilarapproachtoshowthatθhasa\nbottom-up/feedforwardroleafterstimulusoffset.\nWenextmappedthedistributionofθandγacrosslayersandcorticalregions.Wecalculatedtheresidual\npower,givenbyEq.(5),astheratioofthepowerabovethe1/f componentforaspecificfrequencyband\ntothetotalpowerofthesignal.Duringtheflash_250msstimulus,theresidualθpowerprogressively\nincreasesfromsuperficialtodeeplayers(medianKendall’sτ =0.37,allp<10−7),asshowninFig.2f.\nTheprominenceofθpowerindeeplayersmayreflectafeedbackorigingiventheanatomical\ndistributionoffeedbackconnectivityandtheassociationofdeep-layersomawithfeedback\nprojections[63,28].Totestthisfeedbackhypothesis,wecomparedthedistributionofresidualθpower\ntotheanatomicalhierarchyscoreofeachregion.Thesehierarchyscoreswerecalculatedforvisual\n20\ncorticalareasbySiegleetal.[18],usingthetracer-basedanatomicalconnectomeofHarrisetal.[25].\nAsshowninFig.2h,wefoundthatresidualθformsastrongpositivegradientalongtheanatomical\nhierarchyacrossalldepths,peakinginlayer2/3withKendall’sτ ≈0.5(p<0.01forallcorrelations\nafterBenjamini-Hochbergcorrectionacrossdepths).Wealsofoundsimilarpatternsinθorganization\nacrossdifferentstimuli,includingspontaneous(seeFig.S5)andNatural_Images(Fig.S6).\nLines197-212\nR5.2: (2) The authors interpret the gamma packet as signaling the feedforward inputs, citing stronger gamma\npowerinthesuperficiallayersisrelevanttothefeedforwardconnectivity. However,theseanalysesareagain\nmostlybasedonthenon-task-engagingexperiments. Indeed, inthenaturalimagetaskexperiments, the\ngamma-hierachyscorecorrelationsarehighindeeplayers. Thus,furtherdiscussiononhowgammacould\ncorrespondtothefeedforwardprocessingwouldbepreferable.\nA5.2: Wehavenowincorporatedourfindingsrelatedtohierarchicalgradientsofresidualγ intothemaintextinVisual\ncorticalLFPscompriseθ,γ,and1/f,atlines213-242.\nWehavealsoaddedaschematicsummarizingourinterpretationofθ andγ toFig.4, whichcoverstherole\noffeedforwardγ insuperficiallayers. WediscussthissummaryschematicinthemaintextatInter-arealθ\npropagationswitcheshierarchicaldirectionality:\nTogether,thesecomplementaryfindingsindicatethatθ-scaleactivityplaysadualroleinvisual\nprocessing,flexiblyswitchingbetweenfeedback-oriented(top-down)andfeedforward-oriented\n(bottom-up)modesaccordingtothevisualstimulusandtaskdemands.Figures4eand4fsummarizehow\ntranslaminarandinter-arealθmayswitchbetweenthesetwoprocessingmodesduringtasktrials.\nFirst,asdepictedinFig.4e,wesuggestthatstimulusonsettriggersarapidtransmissionof\nfine-grainedinformationaboutthevisualscene[17]alongafeedforwardpathway.Thisfeedforward\npassappearstodominateγLFPactivityinlayers2/3and4(seeFig.2),withtheearliestsignaturesof\nvisualinputappearinginVISplessthan50msafterstimulusonset[18].Ashorttimelater,\nbetween\n50–150ms,afeedback-orientedθsignalarrivesinsuperficiallayersofthelower-ordervisualareas,\nbeingtransmittedfromhigher-orderregionsalongadeeperpathway[19,20].Thisfeedbackθ\ncarriestop-downinformation,likelytiedtotaskengagement(seeFig.S3),thattravelsbetweenlayers\nineachregiontomeetandmodulatefeedforwardγinthesuperficiallayers.Atstimulusoffset,\nhowever,θplaysadifferentrole,asshowninFig.4f.Offset-evokedθoscillationsmaysynchronizethe\nfeedforwardγpacketsthattravelfromlower-orderregionstohigher-orderregionsinasuperficial\npathway.Theseθoscillations,whichhaveacharacteristicfrequencyof∼4HzinVISp(asshown\ninFig.2),resembletheevokedactivitypreviouslyobservedinlowervisualcortical\nareas[21,22,23,24].Here,wehaveshownhowthisoffset-evokedθoscillationpropagatesacrossthe\nvisualhierarchyandcorticallayers,whichmayservetoperiodicallyshieldneuralactivityinhigher\nregionsfromthefeedforwardsignal,reducingtheinterferenceofbottom-upinformationwith\ncognitiveprocessing.Thesefindingssuggestthatthespatiotemporalstructureofθaffordsitaflexible,\ncontext-dependentfunctionalrolethatadaptstothecomputationaldemandsofvisualtasks.\nLines396-417\nR5.3: Ialsohaveafewminortechnicalcomments/questions. (3)Theauthorshaveused“wavenumber”toevaluate\nthe property of theta as traveling waves, claiming that it can capture single-trial measurement of LFP\npropagation. However, major analyses are done using the average of it, not so much different from the\nconventionalCSDanalysis.\nA5.3: Ourwave-basedanalysishasanumberofcrucialfeaturesthatdistinguishitfromconventionalCSDanalysis,\n21\nbutmostimportantly, thewavenumberprovidesasingle-trialmeasurementofθ propagationthatcannotbe\ncapturedbytheaverageCSD.Bycomputingthewavenumberateachpointintimeandspace,foreachtrial,\nthentakingthemedianofthisstatisticacrosstrials,wearemeasuringtheaveragepropagationoftheθsignal.\nBycontrast,conventionalsingle-trialCSDhasapoorsignal-to-noiseratioduetothenoise-boostingeffectofthe\nsecondderivative. ExtractingmeaningfulinformationfromtheCSDrequiresaveragingacrosstrials,andany\ntravelingpatternsinferredfromtheaverageCSDwouldnotnecessarilyreflectsingle-trialpropagation. Thusthe\ntrial-averagedCSDcannotdisambiguatetruepropagationsthatoccurinasingletrialfromspuriouspropagations\nintheaveragesignal—forexample,localizedactivitythatoccursatdifferenttimepointsandlayersfordifferent\ntrialsmayappearasapropagatingcomponentintheaveragesignal.\nTohighlightthisdistinctionbetweentheCSDandthewavenumber,weshowtheresultsofCSDanalysisondata\nduringthetaskstimulusinFig.S8. Comparedtothewavenumberdescription,thisconventionalCSDanalysis\ncaptureslayer-specificpatternsofactivity,butisnotsensitivetothepropagationofθacrosslayers.\nR5.4: (4) Related to (3), are they really the translaminar traveling wave? Perhaps just interpreting the LFP\nfluctuationsmadebycurrentsinksandsourcesas“traveling”. Youcanalwaysassign“thetaphase”after\nbandpassfilteringoftheoriginalLFPsignalsinthethetaband,andthuscancomputethewavenumberout\nofit.\nA5.4: Weagreethattheprecisemechanismsunderlyingtranslaminarθpropagationremaintobeestablished. However,\nwehaveclarifiedthatwedidnotinferpropagationmerelybyband-passfilteringandassigningphase. Beyond\nfiltering,wequantifiedsystematicphasegradientsalongthelaminaraxisandsummarizethembyawavenumber\nandaphase-coherence/orderparameter. Thesequantitiesareconsistentwithdirectedpropagationandarenot\nexpectedfromastaticconfigurationofcurrentsinksandsources. Phasecoherenceismarkedlyreducedfor\nsurrogatecontrolsshowninFig.3andFig.S3. Asnotedabove,anddemonstratedinFig.S8,conventional\nCSDanalysisthataverageovertrialscancaptureamplitudestructure,butnotslowtravelingpatterns. Thetask-\ndependentmodulationofwavecoherenceanddirectionfurtherstrengthensthecasethatthesearebiologically\nmeaningfultravelingwavesratherthananepiphenomenon. Intherevisedmanuscript,intheDiscussion,we\nhaveemphasizedtheopportunityforfutureexperimentstodeterminewhichcircuit-levelmechanismsunderlie\ntranslaminarθpropagation:\nFutureinvestigationwillhelptoclarifythephysicaloriginsandcausalmechanismsunderlying\nthesenestedspatiotemporalpatterns.Inparticular,furtherworkisneededtoquantifythedegreeto\nwhichlocalfieldpotentialscanmodulateneuronalexcitability[52,64,65,66,67].Additionally,the\ncellularbasisoftravelingfieldpatternsremainspoorlyunderstood[52]:smoothlymigratingspiking\nsequencesoractionpotentials[68,69],distributedpopulationdynamics[70,71],ordipole\neffects[72,52]mayhaveredundantfieldsignatures,anddisentanglingthesemechanismsiscriticalto\nunderstandingtheinfluenceofcircuit-leveldynamicsonlarge-scaleLFPdynamics.Volumeconductionof\nhippocampalorsubcorticalLFPmayalsocontributetotheobservedpropagationofθ,especiallyin\nthedeep–superficialdirection[26,73].However,variationsinθpropagationcoherenceacrossvisual\nareas,switchesinpropagationdirection,andnon-zerophasedelaysbetweenvisualareasthatare\nhighlycorrelatedtotheanatomicalhierarchycollectivelyargueagainstthisvolume-conduction\nhypothesis.Finally,futureexperimentalworkshouldaimtoclarifywhetheronsettop-downθresults\nfromrapidbottom-uptriggersorprecedingattentional,engagement-related,orpredictivestates;this\nmayinvolvemorecomplexbehavioralparadigmsthatincludechange-stimulusomissionsorstimulus\npresentationsthatvaryinduration.\nLines679-695\n[1] B.HarrisandP.Gong,Plotdatafor\"nestedspatiotemporaltheta–gammawavesorganizehierarchical\nvisualprocessing\"(2025).\n[2] B.Harris,brendanjohnharris/spatiotemporalmotifs.jl: v0.2.0(2025).\n[3] B.HarrisandP.Gong,Processeddatafornestedspatiotemporaltheta–gammawavesorganizehierarchical\nvisualprocessing(2025).\n[4] AIBS, Allen Brain Observatory – Neuropixels Visual Behavior [Data set], (2022), http://portal.brain-\nmap.org/circuits-behavior/visual-behavior-neuropixels.\n[5] J.Bezanson,A.Edelman,S.Karpinski,andV.B.Shah,Julia: Afreshapproachtonumericalcomputing,\nSIAMreview59,65–98(2017).\n[6] B.Harris,brendanjohnharris/allenneuropixels.jl: v0.2.0(2024).\n[7] B.Harris,brendanjohnharris/timeseriestools.jl: v0.5.0(2024).\n[8] B.HarrisandS.C.Surace,brendanjohnharris/normalization.jl: v0.7.1(2024).\n[9] S.DanischandJ.Krumbiegel,Makie.jl: Flexiblehigh-performancedatavisualizationforjulia,Journalof\nOpenSourceSoftware6,3349(2021).\n[10] J.Gehring,D.Widmann,D.Kleinschmidt,R.Finnegan,B.Kamiski,C.Bowers,H.Perez,M.Molignano,\nM.Bouchet-Valat,P.K.Mogensen,T.Kelman,andN.Ignatiadis,juliangehring/bootstrap.jl: Bootstrap\nv2.4.0(2023).\n[11] G.Datseris,J.Isensee,S.Pech,andT.Gál,Drwatson: theperfectsidekickforyourscientificinquiries,\nJournalofOpenSourceSoftware5,2673(2020).\n[12] M. Vinck, M. van Wingerden, T. Womelsdorf, P. Fries, and C. M. A. Pennartz, The pairwise phase\nconsistency: Abias-freemeasureofrhythmicneuronalsynchronization,NeuroImage51,112–122(2010).\n[13] J. Aru, J. Aru, V. Priesemann, M. Wibral, L. Lana, G. Pipa, W. Singer, and R. Vicente, Untangling\ncross-frequencycouplinginneuroscience,CurrentOpinioninNeurobiology31,51–61(2015).\n[14] S.R.ColeandB.Voytek,Brainoscillationsandtheimportanceofwaveformshape,TrendsinCognitive\nSciences21,137–149(2017).\n[15] R. T. Canolty and R. T. Knight, The functional role of cross-frequency coupling, Trends in Cognitive\nSciences14,506–515(2010).\n[16] A. Hyafil, A.-L. Giraud, L. Fontolan, and B. Gutkin, Neural cross-frequency coupling: Connecting\narchitectures,mechanisms,andfunctions,TrendsinNeurosciences38,725–740(2015).\n[17] L. S. Meyerolbersleben, A. Sirota, and L. Busse, Anatomically resolved oscillatory bursts orchestrate\nvisualthalamocorticalactivityduringnaturalisticstimulusviewing,bioRxiv10.1101/2024.08.21.608936\n(2024).\n[18] J.H.Siegle,X.Jia,S.Durand,S.Gale,C.Bennett,N.Graddis,G.Heller,T.K.Ramirez,H.Choi,J.A.\nLuviano,P.A.Groblewski,R.Ahmed,A.Arkhipov,A.Bernard,Y.N.Billeh,D.Brown,M.A.Buice,\nN.Cain,S.Caldejon,L.Casal,A.Cho,M.Chvilicek,T.C.Cox,K.Dai,D.J.Denman,S.E.J.deVries,\nR.Dietzman,L.Esposito,C.Farrell,D.Feng,J.Galbraith,M.Garrett,E.C.Gelfand,N.Hancock,J.A.\nHarris, R.Howard, B.Hu, R.Hytnen, R.Iyer, E.Jessett, K.Johnson, I.Kato, J.Kiggins, S.Lambert,\nJ.Lecoq,P.Ledochowitsch,J.H.Lee,A.Leon,Y.Li,E.Liang,F.Long,K.Mace,J.Melchior,D.Millman,\nM. Oliver, M. Pachitariu, J. Perkins, M. Reding, D. Reid, M. Robertson, K. Ronellenfitch, S. Seid,\nC.Slaughterbeck,M.Stoecklin,D.Sullivan,B.Sutton,J.Swapp,C.Thompson,K.Turner,W.Wakeman,\nJ. D. Whitesell, D. Williams, A. Williford, R. Young, H. Zeng, S. Naylor, J. W. Phillips, R. C. Reid,\nS.Mihalas,S.R.Olsen,andC.Koch,Surveyofspikinginthemousevisualsystemrevealsfunctional\nhierarchy,Nature592,86–92(2021).\n[19] C.Y.L.Huh,J.P.Peach,C.Bennett,R.M.Vega,andS.Hestrin,Feature-specificorganizationoffeedback\npathwaysinmousevisualcortex,CurrentBiology28,114–120.e5(2018).\n[20] S.Weiler,M.Teichert,T.W.Margrie,andA.Cardona,Layer6corticocorticalneuronsareamajorroute\nforintra-andinterhemisphericfeedback,eLife13,RP100478(2025).\n[21] C. L. Zold and M. G. H. Shuler,Theta oscillations in visual cortexemergewith experience to convey\nexpectedrewardtimeandexperiencedrewardrate,JournalofNeuroscience35,9603–9614(2015).\n[22] M.C.Einstein,P.-O.Polack,D.T.Tran,andP.Golshani,Visuallyevoked3–5hzmembranepotential\noscillations reduce the responsiveness of visual cortex neurons in awake behaving mice, Journal of\nNeuroscience37,5084–5098(2017).\n[23] A.Speed,J.DelRosario,C.P.Burgess,andB.Haider,Corticalstatefluctuationsacrosslayersofv1during\nvisualspatialperception,CellReports26,2868–2874.e3(2019).\n[24] E.A.B.Horrocks,F.R.Rodrigues,andA.B.Saleem,Flexibleneuralpopulationdynamicsgovernthe\nspeedandstabilityofsensoryencodinginmousevisualcortex,NatureCommunications15,6415(2024).\n[25] J. A. Harris, S. Mihalas, K. E. Hirokawa, J. D. Whitesell, H. Choi, A. Bernard, P. Bohn, S. Caldejon,\nL. Casal, A. Cho, A. Feiner, D. Feng, N. Gaudreault, C. R. Gerfen, N. Graddis, P. A. Groblewski,\nA.M.Henry,A.Ho,R.Howard,J.E.Knox,L.Kuan,X.Kuang,J.Lecoq,P.Lesnar,Y.Li,J.Luviano,\nS.McConoughey,M.T.Mortrud,M.Naeemi,L.Ng,S.W.Oh,B.Ouellette,E.Shen,S.A.Sorensen,\nW. Wakeman, Q. Wang, Y. Wang, A. Williford, J. W. Phillips, A. R. Jones, C. Koch, and H. Zeng,\nHierarchicalorganizationofcorticalandthalamicconnectivity,Nature575,195–202(2019).\n[26] Y. Senzai, A. Fernandez-Ruiz, and G. Buzsáki, Layer-specific physiological features and interlaminar\ninteractionsintheprimaryvisualcortexofthemouse,Neuron101,500–513.e5(2019).\n[27] D.Mendoza-Halliday,A.J.Major,N.Lee,M.J.Lichtenfeld,B.Carlson,B.Mitchell,P.D.Meng,Y.Xiong,\nJ.A.Westerberg, X.Jia, K.D.Johnston, J.Selvanayagam, S.Everling, A.Maier, R.Desimone, E.K.\nMiller, and A. M. Bastos, A ubiquitous spectrolaminar motif of local field potential power across the\nprimatecortex,NatureNeuroscience,547–560(2024).\n[28] T.vanKerkoerle,M.W.Self,andP.R.Roelfsema,Layer-specificityintheeffectsofattentionandworking\nmemoryonactivityinprimaryvisualcortex,NatureCommunications8,13804(2017).\n[29] M.P.Morton,S.Denagamage,I.J.Blume,J.H.Reynolds,M.P.Jadi,andA.S.Nandy,Brainstateand\ncorticallayer-specificmechanismsunderlyingperceptionatthreshold,eLife12,RP91722(2024).\n[30] M.N.OudeLohuis,A.C.Canton,C.M.A.Pennartz,andU.Olcese,Higherordervisualareasenhance\nstimulusresponsivenessinmouseprimaryvisualcortex,CerebralCortex32,3269–3288(2021).\n[31] B.Schuman,S.Dellal,A.Prönneke,R.Machold,andB.Rudy,Neocorticallayer1: Anelegantsolutionto\ntop-downandbottom-upintegration,AnnualReviewofNeuroscience44,221–252(2021).\n[32] S.Huang,S.J.Wu,G.Sansone,L.A.Ibrahim,andG.Fishell,Layer1neocortex: Gatingandintegrating\nmultidimensionalsignals,Neuron112,184–200(2024).\n[33] E.BarzegaranandG.Plomp,Fourconcurrentfeedforwardandfeedbacknetworkswithdifferentrolesin\nthevisualcorticalhierarchy,PLOSBiology20,e3001534(2022).\n24\n[34] M. Gao, S. Lim, and A. A. Chubykin, Visual Familiarity Induced 5-Hz Oscillations and Improved\nOrientationandDirectionSelectivitiesinV1,JournalofNeuroscience41,2656–2667(2021).\n[35] Y.Tang,C.Gervais,R.Moffitt,S.Nareddula,M.Zimmermann,Y.Y.Nadew,C.J.Quinn,V.Saldarriaga,\nP.Edens,andA.A.Chubykin,Visualexperienceinduces4–8hzsynchronybetweenv1andhigher-order\nvisualareas,CellReports42,10.1016/j.celrep.2023.113482(2023).\n[36] J.M.Levy,C.L.Zold,V.M.K.Namboodiri,andM.G.HussainShuler,Thetimingofreward-seeking\naction tracks visually cued theta oscillations in primary visual cortex, Journal of Neuroscience 37,\n10408–10420(2017).\n[37] S. Ni, B. Harris, and P. Gong, Distributed and dynamical communication: a mechanism for flexible\ncortico-corticalinteractionsanditsfunctionalrolesinvisualattention,CommunicationsBiology7,550\n(2024).\n[38] D.B.NestvogelandD.A.McCormick,Visualthalamocorticalmechanismsofwakingstate-dependent\nactivityandalphaoscillations,Neuron110,120–138.e4(2022).\n[39] P.Riyahi,M.A.Phillips,N.Boley,andM.T.Colonnese,Experiencedependenceofalpharhythmsand\nneural dynamics in the mouse visual cortex, Journal of Neuroscience 44, 10.1523/JNEUROSCI.2011-\n22.2024(2024).\n[40] O. Jensen and A. Mazaheri, Shaping functional architecture by oscillatory alpha activity: Gating by\ninhibition,FrontiersinHumanNeuroscience4,10.3389/fnhum.2010.00186(2010).\n[41] M.Halgren,I.Ulbert,H.Bastuji,D.Fabó,L.Erss,M.Rey,O.Devinsky,W.K.Doyle,R.Mak-McCully,\nE.Halgren, L.Wittner, P.Chauvel, G.Heit, E.Eskandar, A.Mandell,andS.S.Cash,Thegeneration\nand propagation of the human alpha rhythm, Proceedings of the National Academy of Sciences 116,\n23772–23782(2019).\n[42] J.J.FoxeandA.C.Snyder,Theroleofalpha-bandbrainoscillationsasasensorysuppressionmechanism\nduringselectiveattention,FrontiersinPsychology2,10.3389/fpsyg.2011.00154(2011).\n[43] Y. Bagherzadeh, D. Baldauf, D. Pantazis, and R. Desimone, Alpha synchrony and the neurofeedback\ncontrolofspatialattention,Neuron105,577–587.e5(2020).\n[44] O.Jensen,Distractorinhibitionbyalphaoscillationsiscontrolledbyanindirectmechanismgovernedby\ngoal-relevantinformation,CommunicationsPsychology2,36(2024).\n[45] A.Aggarwalter A.Aggarwal, C.Brennan, J.Luo, H.Chung, D.Contreras, M.B. Kelz,andA.Proekt,Visualevoked\nfeedforward-feedbacktravelingwavesorganizeneuralactivityacrossthecorticalhierarchyinmice,Nature\nCommunications13,4754(2022).\n[46] D. Shin, K. Peelman, J. D. Rosario, and B. Haider, Narrowband gamma oscillations propagate and\nsynchronizethroughoutthemousethalamocorticalvisualsystem(2022),type: article.\n[47] A. B. Saleem, A. D. Lien, M. Krumin, B. Haider, M. R. Rosón, A. Ayaz, K. Reinhold, L. Busse,\nM.Carandini,andK.D.Harris,SubcorticalSourceandModulationoftheNarrowbandGammaOscillation\ninMouseVisualCortex,Neuron93,315–322(2017).\n[48] M.S.Hoseini,B.Higashikubo,F.S.Cho,A.H.Chang,A.Clemente-Perez,I.Lew,A.Ciesielska,M.P.\nStryker,andJ.T.Paz,GammarhythmsandvisualinformationinmouseV1specificallymodulatedby\nsomatostatin+neuronsinreticularthalamus,eLife10,e61437(2021).\n[49] N.Meneghetti,C.Cerri,E.Tantillo,E.Vannini,M.Caleo,andA.Mazzoni,Narrowandbroadγbandspro-\ncesscomplementaryvisualinformationinmouseprimaryvisualcortex,eNeuro8,10.1523/ENEURO.0106-\n21.2021(2021).\n[50] A.Fernandez-Ruiz,A.Sirota,V.Lopes-dosSantos,andD.Dupret,Overandabovefrequency: Gamma\noscillationsasunitsofneuralcircuitoperations,Neuron111,936953(2023).\n25\n[51] R.G.TownsendandP.Gong,Detectionandanalysisofspatiotemporalpatternsinbrainactivity,PLOS\nComputationalBiology14,e1006643(2018).\n[52] G.Buzsáki,C.A.Anastassiou,andC.Koch,Theoriginofextracellularfieldsandcurrents–eeg,ecog,lfp\nandspikes,NatureReviewsNeuroscience13,407–420(2012).\n[53] A. Speed, J. Del Rosario, N. Mikail, and B. Haider, Spatial attention enhances network, cellular and\nsubthresholdresponsesinmousevisualcortex,NatureCommunications11,505(2020).\n[54] Z.Ye,M.S.Bull,A.Li,D.Birman,T.L.Daigle,B.Tasic,H.Zeng,andN.A.Steinmetz,Brain-wide\ntopographiccoordinationoftravelingspiralwaves,bioRxiv10.1101/2023.12.07.570517(2023).\n[55] G.Chen,M.J.Rasch,R.Wang,andX.-h.Zhang,Experience-dependentemergenceofbetaandgamma\nbandoscillationsintheprimaryvisualcortexduringthecriticalperiod,ScientificReports5,17847(2015).\n[56] M.P.Zimmerman,S.T.Kissinger,P.Edens,R.C.Towers,S.Nareddula,Y.Y.Nadew,C.J.Quinn,and\nA.A.Chubykin,Originofvisualexperience-dependentthetaoscillations,CurrentBiology35,87–99.e6\n(2025).\n[57] S.Arroyo,C.Bennett,andS.Hestrin,Correlationofsynapticinputsinthevisualcortexofawake,behaving\nmice,Neuron99,1289–1301.e2(2018).\n[58] J.D.Semedo,A.I.Jasper,A.Zandvakili,A.Krishna,A.Aschner,C.K.Machens,A.Kohn,andB.M.\nYu,Feedforwardandfeedbackinteractionsbetweenvisualcorticalareasusedifferentpopulationactivity\npatterns,NatureCommunications13,1099(2022).\n[59] R.JordanandG.B.Keller,Opposinginfluenceoftop-downandbottom-upinputonexcitatorylayer2/3\nneuronsinmouseprimaryvisualcortex,Neuron108,1194–1206.e5(2020).\n[60] G.Matteucci,M.Riggi,andD.Zoccolan,Atemplate-matchingalgorithmforlaminaridentificationof\ncorticalrecordingsitesfromevokedresponsepotentials,JournalofNeurophysiology124,102–114(2020).\n[61] L.A.Zhang,P.Li,andE.M.Callaway,High-resolutionlaminaridentificationinmacaqueprimaryvisual\ncortexusingneuropixelsprobes,bioRxiv10.1101/2024.01.23.576944(2024).\n[62] K. S. Rockland, Notes on visual cortical feedback and feedforwardconnections, Frontiers in Systems\nNeuroscience16,10.3389/fnsys.2022.784310(2022).\n[63] T. van Kerkoerle, M. W. Self, B. Dagnino, M.-A. Gariel-Mathis, J. Poort, C. van der Togt, and P. R.\nRoelfsema,Alphaandgammaoscillationscharacterizefeedbackandfeedforwardprocessinginmonkey\nvisualcortex,ProceedingsoftheNationalAcademyofSciences111,14332–14341(2014).\n[64] D.A.PinotsisandE.K.Miller,Beyonddimensionreduction: Stableelectricfieldsemergefromandallow\nrepresentationaldrift,NeuroImage253,119058(2022).\n[65] D.A.Pinotsis,G.Fridman,andE.K.Miller,Cytoelectriccoupling: Electricfieldssculptneuralactivity\nandtunethebrainsinfrastructure,ProgressinNeurobiology226,102465(2023).\n[66] M.Wischnewski,H.Tran,Z.Zhao,S.Shirinpour,Z.J.Haigh,J.Rotteveel,N.D.Perera,I.Alekseichuk,\nJ.Zimmermann,andA.Opitz,Inducedneuralphaseprecessionthroughexogenouselectricfields,Nature\nCommunications15,10.1038/s41467-024-45898-5(2024).\n[67] S. van Bree, D. Levenstein, M. R. Krause, B. Voytek, and R. Gao, Processes and measurements:\na framework for understanding neural oscillations in field potentials, Trends in Cognitive Sciences\n10.1016/j.tics.2024.12.003(2025).\n[68] G.Plomp,C.M.Michel,andC.Quairiaux,Systematicpopulationspikedelaysacrosscorticallayerswithin\nandbetweenprimarysensoryareas,ScientificReports7,10.1038/s41598-017-15611-2(2017).\n[69] B.Teleczuk,N.Dehghani,M.LeVanQuyen,S.S.Cash,E.Halgren,N.G.Hatsopoulos,andA.Destexhe,\nLocalfieldpotentialsprimarilyreflectinhibitoryneuronactivityinhumanandmonkeycortex,Scientific\nReports7,40211(2017).\n[70] G.B.HangandY.Dan,Asymmetrictemporalintegrationoflayer4andlayer2/3inputsinvisualcortex,\nJournalofNeurophysiology105,347–355(2011),pMID:21068267.\n[71] Y.Orsher,A.Rom,R.Perel,Y.Lahini,P.Blinder,andM.Shein-Idelson,Sequentiallyactivateddiscrete\nmodulesappearastravelingwavesinneuronalmeasurementswithlimitedspatiotemporalsampling,eLife\n12,10.7554/elife.92254.3(2024).\n[72] U.Mitzdorf,Propertiesoftheevokedpotentialgenerators: Currentsource-densityanalysisofvisually\nevoked potentials in the cat cortex, International Journal of Neuroscience 33, 33–59 (1987), pMID:\n3610492.\n[73] J.J.Orczyk,A.Barczak,J.Costa-Faidella,andY.Kajikawa,Crosslaminartravelingcomponentsoffield\npotentialsduetovolumeconductionofnon-travelingneuronalactivityinmacaquesensorycortices,Journal\nofNeuroscience41,7578–7590(2021).\nNested spatiotemporal theta–gamma waves\norganize hierarchical visual processing\nBrendanHarris,PulinGong\nNatureCommunications\nR:Reviewers’Comment, A:Authors’Comment\nWethankthereviewersfortheirconstructivefeedback. Wearepleasedthatthereviewersfoundourrevisions\nsatisfactory and that the main concerns from the initial review have been addressed. Below we provide a\npoint-by-pointresponsetotheremainingminorcomments,quotingtherevisedmanuscriptwhererelevant(new\ntextingreen). Linenumbersrefertothetracked-changesversionofthemanuscript.\nReviewer#1\nR1.1: Iappreciatetherevisions. Theauthorshaverespondedtothemainpointsfromthefirstround. Inparticular,\nthey now connect θ activity to spiking in a way that is behaviorally meaningful, they separate aperiodic\n(1/f) from oscillatory effects more clearly, and they give a more coherent picture of θ dynamics across\nonset/offsetandlayers. Thismovesthepaperintherightdirection. Myinitialconcernswerethreefold: (1)\nthefunctionalityofspike–LFPcoupling,(2)theroleofthe1/fcomponent,and(3)theorigin/dynamicsofθ.\nOnallthree,thenewversionisimproved. Icanfollowhowspike–θcouplingrelatestohit/missdifferences;\nthe1/fanalysisistreatedasinformativebutnotthemaindriver;andthetwo-modeθstory(onsetvs. offset)is\nstatedmoreexplicitly. Thatsaid,themanuscriptisstillhardtoread. Theanalysesarecomplexanddensely\npresented. Thisisnotaproblemperse,butitdoesmaketheentrypointsteepevenforreadersfamiliarwith\nthefield. Myrecommendationistoacceptafterminor,editorialrevisionsaimedataccessibilityratherthan\nadditionalanalysis.\nA1.1: Wethankthereviewerfortheirconstructivefeedback. Wehavenowaddressedtheminorsuggestionsbelowand\nmadeeditorialrevisionstoimprovethereadabilityofourmanuscriptforabroaderaudience.\nR1.2: Minorsuggestionsaimingtoimprovereadability: 1)Addabriefroadmapparagraphorschematicupfront\nthatsays,inplainlanguage,whatθdoesatonsetandoffset,whereinthelaminaethishappens,andhowthis\ntiestobehavior.\nA1.2: Assuggestedbythereviewer,wehavenowupdatedthefinalparagraphoftheIntroductiontoincludeabrief\nroadmapofourmainfindingsontravelingθ waves. Specifically,wehaveprovidedadditionaldetailonthe\nlaminarstructureofθpropagationatonsetandoffset,itsbehavioralrelevance,aswellasitslinkstoneuronal\nspikingandγ activity(lines74-85). Wehavealsoaddedanannotationtoourfirstschematicfigure(Fig.1),\nlabeling‘superficial’,‘middle’,and‘deep’layerstohighlightthetranslaminarpropagationofθwaves.\nR1.3: 2)Forthespike-θactivityresults,reportsimpleeffectsizesforphaseshiftsandfiring-ratedifferences. This\nmighthelpreaderstoassessthemagnitudeofthechanges.\nA1.3: Assuggested,wehavenowreportedeffectsizesinthemaintextforspike–θcoupling: bothpreferredphase\ndifferences(lines476-489)andfiringratedifferences(lines490-504). Theenhancedpeakfiringratein\nsuperficiallayersduring‘hit’trials,at∼90msafterstimulusonset,hasarankbiserialcorrelationfrom0.37\n(VISp;moderate)to0.68(VISrl;large). Forthechangeinpreferredphaseinthesuperficialandmiddlelayers\nofintermediateareas,themaximumrankbiserialcorrelationis0.85(VISl,layer4).\neffectsrelativetotheoscillatoryreadouts,plusthebandusedforthefit. Addonecompactcontrolinthe\nSupplement(e.g.,causalfilteringorbipolar/CSDonasubset)toaddresstheusualconcernsaboutzero-phase\nfiltersandvolumeconduction. Bottomline: myearlierconcernsareaddressed. Themainremainingissue\nispresentation. Withthesmallchangesabove,Ifeelthepaperwillbemucheasiertofollowforabroader\naudiencewithoutdilutingthetechnicalcontent.\nA1.4: Tomakeourpaperaccessibletoabroaderaudience,wehavenowpresentedthe1/f resultmoresuccinctlyand\nclarifythefrequencybandusedforfittingtheFOOOFmodel(lines197-206).\nAssuggestedbythereviewer,wehavealsoaddedacausalfiltercontrol(forward-only,8th-orderButterworth),\nwhichshowshighlysimilarpatternsoftranslaminarandinterarealθpropagationdelayedby∼150msduetothe\nexpectedfilterphaseoffset(Fig.S15). Thiscontrolindicatesthepropagationpatternsarenotartifactsofthe\nzero-phasefiltering. Furthermore,thesimilartimingofpeaksinθpropagationcoherenceforbothawideθband\n(3–10Hz;Fig.3)andnarrowerband(3–5Hz;Fig.S13)suggeststhatthebidirectionalpatternsarenotdueto\nbandpassringing.\nR1.5: Reviewer #1 (Remarks on code availability): I conducted a static audit of the Code Ocean capsule. The\nrepositoryiswell structured,andtheREADMEcoversinstallation,runorder,andhowtoreproducethemain\nanalysesandfigures. TheJuliaenvironmentispinnedviaProject.tomlandManifest.toml,withtheJulia\nversionstated,whichshouldkeepdependencyresolutionstable. Therearescript-basedentrypoints(notonly\nnotebooks),andthefigure-generationstepsaredescribed. Dataavailabilityisdocumentedandconsistent\nwithNaturePortfolioguidance. Ididnotrunthefullworkflowduetotimeconstraints,butthedefinedrun\npipelineandversionedcomputeenvironmentgivemereasonableconfidencethatotherscanreproducethe\nresults. Overall,thecodeisausablecommunityresource: licensed,documented,portable,andwithclear\nreproductionsteps. Minorsuggestions: makerandomseedsandthreadingdefaultsexplicitintheREADME,\nandincludeaverysmall\"quick-start\"target(e.g.,onefigureonatinysample)thatfinishesinafewminutes.\nA1.5: We thank the reviewer for evaluating our code. We have now added two notes to the README file of the\naccompanyingsoftwarerepository,specifyinghowrandomseedsaresetandthedefaultparallelizationbehavior.\nWehavealsoincludedashortquick-startsectiontothebeginningoftheREADME,whichguidesusersthrough\ntheessentialstepstoloaddataforasingletrialandgenerateFig.1ofthemanuscript.\nReviewer#2\nR2.1: Wethanktheauthorsforthoroughlyaddressingmypreviouscomments. Themanuscriptimprovedsignifi-\ncantly.\nA2.1: Wethankthereviewerfortheirconstructivefeedbackthroughoutthereviewprocess.\nReviewer#3\nR3.1: Theauthorshaveaddressedallthemajorandminorconcerns,andtheclarityofthemanuscripthasimproved.\nTheadditionalanalysesprovidemoreconfidenceinthefindings,andtheresultsareinteresting,novel,and\nimpactful. Westillhavesomeminoradditionalqueriesandsuggestions:\nA3.1: We thank the reviewer for their constructive feedback on our revisions. We are pleased that the additional\nanalyseshavestrengthenedourfindingsandimprovedthemanuscript’sclarity. Below,weaddresseachofthe\nremainingqueriesandsuggestions.\nbegoodtoindicatethehierarchyoftheareasontheright.\nA3.2: Assuggestedbythereviewer,wehaveaddedasmallschematictoFig.4eandFig.4ftoindicatethehierarchical\norderoftheareasattheright.\nR3.3: The order parameter change analysis (figure 3, lines 293-305) shows this value starts decreasing before\nstimulusonset. Pleaseclarifywhythismightbeobserved. Forexample,couldthisbeatrivialeffectofthe\nanalysismethod,duetoacontinuouswavepropagationpatterninsomeareas,orreflectentrainmentoflocal\nnetworkstostimuluspresentations?\nA3.3: We have now clarified, at lines 242-268, that the pre-stimulus dip in the order parameter is an expected\neffectofnon-causalleakagefromthebandpassmethodusedtoextracttheθcomponent. Thiszero-phasefilter\npreservespeaktiming,butbroadenssharppeaksandtroughsforwardandbackwardintime. Toexcludespurious\nphaserelationships,werepeatedtheanalysiswithacausal,forward-only8th-orderButterworthfilter,asshown\ninFig.S15. Thecausalfilterresultsinsimilartranslaminarandinter-arealθasthenon-causalfilter,butshifted\nbyanexpecteddelayof∼150ms,indicatingthattheobservedbidirectionalpropagationisnotafilterartifact.\nAdditionally,filterchoicedoesnotinfluencethetimingofθpropagationcoherence;peaksinthewavenumber\nandorderparametersoccuratsimilartimesforbothwider(3–10Hz;Fig.3)andnarrower(3–5Hz;Fig.S13)θ\nbands.\nR3.4: InFigure4,whilethefunctionalandanatomicalhierarchiesshowdistinctpatternsofconnectivity,propa-\ngationcoherencepatternacrossthesedifferenthierarchiesisquitesimilarintime. Doesthissuggestthat\nthepatternofwavepropagationisstructureinvariant? Isthispatternpresentinspontaneousand/orflashed\nimages?\nA3.4: Thanksforraisingthispoint. Wehavenowclarifiedinthemaintextthatthesimilarpatternsofpropagation\nalong the anatomical and functional hierarchies is expected, since these hierarchies are strongly correlated\n(Siegle et al.). Thus, the resemblance does not imply structure invariance; rather, both hierarchies impose\na comparable ordering of areas, yielding similar temporal profiles of θ propagation. Regarding condition\ndependence,thebidirectionalinter-arealpropagationduringthetaskcondition(Fig.4)isabsentunderflash\nstimulationwithouttaskengagement: repeatingthisanalysisforflashedimagesrevealednosignificantepochs\nofinter-arealpropagation,asshowninanewsupplementaryfigure(Fig.S16).\nR3.5: It is claimed that PAC is strongest in VISp superficial layers during stimulus presentation, but only the\ncomodulogramforVISlinthespontaneousconditionisshown. ThisisconfusingbecausewhileVISLshows\ncleartheta-gammacouplingintheabsenceofvisualstimuli,neitherVISpnorVISlshowclearcouplingin\nanyothercondition. Basedonthecontourmaps,itseemsthatthestrongestcouplingwouldbearound4\nHzand25-20Hz. FigS10alsoshowsthatothersecondaryvisualareashavetheta-gammacouplinginall\nconditions. Forconsistency,wesuggestchangingthecomodulogrampanelinFig5.\nA3.5: We thank the reviewer for pointing out this potential source of confusion. We agree that a single example\ncomodulogramcannotcapturethefullrangeofPACpatternsobservedacrossareasandstimulusconditions.\nToimprovetheclarityofourmanuscript,wehaveremovedtheexamplecomodulogramfromFig.5,focusing\nonFigs.5cto5fthatpresentourmainresultsonspatialandspatiotemporalvariationinPACstrength. For\ncompleteness,wehavenowhighlightedalltypesofcomodulogrampatternsobservedacrossareasandstimuli(at\nlines415-426ofthemaintext)andwedirectreaderstoFig.S9forthefullsetofarea-andcondition-specific\ncomodulograms.\nR4.1: Ico-reviewedthismanuscriptwithoneofthereviewerswhoprovidedthelistedreports. Thisispartofthe\nNatureCommunicationsinitiativetofacilitatetraininginpeerreviewandtoprovideappropriaterecognition\nforEarlyCareerResearcherswhoco-reviewmanuscripts.\nA4.1: Wethankthereviewerfortheircontributionstothereviewprocessandfortheirco-reviewefforts.\nReviewer#5\nR5.1: Theauthorsmadeasignificantefforttoaddressmyconcerns/comments. Belowisabriefcommentrelatedto\nanotherreviewer’scomment. RegardingA1.4,theauthorshavemadeanefforttoruleoutthepossibilityof\nvolumeconductionfromthehippocampus. However,thepossibilityisnotcompletelyruledoutbythetwo\nreasonsIexplainbelow. Havingsaidthat,asReviewer1suggests,itisdifficulttoruleoutthispossibility\ncompletely. Thus,pleasetakethefollowingcommentsjustasanote,notnecessarilyimpactingthedecision\ninhandlingthemanuscript. (1)Differentphaseoffsetsindifferentcorticesindifferenthierarchiescanbe\ncausedbythenatureofthethetaoscillationtravelinginthedorso-ventralaxisinthehippocampus. Thus,\ndifferentphasesacrosscorticalareasmaybesimplyexplainedbythephaseoffsetatthelevelofthesubcortical\nsource(e.g.,thehippocampus)(2)Toeliminatethispossibility,itisessentialtoverifythecross-layerphase\nrelationshipwithinthesamecortex,ratherthanacrosscorticalareas,assuggestedbyReviewer1.\nA5.1: We thank the reviewer for this comment, and have now incorporated this point into the Discussion section\nof the revised manuscript (lines 630-646). Specifically, we have emphasized that the cross-layer phase\nrelationshipsobservedwithinindividualareasareinconsistentwithθconductedfromthehippocampus. We\nhavealsohighlightedthatfutureexperimentsshouldjointlyrecordtravelinghippocampalandcorticalθpatterns\ntodefinitivelyaddressthevolume-conductionhypothesis.\n4"}]

Decision letter