<figure><img src="image_1.png" /></figure>

HAL
open science

# A hierarchical multiscale model of forward and backward alpha-band traveling waves in the visual system

Jakob C B Schwenk, Andrea Alamia

► To cite this version:

Jakob C B Schwenk, Andrea Alamia. A hierarchical multiscale model of forward and backward alpha-band traveling waves in the visual system. PLoS Computational Biology, 2025, 21 (8), pp.e1013294. 10.1371/journal.pcbi.1013294. hal-05365775

HAL Id: hal-05365775

https://hal.science/hal-05365775v1

Submitted on 14 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

# A hierarchical multiscale model of forward and backward alpha-band traveling waves in the visual system

Jakob C. B. Schwenk*, Andrea Alamia

Centre de Recherche Cerveau et Cognition (CerCo), CNRS, Université de Toulouse, Toulouse, France

* jakob.schwenk@cnrs.fr

<figure><img src="image_4.png" /></figure>

## Abstract

Recent studies have shown that cortical low-frequency oscillations are often organized as traveling waves. The properties of these waves have been linked to both sensory processing and cognitive functions. In EEG recordings, alpha-band (~10Hz) traveling waves propagate predominantly along the occipital-frontal axis, with forward waves being most prominent during visual processing, while backward waves dominate at rest and during sensory suppression. While a previous study has proposed a functional model to explain their generation and propagation, a biologically plausible implementation is lacking. Here, we present a multi-scale network model with mean-field dynamics that, building on known cortical connectivity, reproduces the dynamics of alpha-band traveling waves observed in EEG recordings. We show that forward and backward waves can arise from two distinct cortical sub-networks that are connected in infragranular layers at each area. At rest, the network generates spontaneous backward waves and switches to a forward state upon sensory stimulation, reproducing the dynamics observed in EEG recordings. We then show that a cortico-thalamic pathway through the pulvinar can bias the dynamics to the forward state and that pulvinar engagement leads to spontaneous forward waves at rest. This is in line with previous studies suggesting a key role for the pulvinar in directing cortical information flow. In summary, our model provides a biologically plausible architecture for modeling the dynamics of macroscale traveling waves. It bridges the gap between scales by connecting laminar activity to scalp-level spatial patterns, providing a biologically grounded and comprehensive view of the spatial propagation of alpha-band traveling waves.

<figure><img src="image_5.png" /></figure>

OPEN ACCESS

**Citation:** Schwenk JCB, Alamia A (2025) A hierarchical multiscale model of forward and backward alpha-band traveling waves in the visual system. PLoS Comput Biol 21(8): e1013294. https://doi.org/10.1371/journal.pcbi.1013294

**Editor:** Daniele Marinazzo, Ghent University, BELGIUM

**Received:** February 12, 2025

**Accepted:** July 2, 2025

**Published:** August 11, 2025

**Copyright:** © 2025 Schwenk, Alamia. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data availability statement:** The custom code used to generate the model and run the simulations described in this paper is publicly available online (github.com/jcbschwenk/multiscale-tw-model). The human EEG data used to validate the model behavior is available at osf.io/m9s2j (DOI 10.17605/OSF.IO/M9S2J).

## Author summary

The human EEG shows strong alpha (7–13 Hz) activity over visual areas, which has recently been shown to propagate as traveling waves across the scalp.
**Funding:** This project was funded by the European Union under the European Union's Horizon 2020 research and innovation program (grant agreements No. 101075930 to AA). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

The direction of these waves strongly correlates with visual input and attention. While invasive studies have also reported traveling waves at the cortical level, it remains unknown how these two levels of observation are linked. We present a physiologically plausible model that bridges this gap by linking scalp-level wave patterns to mean-field activity at the laminar level in the cortex. In addition to explaining previously published EEG results, our model makes key predictions to be tested in future experiments. Overall, this study contributes to the interpretability of EEG investigations of traveling waves.

## Introduction

Oscillatory activity is a ubiquitous feature of neural processing throughout the brain and has been linked to a broad array of sensory, cognitive, and motor functions [1]. The most dominant temporal frequency band in the visual system is the alpha rhythm (approx. 7 – 13 Hz), which has been linked to several cognitive processes. For example, previous studies considered alpha oscillations as a rhythm mainly involved in modulating top-down inhibition in sensory-specific cortical areas [2–5]. However, due to its strong association with sensory input and phasic influence on perception, alpha oscillations have also been proposed as an internal temporal reference frame for neural processing in the visual system [6–8]. More recently, as studies on neural oscillations shift their focus also to consider their spatial dimension, it was demonstrated that alpha oscillations propagate as traveling waves, as measured at the cortical surface [9], as well as using EEG [10–12]. At the scalp level, the propagation direction of these waves is concentrated mainly on the anterior-to-posterior axis, with waves traveling either forward (FW, i.e., towards anterior sensors) or backward (BW, i.e., in the opposite direction). At rest, as well as during attentional suppression of visual input, the dominant direction is BW, while FW waves are more prominent during visual stimulation and attention [10–14]. These findings and the previous literature on alpha oscillations suggest that alpha-band traveling waves may play a functional role in the visual system. Yet the cortical spatiotemporal patterns of activity that correspond to FW and BW waves at the scalp level remain unclear. On the one hand, from an experimental point of view, some studies have attempted to infer the cortical sources using M/EEG signals [e.g., 11,15], but the precision of these methods remains limited and, ultimately, invasive, multiscale recordings (i.e., simultaneous from the scalp and cortex) may be needed to gain a full understanding of the waves’ cortical origins. An increasing number of studies reporting traveling waves (including in the alpha band) from intracortical and cortical surface recordings [9,16–18] further highlight the need to bridge the gap between this literature and the EEG findings. From a computational perspective, a relatively simple hierarchical network model could predict the FW/BW waves’ direction in EEG recordings [19]. More specifically, in this model, based on the theoretical framework of predictive coding [20], each area continuously predicts the activity of the hierarchically lower area. Here, the difference between the current prediction and the actual activity determines the prediction error,
which is sent back to update the next prediction. This dynamic between regions generates alpha traveling waves, which arise naturally from the resulting inter-areal feedback loop with plausible delays. The FW or BW direction of propagation is determined by the presence (FW) or absence (BW) of the input to the first stage of the hierarchy. While this model offers a computational framework for investigating traveling waves at the scalp level, it fails to provide a physiologically plausible architecture that describes the mechanisms involved in generating alpha traveling waves at the level of cortical neural circuits.

In this study, we aimed to bridge this important gap in the literature and provide a biologically plausible model to investigate the cortical generation and propagation of alpha-band traveling waves observable at the scalp level. In implementing our model, we attempted to integrate a hierarchical architecture, as proposed by previous computational studies and theoretical work, into the known connectivity of the visual cortex and its features. First, we designed a model such that the interlaminar and interareal circuits generate forward and backward waves through different cortical pathways and validated the results of our model with human EEG recordings. We then explored how the thalamus could modulate these waves. Specifically, we considered the pulvinar, which is a critical, higher-order thalamic nucleus with projections both from and to the cortex [21], and also a known generator of alpha-band oscillations [18]. Importantly, several studies point to the pulvinar as being involved in the modulation of cortico-cortical communication, including through oscillatory dynamics [22–26]. This involvement suggests the pulvinar may also be crucial in generating alpha-band traveling waves, and we explored its effect on cortical waves in our implementation.

All in all, the resulting network allows us to accurately describe the spatial propagation of alpha traveling waves along the cortical hierarchy, faithfully replicating the propagation patterns observed in EEG recordings. Our results reveal that large-scale FW and BW traveling waves can arise from specific and distinct cortico-cortical connections, that can be modulated by the pulvinar and that their propagation direction is determined by the presence of sensory input.

## Results

We designed a model reproducing the spatial patterns of alpha traveling waves, as observed in EEG recordings, during rest, and in response to visual stimulation. Our implementation is grounded in a few well-known visual system features. First, the visual cortex is organized hierarchically, with largely separate feedforward and feedback pathways [27–29]. Another prominent feature is its laminar organization, reflecting both in the interareal flow of feedforward and feedback signals and in the circuit within each area where these signals are combined [30,31]. In addition, there is also strong evidence for a laminar separation of oscillatory sources, as alpha generators have been specifically localized in deep (infragranular) layers [32–34]. We considered all of this evidence to constrain the cortical connectivity in our model. In the following, we will describe the model’s structure, characterize its general response behavior, and compare it to experimental data from human EEG recordings during rest and visual stimulation. We then continue with a more detailed analysis of the network’s pathways that generate FW and BW waves. Finally, we explore how the model’s dynamics change by adding a cortico-thalamo-cortical pathway through the pulvinar, investigating how thalamic engagement could reproduce the effects of attentional modulation.

## Model architecture

Our model represents a generalized description of hierarchical processing streams in the visual system without modeling specific cortical areas. The general architecture is illustrated in Fig 1a. It describes a set of $N_{Cx}$ hierarchically connected cortical areas with identical intra- and interlaminar connectivity. Each area comprises a set of nodes with mean-field dynamics across three laminar compartments: the supragranular layers 2/3, the input layer 4, and the infragranular layers 5/6. Incoming feedforward signals enter each area at $L_4$ ($L_4_x$), where initial response adaptation is applied through a local inhibitory population ($L_{4_{IN}}$), and are then passed to the supragranular layers ($SG_x$) within the same area. From here, projections extend as feedforward connections to the next area and locally to infragranular layers ($IG_B$). The infragranular
<figure><img src="image_6.png" /><figcaption>**Fig 1.** Model overview and basic response behavior. A: Modular architecture. The model comprises a laminar network in a hierarchical stream of cortical areas Cx1 to CxN. Each area sends feedforward signals from supragranular layers and receives feedback to both laminar compartments from the infragranular layer of the higher area. Nodes in infragranular layers are pacemakers generating a continuous alpha rhythm. B: Simulated responses and wave fits for a single model run with NCx=3 areas during rest and visual stimulation; (1) Stimulus profile (1.2 nA DC step function); (2, 3) cortical activity in supra- and infragranular layers (colors represent different cortical areas); (4) average of all layers' activity in each cortical region at the source- and (5) at the scalp-level (via forward projection to simulate an EEG signal); (6) traveling wave direction fit: the time-course shows the angle of the best fit to the EEG phase gradient (see Methods for details); shaded areas mark periods of a stable classification into either BW (red) or FW (blue) state; (7) topographical plots showing the mean EEG phase gradients (across 50 trials) corresponding to the classified BW and FW periods. C: Left: Response curve showing the probability of a FW state classification during a 3 sec DC stimulation as a function of input current. The dashed line marks the estimated threshold current (center of the sigmoid) at which the system switches wave states. Right: Temporal response functions for FW and BW responses to suprathreshold stimulation; the stimulus sequence used was a random step function with a fixed autocorrelation of 100 ms (~ 1 cycle). No rhythmicity was introduced by the stimulus sequence. Chance-level correlations (from random trial shufflings) are shown with 95% confidence intervals in the shaded area.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1013294.g001

node in each area ($IG_{IB}$) is modeled as a population of intrinsically bursting neurons that generate a steady oscillation at ~9 Hz at baseline, based on neurophysiological evidence demonstrating the presence of strong alpha pacemakers in infragranular layers [35,36]. Even though the generation of alpha rhythms may indeed be more strongly driven by larger, network-scale interactions (see, e.g., [37] for an implementation of Wilson-Cowan dynamics in a laminar model, and [38] for a thalamic model) we opted for a simplified model to focus our investigation on the traveling waves dynamics, i.e., alpha propagation at the cortical level. We will discuss some limitations of this approach below. The $IG_{IB}$ nodes are unidirectionally connected through excitatory feedback between areas. For its role in generating BW-directed traveling waves in the network (see below), we label this stream of backward-coupled $IG_{IB}$ nodes the *BW-pathway*.

Feedforward projections between areas extending from supragranular layers are inhibited through feedback from the subsequent higher area's infragranular layers ($IG_{IB}$). In isolation, this circuitry constitutes a closed feedback loop that effectively reduces feedforward signals to the residuals between areas: we based this implementation on predictive coding principles which, as developed by Alamia & VanRullen [19] with a much simpler approach, explains the emergence of
traveling waves through feedback-loop connections. We demonstrate the relationship between our model and previous work in S1 File. The circuit described by this feedback loop, which carries the feedforward sensory input (as residuals or prediction-errors), will be denoted as the *FW-pathway*. We will discuss in greater detail how each pathway contributes to the oscillatory dynamics of the network below.

In the following, we will present a version of the model with $N_{Cx}=3$ cortical areas, unless stated otherwise. We chose this mainly for model parsimony, and because this model (given identical neural dynamics between areas) best approximated the scalp-level phase gradients observed in real EEG data (see below). However, our main findings generalize to versions with $N_{Cx} > 3$, and our model is not intended to replicate a specific pathway or set of areas.

## Traveling wave dynamics in the cortex

We first investigated the behavior of the network to sensory stimulation, which we model by applying an external current to the lower end of the hierarchy, thus simulating a constant visual stimulus.

**Fig 1b, panels 1–4** summarize the network’s response during rest and to 1 second of direct current (DC) stimulation. In the absence of sensory stimulation, the intrinsically bursting nodes in the infragranular layers ($IG_{IB}$) drive the network activity. The intrinsic rhythmic activity in these nodes synchronizes between areas due to their backward coupling in the network. The fixed inter-areal delay ($\Delta T = 12$ ms) leads to a consistent phase-gradient corresponding to a backward (BW) traveling wave, propagating from higher to lower areas, which is visible in the map of mean-field activity (Fig 1b, panel 4, showing averages across all laminar compartments per cortical area). At the onset of sensory stimulation, the rhythmic activity persists but switches to an opposite phase gradient, now corresponding to a forward (FW) wave, i.e., propagating from lower to higher areas. The activity traces (Fig 1b, panels 2 and 3) show that the rhythmicity during stimulation is carried by both superficial and deep layers, in contrast to the baseline state in which superficial activity is largely absent. Note here that rhythmicity in both states (BW and FW) is generated within the network and never introduced by the stimulus. In a section below, we explore the extent to which the oscillation in the FW state is differentially determined by intrinsic $IG_{IB}$ rhythmicity vs. the inter-areal feedback loop between supragranular nodes.

## Model dynamics at the scalp level

Next, we aimed at modeling the dynamics of traveling waves at the level of the scalp. To this end, we evaluated our model’s mean-field output by simulating an EEG signal using a standard forward model (with sources for $C_{x1}$, $C_{x2}$, and $C_{x3}$ positioned bilaterally in occipital, parietal, and frontal regions, respectively). This additional step allowed us to quantify traveling wave direction and directly compare our simulations with experimental data. Fig 1b, panels 5 and 6 summarize the results of this procedure for the same simulated trial shown in the panels above. The pattern of BW and FW waves remains visible in the spatiotemporal map of the simulated EEG (shown for electrodes on the mid-line, Fig 1b, panel 5). We used an iterative plane fitting procedure (based on [9], see Methods) to obtain the one direction of propagation that best explains the spatial phase gradient for each time point (Fig 1b, panel 6). The final classification into FW and BW states based on this value (shaded red and blue areas, with unshaded areas corresponding to states not classified as either FW or BW) matches the stimulus-dependent reversal of states that is visible in the mean-field activity (i.e., at the source level in the model, Fig 1b, panel 4). Finally, the topography plots in Fig 1b, panel 7 show the mean alpha-phase gradients (averaged across 50 model runs) corresponding to each state, limited to the electrode-ROI used for the fit, revealing a clear gradient of phases, reversed in the two states. In the following, we will primarily use the final classification output (i.e., FW or BW states) to explore the behavior of our model. It should be noted here that the two states are largely, but not strictly, complementary, as there may be cases where fits could not be classified into either category (i.e., the fit produced a correlation coefficient below a given null distribution, see Methods).
The single-trial example in Fig 1b already shows that the network's wave state depends largely on the stimulus current injected, which is in line with experimental observations and previous modeling efforts [19]. To quantify the relationship between forward waves and visual stimulation, we mapped the probability of a FW state (across runs of 3 seconds of DC stimulations) to the stimulus amplitude, which revealed a sigmoidal response curve (Fig 1 c, left). This confirms that the network state depends systematically on the amplitude of the input current and allows the estimation of the threshold current at which the system switches wave states (at 0.88 nA in our simulations, dashed line).

Next, we aimed to characterize the network's generalized impulse response, i.e., the model's response to a brief (impulse) input. Specifically, we were interested in the time constants of the reversal response. To estimate this, we presented the model with suprathreshold random stimulation. Specifically, the stimuli were temporal sequences alternating randomly (non-rhythmic) between ON- and OFF-stimulus periods, with a fixed stimulus current for the ON-periods (above the reversal threshold) and an update rate of 10 Hz, resulting in an autocorrelation of 100 ms (corresponding to approx. one cycle of the network's oscillation, i.e., the minimum resolution imposed by our wave fitting procedure). The resulting response functions for the classified wave states (obtained by cross-correlation with the stimulus sequence; Fig 1 c, right) exhibit a sharp response onset and a decay over a few oscillatory cycles, returning to chance-level around 700ms. This shape demonstrates two crucial properties of the model: 1) the reversal from the BW- to the FW-state is immediate, i.e., the first cycle after stimulus onset already propagates as a FW wave - notably, this mirrors the spatiotemporal pattern of physiological evoked responses to visual stimulation; and 2) the return to the baseline is comparably slow, indicating the network's tendency to remain in the FW state for a few cycles once it has been evoked.

## Comparison of the model to human EEG data

After establishing the main properties of our model, we then compared its behavior to experimental observations, using EEG data published previously [14]. In that study, human subjects were presented with a static luminance patch presented centrally in regular stimulus-ON/OFF sequences (of 5 seconds each). The EEG results showed a consistent increase of FW- and a decrease of BW-wave strengths during the stimulus-ON period (for a detailed report on the quantification of waves in that analysis, cf. their Methods). We considered this dataset for its relatively straightforward design in assessing the effect of visual stimulation on alpha waves. However, similar alpha waves dynamics, in which the direction of propagation depends on the presence of the stimulus, were obtained in other datasets with larger sample size [13,19]. We first replicated this pattern of findings by applying the same phase-gradient analysis that we use to evaluate our model output (cf. Fig 1b, panels 6 and 7, and Methods). The left panel in Fig 2a shows the mean probability of observing a stable FW or BW state over time across trials and subjects. Here and in all other panels, the grey shaded area represents the time-period the stimulus is on. As shown in the right panel of Fig 2a, our model captures the empirical response pattern well, not taking into account the obvious difference in absolute probability values. Note that the aim of our simulations was not to reproduce the exact signal-to-noise (SNR) ratio of the actual EEG data but the relative change over time, and the difference in scaling for both classification probabilities (FW and BW) is an expected consequence of such a difference in SNR. As expected from the impulse response function, our network exhibits the same reversal from a BW wave state during baseline to a FW wave state upon stimulus onset, as present in the experimental data. The transition times for stimulus on- and offset were comparable between experimental and simulated data, accounting for a fixed processing delay in the physiological data that offsets the temporal response of the model of around 300ms with respect to the experimental data, possibly due to processing in the retina and subcortical nuclei (time to first peak at stimulus onset - data: 363 ms, model: 61 ms; return to baseline (corrected for linear slope) - data: 744 ms, model: 421 ms). Thus, in summary, our model can qualitatively reproduce the population-average dynamics of traveling wave direction to visual stimulation.

While the aim of our modelling was on wave direction, we also examined the propagation velocities (spatial frequencies) of the waves generated in our model. Indeed, the phase velocity predicted for the BW state matches those obtained for the real EEG data (~3.7 rad phase difference between electrodes Oz and Fz), while FW wave velocity was
<figure><img src="image_7.png" /><figcaption>**Fig 2.** Comparison of our model to experimental EEG data (from [14]). The stimulation consists of a 5 sec static input (shaded area) followed by a rest period. A: Left: Grand average time-courses of mean FW/BW probabilities from the fit on the experimental data. Shaded areas represent +/-1 SD to indicate between-subject variability. Right: Probability time-courses obtained from 100 model runs (trials). Here, shaded areas represent the 95% CI based on trial-variability. Stimulus amplitude was 0.88 nA (threshold value obtained from the response curve, cf. Fig 1c). B: Experimental data: responses from two subjects, shown as the log-ratio between wave states over time. Variability in the response shape is visible in the baseline offset and the amount of response adaptation after the initial peak. C: We explore two parameter variations to model the variability observed between subjects. FW/BW ratio are shown as a function of different parameter values. Left panel: the weight of inhibitory recurrent feedback on the input in Layer 4 (wL4IN→L4X); Right panel: the weight of excitatory feedback between infragranular pacemakers (wIGIB→IGIB). All weights were varied together for all areas in the network.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1013294.g002

overestimated, most likely as a result of our non-realistic EEG projection (S2A Fig). At the source level, both supra- and infragranular phase gradients (respectively, during stimulus ON and OFF states) extrapolate linearly when more areas are added to the hierarchy (S2B Fig). While the model version we used for our EEG projections (N<sub>Cx</sub>=3) fits the real EEG data best, this should not be taken as an indication for the actual number of areas involved in generating the waves, as our model assumes identical neural time-constants for all areas, which is not true in the visual hierarchy [39,40].

A closer analysis of the data from Pang et al. [14] reveals considerable inter-individual differences in the wave direction responses. Fig 2b shows the response time courses for two subjects, represented here as the log ratio of FW and BW state probabilities. These illustrate the two primary sources of variability in the data we consider relevant for our model: 1) the baseline ratio between the two states (i.e., the individual strength of the bias towards BW waves in the absence of stimulation), and 2) the temporal evolution, i.e., shape of the response within the 5 sec of the stimulus-ON state. We turn to the former in more detail below, where we show that the extension of the cortical model by a pathway through the pulvinar reproduces such bias in the network resting state. Concerning the variability in our model's response shape, we explore its dependence on two key parameters. The first is the amount of suppression exerted on the feedforward signal through neural adaptation at each area. In our model, adaptation is achieved through recurrent coupling to a local inhibitory population in Layer 4 (i.e., L4<sub>IN</sub> in Fig 1a). This allows the strength of input adaptation to be varied directly through the weight of the inhibitory connection. The left panel in Fig 2c shows the effect of varying this parameter on the network response (expressed as the log-FW/BW ratio, as in Fig 2b), given the same stimulus pattern presented before (i.e., 5 seconds of DC stimulation). As expected, the evoked FW response is increasingly suppressed with stronger adaptation, while
the first cycle of the response remains unaffected. The dynamics of this rapid decay follow directly from the time constant of the Layer 4 recurrent inhibition (neural integration time + synaptic time constants).

The second parameter we varied is the strength of the excitatory infragranular feedback connection (IG<sub>IB</sub> in Cx<sub>i</sub> to IG<sub>IB</sub> in Cx<sub>i-1</sub>) to determine whether the network's wave state can also be biased through modulation of the BW pathway (higher weights correspond to stronger feedback). The results are shown in the right panel in Fig 2c. The synaptic connection between IG<sub>IB</sub> nodes is excitatory, thus higher weights correspond to stronger feedback. The effect of feedback weight largely mirrors the one obtained for the input adaptation strength, with stronger feedback leading to a weakened sustained FW response. This result illustrates the continuous competition between the FW and BW pathways, which will be discussed below.

These results demonstrate that the response dynamics can already be flexibly modulated through the variation of two synaptic weights in the network, possibly explaining some inter-subject variability in EEG recordings. Naturally, other parameters may play a role in modulating the pattern and strength of traveling waves, such as the weight of the feedforward connection (SG<sub>x</sub> → L4<sub>x</sub>), or the di-synaptic inhibitory feedback to superficial layers (IG<sub>IB</sub> → SG<sub>IN</sub> → SG<sub>x</sub>). Besides explaining inter-subject variability, differences in specific parameters between cortical areas may be instrumental for the modeling of specific functional networks, e.g., a systematic increase in neural timescales along cortical hierarchies to match empirical observations [39,40], or in formulating specific predictions about alteration in traveling wave patterns in clinical populations [41].

## Subcomponents of the network: FW and BW pathways

Next, we investigate the mechanisms that drive the generation and propagation of FW and BW waves in our model. As described above, the network architecture can be separated into two components: 1) the BW-pathway, a set of backward coupled nodes with intrinsic rhythmicity in the infragranular layers (IG<sub>IB</sub>), and 2) the FW-pathway, an inter-areal closed loop in which signals are passed forward from superficial layers into Layer 4 (SG<sub>x</sub> → L4<sub>x</sub>), then locally through superficial to deep layers, and finally as feedback into superficial layers of the lower area (IG<sub>IB</sub> → SG<sub>IN</sub>). This implementation is inspired by previous modeling and theoretical work considering a predictive coding perspective [19,42]. In both pathways, rhythmic activity emerges in response to non-rhythmic (constant) external stimulation. In the following, we will examine how the rhythmicity of the network (specifically, its dominant temporal frequency) is influenced by the mechanisms generating oscillations in either pathway. Fig 3a shows the distribution of spectral power in supra- and infragranular nodes at baseline (gray lines) and during sensory stimulation (colored). As illustrated by the single-trial example in Fig 1b, the intrinsically bursting property of the IG<sub>IB</sub> nodes leads to a constant stable rhythm in the infragranular layers that increases in frequency only slightly during stimulation (see below). In contrast, rhythmic SG<sub>x</sub> activity only emerges during stimulation, matching the frequency of the IG<sub>IB</sub> pacemaker. Mapping the peak frequency of activity in both compartments against the input current (Fig 3b) confirms this pattern. Fig 3b also reveals that the two frequencies are tightly coupled, showing a slight increase in frequency in both IG<sub>IB</sub> and SG<sub>x</sub> following an increase in input current above threshold (0.62 Hz and 1.30 Hz increase over the suprathreshold range for IG<sub>IB</sub> and SG<sub>x</sub>, respectively). To disentangle the influences of the two pathways on the network's temporal frequency, we leverage the fact that the mechanisms generating oscillatory activities within each pathway are different. While rhythmicity in the BW pathway is a property of the nodes and thus independent of connectivity, it emerges through a feedback loop in the FW-pathway. In the latter, the oscillation frequency depends (in addition to neural integration times), crucially, on the inter-areal delay (for a detailed derivation of this, we refer the reader to [19]). Building on this dependency specific for the FW, but not the BW pathway, we use a systematic variation of the inter-areal delay to disentangle the two pathways during visual stimulation. To identify the interaction between the two pathways, we compare the changes in the dynamic between two different model configurations while varying this parameter (as summarized in Fig 3c and 3d). In the first configuration, we test the full model with both pathways in place (left panels in Fig 3c and solid lines in Fig 3d). In the second configuration, the FW-pathway has been isolated by removing
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Fig 3.** Isolation of FW- and BW-pathways in the model. A: Mean power spectra of all $IG_{IB}$ (top) and $SG_x$ activity (bottom) during stimulus On and Off periods. B: Peak frequency for both node types as a function of input current. The two frequencies are tightly coupled for stimulation above the threshold. C: Variation of inter-areal delay in two configurations. The left panels show the full model; in the right panels, the FW pathway has been isolated by reducing the infragranular nodes to a passive relay (without intrinsic rhythmicity). The bottom panels show the mean power spectra for $SG_x$ as a function of delay. D: Peak frequencies for the data presented in C. The two pathways diverge for larger delays in the full model, while the intrinsic frequency of the FW pathway is consistently lower. Shaded areas in A, B & D show 95% quantiles for the variability between model runs.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1013294.g003

both the intrinsic rhythmicity and the inter-areal feedback in the BW-pathway (right panels and dashed line). As predicted, in both model versions, increasing the inter-areal delay caused a decrease in the dominant frequency of $SG_x$ activity (Fig 3c, lower panels). In the full model, this decrease leads to a divergence between $SG_x$ and $IG_{IB}$ peak frequencies (Fig 3d) for delays larger than approximately 20 ms. In this diverging regime, oscillatory power in the $SG_x$ nodes also decreases due to the loss of synchronization between the two pathways. This desynchronization removes the IG pacemakers' facilitating influence on the FW-pathway activity. For larger delays, $SG_x$ oscillatory power thus returns to the same level it can sustain without the BW pathway, i.e., in the FW-isolated configuration (Fig 3c, note the difference in scales between colormaps). These two observations (i.e., the divergence in temporal frequency and the loss of oscillatory power) suggest at least some degree of independence between the two pathways, even when they are interconnected in the full model. Importantly, however, the intrinsic frequency of the isolated FW-pathway is consistently lower than that observed in the full model ($SG_x$ peak frequency difference at $\Delta T = 12$ms: 2.69 Hz). This is true for all tested delays, i.e., both in the non-diverging and the diverging regime. This offset in frequency directly shows that the intrinsic rhythmicity of the $IG_{IB}$
pacemakers largely dominates the frequency of the FW pathway. It is important to note that all of the results shown in Fig 3c and 3d represent oscillatory activity during the FW wave state. Thus, even though the temporal frequency of the network follows the infragranular pacemaker (i.e., the BW pathway) and is more or less stable between stimulus ON and OFF states, the spatial propagation of the activity (i.e., the relative phase between areas) reverses when the FW pathway is activated. We confirmed that this reversal of wave direction remains stable also at the scalp-level across the range of delays in the non-diverging regime (up to 20 ms) (S1 Fig).

In sum, we showed here that the two pathways have independent mechanisms to generate oscillations. However, in the fully connected model and for a fixed delay of 12 ms, the infragranular pacemaker largely determines the temporal frequency of the network. In contrast, the involvement of the FW pathway only determines the directionality of the phase gradient.

## The pulvinar pathway as a mechanism to bias traveling wave direction in the cortex

After demonstrating that the proposed cortical model can explain the dynamics of alpha traveling waves at the EEG level, we next investigate the role of the cortico-thalamic network, specifically the pulvinar, to explore how its connectivity affects wave dynamics in the cortex. The cortico-pulvinar circuitry (including the thalamic reticular nucleus, TRN) comprises a complex set of connections involved in seemingly diverse functional roles (for overviews, see, e.g., [21,43–46]). A detailed analysis of how each part of that circuitry could modulate wave dynamics is beyond the scope of this study. Here, we focus on the main Cx-Pul-Cx pathway that connects any two areas connected directly in the cortex (i.e., replication principle [21]). Furthermore, within this pathway, we focus on the driving feedforward projections extending from deep layers to the pulvinar and then to Layer 4 of the higher cortical area [47] (while noting that both feedback, as well as non-driving, modulator projections also exist between some cortical regions and the pulvinar; see Discussion). Importantly, we also do not consider the cortico-thalamic system as an additional, independent generator of alpha oscillations, as suggested by some experimental findings [25,26] and modeling [48]. As before, we aim to keep our model as simple and general as possible by reducing connectivity to what is likely shared between most hierarchically connected visual areas.

Fig 4a shows the schematic of our model extended by the trans-pulvinar pathway (in orange). In the logic of the FW- and BW-pathways described above, the connection through the pulvinar may be seen as running parallel to the FW-pathway in that it provides a second feedforward connection into Layer 4 of each area. The critical distinctions from the FW-pathway are 1) that it is unidirectional, i.e., it does not support the generation of its own oscillation through a feedback loop, and 2) that it is active in the absence of external input by the intrinsic activity propagating from the infragranular nodes. Thus, the pulvinar in this configuration provides an excitatory shortcut feeding spontaneous $IG_{IB}$ activity into the FW pathway.

Fig 4b shows how the addition of this pathway affects the scalp-level wave dynamics in our model: for a high pulvino-cortical weight, the response curve for the FW wave state (left panel, same as Fig 1c) is shifted to a higher baseline probability at the lower bound. The response curve for BW waves (right panel) is correspondingly shifted in the opposite direction (while the response threshold, i.e., the sigmoid inflection point, remains similar in both cases). In other words, spontaneous network activity (without any input current) is more likely to generate FW waves when the pulvinar pathway is engaged. This pattern is also visible when continuously mapping the response curve to the weight of the pulvino-cortical connection (Fig 4c). The dashed line in the upper panel marks the estimated threshold input current (50% point of the sigmoid), while the lower panel shows the log-ratio between FW and BW waves at baseline (no stimulation).

Since they are based on scalp-level classification, the response curves in Fig 4b and 4c do not allow conclusions about the activity in the cortex that underlies the FW-biased state. Indeed, the baseline probabilities for both states (y-axis intercepts of orange curves in Fig 4b) sum to a greatly reduced percentage of all fits for the higher pulvinar weight, indicating more ambiguous/unclassified fits. To elucidate this, we compared three states (marked in Fig 4b and 4c): 1) the stable BW wave state (without external input) in the cortical-only model, i.e., with the pulvinar disconnected, 2) the stable
<figure><img src="image_10.png" /><figcaption>**Fig 4.** Modulation of FW waves by the pulvinar. A: Model architecture with the added pulvinar pathway. For every two successive areas, one pulvinar node relays a feed-forward connection from infragranular layers into Layer 4. All cortical connectivity remains the same as before. B: Response curves for FW (left) and BW (right) waves for the Cx-only (grey) and Cx+Pul (orange) configurations. C: Full map of FW state probability as a function of pulvinar engagement (w<sub>Pul</sub> → L<sub>4X</sub>). The dashed line indicates the estimated response threshold. The bottom panel shows the log-ratio of FW and BW probabilities at I<sub>Stim</sub> = 0 nA. D: Spatiotemporal maps of cortical activity associated with each wave state labeled 1-3 in B & C. Each sub-panel shows SG<sub>X</sub>/IG<sub>IB</sub> activity across cortical areas and time, averaged over epochs centered on the peak activity at Cx<sub>2</sub>-IG<sub>IB</sub>. For spontaneous FW waves, results for a 5-area version of the model are shown for comparison. The maps show that the spontaneous FW waves associated with pulvinar engagement originate in Cx<sub>2</sub>, i.e., are shifted downstream compared to evoked FW waves. E: Distributions of mean oscillatory power across the network for 3-area (left) and 5-area versions, comparing the wave states described above. Error bars represent +/-1 SD between model runs. Colored arrows show the direction of propagation.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1013294.g004

stimulus-evoked FW state in the same configuration (cortex-only), and 3) the state without external input in which sponta-
neous FW waves emerge as a result of the added pulvinar pathway (w<sub>Pul</sub> → L<sub>4X</sub> = 2). For each state, we consider the spatial
activity maps (obtained by concatenating firing rates from nodes of the same type across the cortical hierarchy) in SG<sub>X</sub>
and IG<sub>IB</sub> nodes after averaging across epochs centered on the peaks in the rhythmic IG<sub>IB</sub> activity at Cx<sub>2</sub> (shown in Fig 4d).
We chose this reference for its position in the center of the network (Cx<sub>2</sub> for N<sub>Cx</sub> = 3) and because the IG<sub>IB</sub> node is at the
intersection between FW and BW pathways. The maps for BW (1) and evoked-FW states (2) show the expected pattern
of propagation with Cx<sub>3</sub> and Cx<sub>1</sub> activity leading in both laminar compartments, respectively (Fig 4d). In contrast, during
the spontaneous-FW state (3), activity at Cx₂ leads the other areas, while the supragranular nodes in Cx₁ remain inactive. This result illustrates how the spontaneous FW waves are elicited in the Cx+Pul configuration: spontaneous IGᵢB activity in Cx₁ is relayed to Cx₂ via the pulvinar, activating the FW pathway in the same way as the external input does in Cx₁. This results in a spontaneous FW wave spatially originating at Cx₂. The infragranular activity in Cx₁ disrupts the spatial phase gradient, leading to an unstable classification at the scalp level.

The activity pattern of spontaneous FW waves becomes more evident when considering a network with more than three areas. The two maps on the right in Fig 4d show the case for N<sub>Cx</sub> = 5 as an example. Here, the consistency in the FW wave is more clearly visible, while the leading activity remains at Cx₂.

Lastly, the two FW wave states (evoked/spontaneous) also differ in their spatial distribution of activity levels and oscil-
latory power (Fig 4e). As a result of cascading effects in the FW- and BW-pathways, overall activity levels in either state
naturally increase with the direction of propagation (i.e., more oscillatory power at regions trailing in phase; Fig 4d; BW:
IG<sub>IB</sub>, FW: SG<sub>x</sub>). Since physiological firing rates are subject to strong normalization from various sources (and were not the
target of our modeling), differences in the baseline firing rate can be disregarded. Instead, we tentatively interpret oscilla-
tory power as a measure for the spatial distribution of wave amplitude (shown in Fig 4e) (note that we deliberately refrain
from modeling power topographies at the scalp level as these depend significantly on assumptions about the size of
active populations in each area). The power distribution in the evoked FW state peaks at the center of the network for both
versions of the model (3-/5-stage). For spontaneous FW waves, this distribution is skewed towards downstream areas
due to the shift in wave origin from Cx₁ to Cx₂.

In summary, adding the pulvinar feedforward pathway in our model leads to the emergence of FW wave patterns in the baseline state while leaving the response to external stimulation unaltered. These spontaneous FW waves differ from those evoked by stimulation in that their spatial origin is shifted to the first area downstream of pulvinar feedforward activation. In our simulations, where the pulvino-cortical weights w<sub>Pul</sub> →<sub>L4X</sub> are identical across regions, the leading area is the second lowest (i.e., Cx<sub>2</sub>), but a differential modulation of these weights would allow for a flexible control of the spatial distribution of the wave.

Lateralization of cortical wave states

In the previous section, we demonstrated that our network's wave state (FW/BW) can be biased through the pulvinar
pathway. Interestingly, a possible cognitive function associated with this modulation type is the control of attentional allo-
cation. Empirical observations have suggested that traveling wave direction at the scalp level is lateralized as a function
of visual attention [13], similar to alpha-band oscillations [49-51]. As we will demonstrate in the following, our model can
reproduce such lateralization of wave states consistent with hemifield-directed attention through modulation of pulvinar
engagement. Specifically, we adapted our simulations and wave fitting procedure to differentiate between hemispheres
(Fig 5a). The novel network consists of two independent thalamo-cortical streams, each following the same architecture
and connectivity as above. In this configuration, we did not include direct cortical or subcortical interconnections between
the two hemispheres. Each hemisphere receives its input current into the first area, simulating left and right visual field
stimulation. After projecting all cortical sources together to the EEG, we estimate wave states in two separate regions of
interest (ROIs) lateralized to either hemisphere in the sensor space (Fig 5a). We first tested the effect of lateralization in
this model by passing the input current to only one hemisphere. Fig 5b shows the probability time courses for FW state
classification in both ROIs across a 1 sec left-hemisphere DC pulse stimulation. The response clearly shows the expected
lateralization to the stimulated side. Notably, the opposite ROI does not show a residual response (an intuitive expectation
one may have from the spatial mixing of signals in the EEG projection) because the right hemisphere remains in a sta-
ble BW state during stimulation. This demonstrates, as a pre-requisite for the modeling of hemispheric modulation, that
simultaneous waves traveling in opposite directions between the two hemispheres can be picked up reliably from a single
simulated EEG signal. However, the effect of signal mixing is evident in an overall reduction of classification consistency
a

<figure><img src="image_11.png" /></figure>

b

<figure><img src="image_12.png" /></figure>

c

<figure><img src="image_13.png" /><figcaption>**Fig 5.** Lateralization of cortical waves and hemispheric modulation by the pulvinar. A: Illustration of the multi-scale model with a separate cortical stream per hemisphere. Waves states are estimated from separate ROIs based on a single EEG projection. B: FW state probabilities for the two ROIs in response to a 1 sec left hemisphere stimulation. The responses clearly show the expected lateralization of wave states. The topographical plot shows the root-mean-square (RMS) within the first 100 ms after stimulus onset to illustrate the lateralization in the time-domain. C: Response curves for each ROI with equal (left panel) and left-lateralized (right panel) pulvinar engagement ($w_{Pul→L4X}$). The response shift shows that the effect of hemispheric bias from pulvinar modulation can be measured at the scalp level.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1013294.g005

(here, $I_{Stim}$ = 1.2 nA, compare with Fig 1c). For reference, the topographical plot on the right shows the lateralization of the same signal in the time-domain (root-mean-square (RMS) within the first 100 ms after stimulus onset).

Lastly, after establishing that our model can generate lateralized wave states (simultaneous opposing waves), we also investigated whether it can reproduce hemispheric modulation of waves traveling in the same direction. For this, we mapped responses in the two ROIs separately while passing the same shared input current to both hemispheres. Importantly, we varied the weight from pulvinar to cortex selectively for the left hemisphere, while keeping the weight for the right fixed at a low initial value (Fig 5c). The results reveal that the response curve for the left ROI is selectively modulated, confirming the shift in baseline of the FW state probability with increased pulvino-cortical weight. Additionally, due to the lower SNR in the responses (as compared to the non-lateralized version), the modulation persists even in the supra-threshold state (high input current), in which a FW wave is evoked in both hemispheres. This simulates the evaluation of EEG recordings better than the more ideal response curves in Fig 4b. Overall, our results show that a lateralized modulation of pulvino-cortical weights can lead to a hemispheric bias in FW wave consistency, similar to those obtained for selective visual attention. As shown experimentally [13], this lateralization can be reliably detected at the scalp-level, despite mixing signals from both hemispheres in the EEG.

Discussion

We constructed a multiscale network model that reproduces the pattern of alpha-band travelling waves at the level of the scalp, building on physiologically plausible neural dynamics. The modulation of FW/BW propagation along the cortical hierarchy by stimulus input qualitatively matches that observed in human EEG data. We have demonstrated that the generation of FW and BW waves in our model relies on distinct cortical sub-components (FW- and BW-pathway, respectively).
Lastly, we have further shown that the addition of a redundant feed-forward pathway through the pulvinar biases the network's dynamics towards the FW-wave state in the resting (stimulus OFF) state. This type of control can be used to simulate hemispheric modulation of traveling wave direction, as measured via EEG. All in all, starting from a few established anatomical principles of cortical organization [21,42,52], we describe an architecture that accurately replicates the pattern of traveling waves observed at the scalp level. Importantly, the aim of our work was not to explore the complete possible parameter space of the model but to propose a relatively simple framework of cortical and laminar organization with reliable explanatory power.

## From source to sensor space: Linking scalp-level traveling waves to laminar circuits

A number of EEG studies in recent years have found evidence that traveling waves at the scalp-level are linked to cognitive and sensory functions [12–14,53–55]. Similarly, traveling waves at different scales have been observed in a variety of paradigms, both intracortically and in recordings from the cortical surface [9,56–58]. It remains unclear how the findings at these two levels (scalp and cortical) relate to each other. The source-level correlates of the scalp waves are ambiguous. This is in part due to the well-established forward problem (for which standardized solutions exist [59]). However, more importantly, phase gradients at the scalp can equally result from intracortical traveling waves and from sequential activation of phase-lagged stationary sources [15,60]. Lastly, establishing a direct link with (simultaneously recorded) cortical activity is methodologically challenging.

Our multiscale model constitutes a first effort towards bridging this gap. It describes a biologically plausible laminar circuit to explain the traveling wave dynamics observed at the scalp-level using EEG recordings. It also generates predictions about the cortical activity expected during FW and BW wave states, which could be tested experimentally in future studies combining multi-scale measurements.

A central prediction made by our model is that the scalp-level propagation of waves corresponds to a propagation along the hierarchy, i.e., a sequential activation of cortical areas that each oscillate as stationary sources with phase lags between areas. This is consistent with the conclusion drawn by other authors that have identified long-range cortico-cortical fibers as the likely substrate of waves measured by the EEG based on conduction velocities [11,61] (see also the dependence of wave stability on interareal delays in [19]). Notably, this is not equivalent to a continuous wave in the cortex, as our model does not make predictions about cortical tissue between two areas in the hierarchy, or non-activated regions within an area. Yet, a continuous propagation seems more likely, given experimental evidence of continuous traveling waves at meso- and macro-scales [9,57]. However, several different mechanisms of propagation would likely interact to generate a continuous wave from the successive activation of areas predicted by our model (see [61] for an overview of the neural mechanisms driving propagation at different scales). Additionally, a recent study showed that scalp-level waves elicited during visual stimulation can sufficiently be explained by just two cortical sources [15]. Our model is in line with this in principle (given that it generalizes to all cases with more than one area), however, our hypothesis would be that the number of sources needed to explain the wave depend on the particular functional pathway activated. This could readily be tested in future EEG/MEG experiments.

Another important prediction from our model is that FW and BW wave states are associated with different laminar distributions of alpha activity. Specifically, alpha power in superficial layers should be stronger during the FW state (cf. Fig 3a), whereas rhythmic activity in the deep layers is continuous (and only switches in directionality of the phase offset between states). While a dependence on global wave state remains to be tested experimentally, there is already a large body of evidence showing that different neural rhythms have distinct laminar distributions. Generally speaking, low-frequency activity (including alpha) seems to be stronger in infragranular layers, while higher frequencies are more prominent in supragranular layers [33,34,62]. Pacemaker neurons in Layer 5 have been identified as a potential generator of the infragranular alpha [35,36], and the main generator for alpha in our model (IG<sub>Ib</sub>) reflects these findings. However, other studies report evidence for separate generators in both deep and superficial layers [63–65]. In line with this, the network
architecture in our model allows for the generation of a second alpha rhythm through interareal feedback to supragranular layers. It will be crucial for our understanding of cortical traveling waves to investigate how potential different alpha generators interact across different cortical laminae.

While using regular pacemaker neurons to generate alpha in our model allowed us to investigate wave propagation in isolation, it remains a somewhat artificial substitute. Several biophysical models have been put forward to explain the generation of the alpha rhythm (or rhythms). Among these, some identify possible sources in local cortical circuits [66–68], or cortico-cortical connectivity [69], while others highlight the contribution of the thalamus and/or thalamocortical loops [38,70,71]. It is possible that the dynamics of traveling wave direction as described by our model would change if a more realistic generator of alpha were used, and more work would be needed to integrate ours with the existing models. A key question to answer in this regard will be to expand the role of the thalamus (including the pulvinar) as a possible driving source of cortical alpha (discussed further below).

As noted earlier, a key limitation to EEG investigations of traveling waves is that any source-level interpretation of the data relies on inference. We attempted to include some of the biases related to this with our multiscale approach by evaluating wave direction directly at the scalp-level. However, given our focus on wave direction, a realistic simulation of the EEG was not within the scope of our model. Indeed, this limitation is reflected in our scalp-level estimates of wave propagation velocities (S2 Fig) (which, unlike wave direction, are strongly influenced by the number of areas, exact positioning of sources, the relative scaling of their activity, etc.). Our model should, therefore, not be taken as a proposed solution to the inverse problem for scalp-level waves, - instead, it provides a generalized architecture that can be used to model exact pathways using more realistic EEG projections in future studies. These investigations may also use a more flexible measure to quantify waves in the EEG (e.g., non ROI based, accounting for spherical distortion) to explore how their spatial topography is influenced by the topology of the underlying cortical network.

Our analyses show that wave direction at the scalp is rather stable across a broad range of interareal delays. For generalizability, we chose to keep this parameter fixed in our main simulations. However, more anatomically realistic delays, as well as variable neural time-scales between areas will be important to include in future models of specific pathways. Interestingly, our findings suggest that the directionality of scalp-level phase gradients may behave similarly across different configurations of activated areas, provided their hierarchical organization is aligned with the anterior-to-posterior axis. For example, a sequential activation of two areas may result in similar dynamics of FW/BW waves whether the sources are located in occipital and parietal, or in parietal and frontal regions, even though the interareal delays differ greatly between the two configurations. While this speaks for the use of (binary) wave direction as a more global parameter of directed communication, it implies that other wave parameters (spatial frequency, topographical distribution) may be more informative to distinguish between different networks at the source level.

## Alpha traveling waves as correlates of forward- and backward-directed communication in the cortex

Our model distinguishes largely between two global states, with waves propagating either FW or BW along the cortical hierarchy. Here, these two states simply reflect the directionality of the phase gradient between areas, but there is broad evidence to support a relationship between the direction of propagation and cortical communication. For example, in EEG recordings, alpha-band FW waves have been linked to active visual processing and attention, whereas alpha-band BW waves are more prominent during rest and attentional suppression [13,14].

In contrast to this, in recent years intracortical recordings have led to the emerging view that communication directions are separated among different frequency 'channels', with higher frequencies (gamma) carrying feedforward and lower frequencies (alpha/beta) carrying feedback signals [72–74]. The two channels interact in that gamma activity (in supragranular layers) is modulated by the phase of (infragranular) alpha [75], in line with the concept of cross-frequency coupling [76].

Our model, supported by experimental data [13, 14, 19], challenges the notion that alpha rhythms represent purely feedback signals. Instead, it suggests that alpha provides a constant global rhythm whose phase gradient can switch
directionality to support communication in either direction. Importantly, the high-frequency (gamma) channel could be
easily incorporated into our model, e.g., by a separate node in supragranular layers generating gamma through a
pyramidal-interneuron gamma (PING) mechanism. The supragranular activity in the current model would translate to an
amplitude modulation of the gamma envelope in that node. Similar accounts of trans-laminar cross-frequency cou-
pling have been put forward in computational models [37,77] and are supported by physiological observations [75,78].
Implemented in our model architecture, the alpha phase modulating the gamma activity would be generated by the
inter-areal feedback connection. Additionally, a second, direct projection to the supragranular layers from deep layers
of the same area is also likely [79,80]. While the cross-frequency interaction was not the target of our modeling, we
consider this integrated model (i.e., one with separate alpha- and gamma-frequency channels) here because the distinc-
tion between alpha-phase-modulated- and intrinsic alpha activity may explain the seemingly contradicting roles of alpha
between EEG and laminar recordings. In other words, the scalp-level alpha during the FW wave state may represent
rather the phase-modulated gamma activity in supragranular layers, while BW alpha would represent infragranular
intrinsic alpha.

Previous studies provide some support for this more flexible role for alpha oscillations (as reflecting cortical com-
munication more generally, as opposed to a feedback-specific channel). Using ECoG in humans, Bahramisharif et al.
[16] found that the amplitude of gamma activity is modulated by the local phase of large-scale alpha traveling waves,
effectively resulting in the spatial propagation of gamma bursts. Chapeton et al. [81], without characterizing traveling
waves specifically, showed experimentally that the communication between regions is optimal when their alpha activity
is coherent and phases are aligned to match the conduction delays between them. A similar framework of alpha-gated
communication has been proposed by Bonnefond et al. [82], and demonstrated in a biologically plausible network model
[24]. Lastly, Halgren et al. [18], using ECoG and sEEG in humans, found traveling alpha waves that propagate towards
the occipital pole during rest, consistent with our prediction that the resting alpha rhythm is backward-directed. They also
report evidence of a reversal of waves to the FW direction after opening of the eyes in a macaque (but only in some of
their human subjects, see their Supplemental Figure S4 [18]). Future studies could establish the dynamics of traveling
wave direction during similar transitions from rest to active states.

Interestingly, Halgren et al. [18] localize the laminar source of alpha selectively in the superficial layers, contrary to
previous studies and our model predictions that assume either distributed or infragranular sources. It is possible that the
discrepancy in the literature partly reflect the use of different analytical measures (current source density (CSD) in Halgren
et al. [18] vs. local-field potential (LFP) in most studies). However, it seems possible that waves that are generated by
infragranular pacemakers down-stream in the visual hierarchy propagate further via short-range connections in superfi-
cial layers. This would be an indication that different laminar distributions are to be expected at the presumed source of a
wave and at remote cortical regions through which it propagates.

A predictive coding interpretation of the FW and BW states

The interlaminar and interareal connectivity in our model is inspired by the hierarchical predictive coding model developed by Alamia & VanRullen [19]. It also follows the microcircuitry proposed by existing laminar models of predictive coding [42,52,83]. Thus, even though our model itself is dynamical and does not include a feature space, it supports the computation and passing of prediction errors between areas. Specifically, in line with previous studies, it predicts that this error computation occurs in superficial layers. An explicit representation of this prediction error could be included in the model by adding a local relay of supragranular activity to a second excitatory population, from which the feedforward projection then extends to the higher area. A supragranular computation of prediction errors is supported by recent electrophysiological studies who found evidence for a neural comparison between top-down and bottom-up signals in Layers 2/3 [84–87]. Complementary to this, Thomas et al. [88] were able to show with 7T fMRI in humans that stimuli with large prediction errors were represented selectively in superficial layers. Interestingly, model implementations of the local circuitry involved
in predictive error-computation exhibit naturally emerging oscillatory dynamics [89], supporting a direct link between neural oscillations and predictive coding.

Recent work has furthermore suggested that different cortical rhythms (alpha/beta and gamma, see above) may sup-
port different functions in predictive processing. Specifically, the ‘predictive routing’ model [90] poses that the processing
of expected stimuli (at high frequencies, in superficial layers) is inhibited by increased infragranular alpha activity, mod-
ulated through stronger alpha-feedback from higher areas. Follow-up work shows that predictive routing could emerge
naturally from biophysically plausible learning constraints [91]. As discussed above in the context of laminar cross-
frequency coupling, our model supports this type of frequency-specific communication, if high-frequency (gamma) activity
is introduced at superficial layers and its amplitude is modulated by infragranular alpha. This framework generally endorses
the interpretation of BW waves as representing feedback signaling of predictions, and FW waves as representing feedfor-
ward signaling of sensory input. However, the predictive routing model suggests that alpha-rhythmicity should decrease
in the absence of expectations (while superficial activity shifts to a gamma-/transient state). It is unclear whether stable
FW waves at alpha would still be expected in this state, and more work will be needed to explore the functional interplay
between transient and rhythmic network states.

Similarly, it remains unclear whether the functional role of alpha traveling waves is global or local. In particular, consid-
ering the retinotopic representation at early stages of the visual hierarchy (V1/V2), it needs to be shown to what extent the
reversal of wave direction is spatially selective, to support the feedback of spatially informative predictions to these areas.
There is already some evidence from the EEG that different wave states can co-occur, triggered by stimuli at separate
positions in the visual field [10,13], and our own simulations with two separated hemispheric streams (Fig 5) also demon-
strate this. However, the low spatial resolution of the EEG limits the conclusions to be drawn. Future studies could expand
our dynamical model to include a feature space, in order to investigate how feature-based predictions (such as about the
position of a stimulus) interact with traveling wave dynamics.

Interestingly, a recent computational model showed that meso-scale waves, traveling horizontally across retinotopic
space within an area, can also carry predictions about future visual input [92]. This corroborates experimental evidence
that similar meso-scale waves occur spontaneously in the cortex and modulate local spiking activity as well as perception
[57]. An interesting venue for future studies will be to explore how the two types of waves (meso-scale within areas vs.
feedforward/feedback between areas) interact, both dynamically and in the potential encoding of predictions and predic-
tion errors.

The role of the pulvinar

Our simulation results showed that the addition of second pathway through the pulvinar biases our network towards the
FW wave state. A detailed analysis revealed that the spontaneous FW waves at rest propagated from the second cortical
area (the first receiving input from the pulvinar), in contrast to stimulus-evoked FW waves, which propagated from the
lowest area.

These results are in line with a number of studies (experimental and modeling) showing that the pulvinar plays a key
role in the control of communication between cortical areas [22–26]. In particular, the modeling effort by Jaramillo et al.
[23] links the pulvinar to diverse cognitive functions as a flexible gating mechanism in the flow of information between
cortical areas. Together with our results, this opens an intriguing connection between the modulation of traveling wave
direction and cognitive functions. There is initial evidence for a similar link in the hemispheric modulation of wave direction
with spatial attention [13]. The decision-making paradigm modeled in Jaramillo et al. [23] may provide another venue in
this direction for future EEG studies on wave dynamics.

Furthermore, considering the role of the pulvinar in generating spontaneous FW waves may offer a potential expla-
nation for the between-subjects variability we observed in experimental data at rest in the ratio between FW and BW
waves. Unspecific attention or general arousal may modulate pulvinar engagement in this state and lead to a higher or
lower probability of FW wave patterns at the scalp-level. Conversely, pulvinar input to the cortex would be expected to be down-modulated during the passive resting state and periods of attentional suppression, in order to prevent spontaneous FW waves to be triggered. The thalamic reticular nucleus (TRN) may serve this kind of gating function [93,94].

Interestingly, our model predicts that the spatial distribution of spontaneous FW waves should be more towards frontal regions, as compared to stimulus-evoked waves. With a suitable resting state paradigm, this could be tested in future studies.

It is important to note that our investigation of the pulvinar pathway leaves out several known properties of the thalamo-cortical connectivity. First, the pulvinar nodes in our model only passively relay activity, while physiological thalamic relay cells form a highly specialized circuit with the thalamic reticular nucleus that together exhibits complex non-linear response properties [95]. Moreover, cortico-thalamic projections also comprise two distinct subtypes (drivers and modulators), and the interaction between these appears to play a role in the generation of alpha oscillations in the pulvinar itself [45,48,96,97]. Lastly, our model does not include the backward-directed connection from one cortical area through the pulvinar to a lower area (see [21]). This pathway is unlikely to assume a modulating role for BW waves equivalent to the one predicted for FW waves in our model, given that pulvino-cortical projections terminate largely in superficial layers [46,98]. However, pulvinar input has been shown to play an important role in refining representations in V1 and gating its output [99]. This makes the pulvino-cortical feedback pathway a likely candidate for carrying prediction signals in the context of predictive coding [100,101]. It will be important to consider this pathway in future implementations that expand the current dynamical framework to a feature space.

## Conclusion

In this work we presented a multiscale mean-field network model to explain and simulate the dynamics of alpha-band traveling waves measured with the EEG. We proposed two distinct inter-laminar pathways for the propagation of forward and backward alpha waves, and we modeled the role of the pulvinar in modulating the cortical pathways. Importantly, we ground the proposed architecture and the results of our simulations in the theoretical framework of predictive coding, in line with previous work. All in all, our model provides a first theoretical framework connecting these scalp-level waves to cortical circuits and could provide a base architecture for future research in several directions.

## Materials and methods

### Ethics statement

Our study includes a re-analysis of previously published data [14]. All participants in that study gave written informed consent, and the experimental protocol was approved by the committee 'Comité de protection des Personnes Sud Méditerranée 1' (ethics approval number N°2016-A01937-44).

### Model dynamics

Our model is based largely on the mean-field model introduced by Jaramillo et al. [23]. We used the same dynamics to model all cortical and thalamic nodes and their connectivity, except for the infragranular pacemaker (IGIB) described by average single-neuron dynamics. In the following, we first describe the regular mean-field dynamics and then extend it to the special case of the IGIB. Exact values (or ranges wherever parameter spaces were mapped) for all model parameters are listed below in Table 1.

Each node i is described by its firing rate $r_i$ varying over time as a fixed function of the input current $I_i$ (FI-curve):

$$r_i = \begin{cases} \frac{\lambda_i l_i - \beta_i}{1 - e^{-\beta_i}(\lambda_i l_i - \beta_i)} & \text{for } l_i \neq \frac{\beta_i}{\lambda_i} \\ \frac{1}{\theta_i} & \text{otherwise} \end{cases} \quad (1)$$
<table><caption>Table 1. List of parameter values used in the model simulations.</caption><thead><tr><th colspan="2">Mean-Field Dynamics</th></tr></thead><tbody><tr><td>λ</td><td>270 Hz/nA</td></tr><tr><td>β</td><td>108 Hz</td></tr><tr><td>θ</td><td>0.154 s</td></tr><tr><td>I<sub>base</sub></td><td>0.33 nA</td></tr><tr><td>σε</td><td>0.025 nA</td></tr></tbody><thead><tr><th colspan="2">Mean-Field Synaptic Parameters:</th></tr></thead><tbody><tr><td rowspan="5">[τ<sub>syn</sub>: γ<sub>syn</sub>]</td><td>Fast Excitatory</td><td>[20 ms; 0.8]</td></tr><tr><td>Slow Excitatory</td><td>[120 ms; 0.5]</td></tr><tr><td>Fast Inhibitory</td><td>[3 ms; 0.8]</td></tr><tr><td>Slow Inhibitory</td><td>[120 ms; 0.1]</td></tr><tr><td>Infragranular Self-Excitatory</td><td>[1 ms; 0.8]</td></tr></tbody><thead><tr><th colspan="2">Spiking Neuron Dynamics:</th></tr></thead><tbody><tr><td>a</td><td>0.0067</td></tr><tr><td>b</td><td>0.2</td></tr><tr><td>c</td><td>-50 mV</td></tr><tr><td>d</td><td>2</td></tr><tr><td>[με; σε]</td><td>[6; 4]</td></tr></tbody><thead><tr><th colspan="2">Inter-Areal Delay:</th></tr></thead><tbody><tr><td>ΔT<sub>Cx→Cx</sub></td><td>12 ms</td></tr><tr><td>ΔT<sub>Cx Pull→Cx</sub></td><td>6 ms</td></tr></tbody><thead><tr><th colspan="2">Local Connections (Intra-areal):</th></tr></thead><tbody><tr><td>L4<sub>x</sub> → SG<sub>x</sub></td><td>w=1</td></tr><tr><td>L4<sub>x</sub> → L4<sub>IN</sub></td><td>Fast Excitatory</td></tr><tr><td>SG<sub>x</sub> → IG<sub>IB</sub></td><td>w=0.7</td></tr><tr><td>IG<sub>IB</sub> → IG<sub>IB</sub></td><td>Infragranular Self-Excitatory</td></tr><tr><td>L4<sub>IN</sub> → L4<sub>x</sub></td><td>w=-1.2</td></tr><tr><td>SG<sub>IN</sub> → SG<sub>x</sub></td><td>Slow Inhibitory</td></tr><tr><td>IG<sub>IB</sub> → IG<sub>IB</sub></td><td>w=-2</td></tr><tr><td>SG<sub>x</sub> → L4<sub>x</sub></td><td>Fast Inhibitory</td></tr></tbody><thead><tr><th colspan="2">Feedforward Connections (cortico-cortical):</th></tr></thead><tbody><tr><td>IG<sub>IB</sub> → SG<sub>IN</sub></td><td>w=2</td></tr><tr><td>IG<sub>IB</sub> → IG<sub>IB</sub></td><td>Fast Excitatory</td></tr></tbody><thead><tr><th colspan="2">Feedback Connections (cortico-cortical):</th></tr></thead><tbody><tr><td>IG<sub>IB</sub> → SG<sub>IN</sub></td><td>w=1.5</td></tr><tr><td>IG<sub>IB</sub> → IG<sub>IB</sub></td><td>Fast Excitatory</td></tr><tr><td>IG<sub>IB</sub> → IG<sub>IB</sub></td><td>w=1</td></tr><tr><td>IG<sub>IB</sub> → IG<sub>IB</sub></td><td>Slow Excitatory</td></tr></tbody><thead><tr><th colspan="2">Connections to/from Pulvinar:</th></tr></thead><tbody><tr><td>IG<sub>IB</sub> → Pul</td><td>w=1</td></tr><tr><td>Pul → L4<sub>x</sub></td><td>w∈ [0.22] (target parameter variation)</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1013294.t001
where r is in Hz, I is in nA, and λ, θ, and β are the FI-curve's slope (excitability) and offset (neural threshold) parameters. The input current $I_i$ to the node is given by:

$$I_i = \sum_j (w_{j \to i} s_{j \to i}) + I_{base} + \varepsilon_i \quad (2)$$

where the first part is the sum of the weighted synaptic input to node i (dynamics described below) $I_{base}$ is a constant base current, and $\varepsilon_i$ is time-varying noise specific to each node. The noise term is an Ornstein-Uhlenbeck process of the form:

$$\frac{d\varepsilon_i}{dt} = \frac{-\varepsilon_i}{\tau\varepsilon_i} + \sigma\varepsilon_i \eta(t) \sqrt{\tau\varepsilon_i} \quad \text{with} \quad \eta(t) \sim N(0, 1) \quad (3)$$

where $\tau\varepsilon$ and $\sigma\varepsilon$ are parameters for the time constant and the strength of the noise, respectively.

A synapse k transmitting activity of node i to another node is defined through a constant weight $w_k$, and a dynamic variable $s_k$ representing the gating of transmitter-mediated current flow, defined by the following dynamics:

$$\frac{d}{dt} s_k(t) = \frac{-s_k(t)}{\tau_k} + \gamma_k (1 - s_k(t)) r_i(t - \Delta T_k) \quad (4)$$

Where t is time, $\Delta T_k$ is the synaptic delay, $\tau_k$ is the time constant of the synapse, $\gamma \in [0, 1]$ is a saturation parameter and $r_i$ is the firing rate of node i as defined in (1).

## Infragranular pacemaker (IGIB)

To model the activity of intrinsically bursting populations in deep cortical layers, we substituted the mean-field response dynamics defined in (1) and (2) by a spiking-neuron model. For each node, we defined a homogeneous population of N=350 Izhikevich spiking neurons [102]. Every neuron is described by two dynamic variables - its membrane potential v and a recovery variable u, which evolve according to the following equations:

$$\frac{dv_n}{dt} = 0.04 v_n^2 + 5v_n + 140 - u_n + I_n \quad (5)$$

$$\frac{du_n}{dt} = a(bv_n - u_n) \quad (6)$$

Here, a and b are parameters determining the dynamic behavior of the membrane potential (see [102]). Once the potential passes a threshold, a spike is elicited and the values of v and u are reset:

$$\text{if } v \ge 30 \text{ mV}: \quad \begin{cases} v \leftarrow c \\ u \leftarrow u + d \end{cases} \quad (7)$$

where c and d are the reset parameters. The values for a, b, c, and d were chosen such that the neurons generate short bursts of spiking activity at rest in regularly paced intervals of approx. 100ms (see list of parameter values in Table 1).

The activity of the spiking-neuron nodes was integrated into the mean-field model by applying fixed conversions at input and output stages (i.e., substituting Eq. 1 and 2 above). For a given spiking node $i_{SN}$, its mean activity $r_{i_{SN}}$ is given by the instantaneous spike rate within the population:

$$r_{i_{SN}}(t) = \frac{1}{N} \sum_{n}^{N} S_n(t) \quad \text{with} \quad S_n = \begin{cases} 1 & \text{if } v_n \ge 30 \text{ mV} \\ 0 & \text{otherwise} \end{cases} \quad (8)$$
where N is the number of neurons in the population and $v_n$ is the membrane potential of neuron n as given by Eqs 5–7. As input, each neuron in the population receives the (scaled) input current to the node, plus a Gaussian noise term (independent for each neuron):

$$I_n = 4I_{[SN]} + \mu_{\epsilon} + \sigma_{\epsilon}\epsilon_n \quad \text{with} \quad \epsilon_n \sim N(0, 1) \qquad (9)$$

where $I_{[SN]}$ is the summed input current to the node as defined in (3), and $\mu_{\epsilon}$, $\sigma_{\epsilon}$ are the mean and standard deviation of the noise. Since the parameters a-d were identical between all neurons, the noise term in (9) is the only source of variability in the population. Note also that the neurons were not connected directly to each other. However, spiking nodes as a unit formed self-connections in the model (at the level of the mean-field input and output, i.e., Eqs 2 and 8). This was designed to retain synchronization of the intrinsic bursting behavior of the population as a whole.

All equations of the model were evaluated using Euler's method in time steps of 1 ms resolution, with the exception of the spiking neurons' membrane potential Eqs 5 and 6 for which a 0.5 ms step was used.

## Cortical model architecture and connectivity

The general architecture and connectivity of the network are depicted in Fig 1. Our model describes a single-stream hierarchy of $N_{Cx}$ (visual) cortical areas. Unless noted otherwise, all simulations use a version of the model with $N_{Cx}=3$. The general response behavior and dynamic reversal of traveling wave direction remained the same when we tested versions with up to 7 areas.

Each cortical area comprises six nodes across three laminar compartments: L4x and L4in in Layer 4, SGx and SGin in supragranular Layers 2/3, and IGib in infragranular Layers 5/6, where the IN subscript denotes inhibitory nodes. We focus our analysis mainly on SGx and IGib activity as the central nodes in the FW- and BW-pathways (see main text).

Our choice of model architecture was motivated in equal parts by the hierarchical predictive-coding model developed by Alamia & VanRullen [19] and physiological evidence of functional connectivity between laminar compartments. The predictive coding framework [20] postulates that feedforward projections mainly carry the difference (prediction error) between the information of the sending area and that sent as feedback from the receiving area (prediction). We describe how our architecture relates to the original model [19] in S1 File. As a reference for laminar connectivity (as well as existing neural implementations of predictive coding), we refer to previous work by Shipp et al. [52,83] and Bastos et al. [42]. Our architecture follows the dominant flow of information in the canonical microcircuit from Layer 4 to superficial and from there to deeper layers. We also assume that inter-areal feedforward connections project mainly to Layer 4 while feedback projections terminate in supra- and infragranular layers [28].

The input current in our model is applied to a separate input layer (ILst) comprising one pair of recurrently connected excitatory and inhibitory nodes (identical to Layer 4 of each cortical area). This stage was included for additional input normalization prior to $Cx_1$, and may be seen as equivalent to LGN.

An identical pair of nodes was included at the other end of the network (ILpr). The excitatory node in this pair was connected to the IGib node in the last cortical area ($Cx_N$) and received a continuous white noise input (varying between 0 and 0.3 nA). This was not intended to model a prior signal (as included in Alamia & VanRullen [19]), but it merely acts as a cortical top-down signal, ensuring stable activation of the BW-pathway.

Unless stated otherwise, cortical synaptic delays were fixed at 0 ms for local (intra-areal) and 12 ms for remote (inter-areal) connections. It should be noted here that the 'net' (e.g., as peak-to-peak) delay for activity sent from one node to another depends on the synaptic delay, synaptic time constant and neural integration times of the receiving and any intermediate nodes.

## Pulvinar pathway

The thalamocortical version of our model includes, in addition to the cortical modules, a single pulvinar node between any two consecutive cortical areas, modeling the functional feedforward connectivity of thalamic relay neurons. The
connection delay between cortex and Pulvinar was fixed at half the cortical inter-areal delay (i.e., 6 ms). That is, the summed delay of the connections Cx-Pul-Cx was the same as that of direct Cx-Cx connections. However, we confirmed that the main pattern of results for the pulvinar simulations were the same with double delays.

## EEG simulation

We used a simple forward model to obtain a scalp-level projection of our model's mean-field output. The primary goal of this was to compare wave dynamics at the laminar level in our model to existing EEG data, using a phase-based wave fitting procedure (see below).

The activity for each cortical node was first low-pass filtered using a 3rd order Butterworth filter (cutoff 20 Hz) and then down-sampled to 100 Hz from the 1kHz sampling rate of the simulations. For each simulated trial, a set of 5 independent noise sources (power spectrum slope $1/\sqrt{f}$) were generated and projected together with the laminar signals. The signal-to-noise ratio (with respect to the root mean square of both signals) for each noise source was drawn randomly from a uniform distribution between 0.4 and 1.6.

The forward model projections were performed using functionalities provided by the fieldtrip toolbox (ft_dipolesimulation; Oostenveld et al., 2011 [103]). Source dipoles were oriented radially and positioned according to their area, using the MNI-coordinates listed in Table 2 (based on the AAL atlas [104]).

For the single-stream version of the model (used everywhere except in the simulations of hemispheric lateralization), the cortical sources were positioned bilaterally in both hemispheres, scaled at half amplitude each. The noise sources were positioned randomly for each trial on a 5 x 5 mm grid covering all positions located inside the skull. Neither the pulvinar nor the stimulus input layer ($IL_{SI}$) were included in the forward projection.

## EEG data collection and analysis

The EEG data used in the comparisons with our model output have been reported in more detail previously [14]. In brief, observers (N=13, six females, mean age 25.57, range 21–31, two left-handed) were presented with a single, iso-luminant disk (7° diameter, centered at 7.5° above fixation) for a duration of 5 sec on each trial, followed by 5 sec of blank screen. The disk's luminance either followed a white-noise random sequence ("Dynamic" condition) or remained uniform ("Static"). We include only the "Static" condition in our comparison, which was composed of a total of 150 trials per subject. Throughout all trials, subjects performed an adaptive visual detection task requiring them to attend the stimulus covertly. The EEG was recorded using a 64-channel active BioSemi system at a sampling rate of 1024 Hz.

All pre-processing steps for the present analyses were identical to the previous study. Signals were recorded on a 64-channel BioSemi system at 1024 Hz and down-sampled offline to 160 Hz temporal resolution. Re-referencing was performed by subtracting the average activity and band-stop (50 Hz) and high-pass (> 1Hz) filters applied. Artifact detection and rejection were performed manually on the epoched data.

<table><caption>Table 2. Coordinates of dipole positions used for the EEG forward model. All values are MNI coordinates based on the AAL atlas [104].</caption><thead><tr><td>Model Area</td><td>Label</td><td>Coordinates</td></tr></thead><tbody><tr><td>Cx<sub>1</sub></td><td>Occipital</td><td>Left: [-8 -76 10]<br/>Right: [8 -76 10]</td></tr><tr><td>Cx<sub>2</sub></td><td>Parietal</td><td>Left: [-24 -61 58]<br/>Right: [24 -61 58]</td></tr><tr><td>Cx<sub>3</sub></td><td>Frontal</td><td>Left: [-5 48 30]<br/>Right: [5 48 30]</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1013294.t002
## Wave fitting procedure

To quantify the direction of wave propagation in both simulated and real EEG signals, we adapted the procedure introduced by Zhang et al. [9] for use with scalp-level data. This method relies solely on the analytic phase of the signal and is therefore well suited to analyze our simulated data (in contrast with the 2D-FFT method used in previous studies [13,19] that also factors in amplitude topography which was not the target of our simulations). Importantly, the goal of this procedure was not to quantify the wave at the source level (which is best described by the phase lag between areas). Instead, it intends to quantify the wave direction that would be observable at the scalp-level given the simulated cortical activity. We confirmed that our algorithm correctly estimates wave direction by simulating perfect sinusoidal FW and BW waves at the source level, using configurations with two and three phase-lagged oscillatory sources.

For a given frequency band of interest (7–13 Hz unless otherwise stated), the continuous phase of the signal at each electrode was obtained using a bandpass FIR filter and the analytic Hilbert transform. The resulting phase values were referenced to the mean phase across electrodes within the electrode ROI (described below) at each timepoint. The following iterative fit procedure was evaluated using the time-averaged relative phase within a sliding window 100 ms. For all analyses, we used a 2D projection of the electrodes' positions on the scalp (10–20 topography layout provided by the fieldtrip toolbox [103]). For a set of electrode positions x, y, the phase values $\hat{\theta}_{x,y}$ predicted by a perfect planar traveling wave are given by:

$$ \hat{\theta}_{x,y} = ax + by + \vartheta \quad (10) $$

where a, b are free parameters and $\theta$ is the reference- (mean) phase. The propagation direction of the wave is given by $\alpha = atan2(b, a)$ and its spatial frequency by $\xi = \sqrt{a^2 + b^2}$.

For each time-point, the model was evaluated iteratively within a given set of values for a and b. This parameter space was defined by 60 steps over the possible space of wave directions $\xi$ (360°) x 30 steps of increasing spatial frequency up to 1 cycle over the maximum distance covered by the electrode space. The best fit was determined by finding the maximum vector length of the summed residuals in circular space:

$$ \bar{r}_{\theta} = \sqrt{\left| \frac{1}{N_{x,y}} \sum_{x,y} \sin(\theta_{x,y} - \hat{\theta}_{x,y}) \right|^2 + \left| \frac{1}{N_{x,y}} \sum_{x,y} \cos(\theta_{x,y} - \hat{\theta}_{x,y}) \right|^2} \quad (11) $$

Further following Zhang et al. [9], the circular correlation $\rho$ between predicted ($\hat{\theta}$) and observed ($\theta$) phases was used as a statistical measure for the absolute goodness of fit:

$$ \rho = \frac{\frac{\sum_{x,y} \sin(\theta_{x,y}-\bar{\theta}) \sin(\hat{\theta}_{x,y}-\bar{\theta})}{\sqrt{\sum_{x,y} |\sin(\theta_{x,y}-\bar{\theta})|^2 \sum_{x,y} |\sin(\hat{\theta}_{x,y}-\bar{\theta})|^2}}}{\sqrt{\sum_{x,y} |\sin(\theta_{x,y}-\bar{\theta})|^2 |\sin(\hat{\theta}_{x,y}-\bar{\theta})|^2}} \quad (12) $$

with $\bar{\theta}$, $\hat{\theta}$ denoting the mean within the ROI of observed and predicted phases, respectively. To determine the significance of a given fit, we compared it to a chance-level distribution of $\rho$ values (using a one-sided threshold of 95%). This distribution was obtained once for each dataset by repeating the fitting procedure 10 times with randomly permuted electrode positions, at reduced temporal resolution, and collating the resulting fit values for all trials and time-points

Finally, each time-point was classified into one of three categories: as “Forward” (“Backward”) state, if the fit was statistically significant and the propagation direction was within 0.5 rad of the FW/BW axis (see below), and as a “Null” state otherwise. For non-lateralized stimulation, activation can be assumed to be symmetrical between hemispheres on average, therefore the reference axis is defined as the Fronto-Occipital axis. Spatial propagation in response to lateralized
stimulation can deviate from this axis [10,13]. For lateralized simulations, we therefore defined the reference as the median propagation direction after combining BW and FW directions into a single distribution through collapsing on the medio-lateral axis.

For all non-lateralized simulations, we chose a rectangular (in 2D-projected coordinates) ROI of 38 electrodes centered on the midline, spanning the occipito-frontal distance between Iz and Fz on the midline and between C3 and C4 on the mediolateral axis at its center. For lateralized simulations, we split our analysis into two separate ROIs, each centered on one hemisphere. To include a sufficient number of electrodes in the fit, these ROIs extended further laterally than the central ROI (spanning the distance from C5 to C6 between both ROIs at the center), excluding the midline, resulting in 20 electrodes per ROI. Our choices for both the central and hemispheric ROIs were based on previous observations [10, 13, 14, 19].

## Supporting information

**S1 File. Derivation of the FW-pathway architecture in our model from the hierarchical model in Alamia & VanRullen [19].**
(DOCX)

**S1 Fig. Effects of varying the inter-areal delay on power and fit wave direction in the simulated EEG.** Data correspond to the full model with rhythmic $IG_{IB}$ activity and $IG_{IB}$ - $IG_{IB}$ feedback intact (compare Fig 3c, left panels). In this configuration, the network enters a diverging regime for delays > approx. 20 ms (dashed line), in which FW- and BW-pathways differ in their temporal frequencies. A: Mean power spectra (averaged between electrodes on the mid-line) during stimulus On and Off periods as a function of inter-areal delay. The temporal frequency of the scalp-level activity remains constant independent of delay in the stimulus Off period (left panel), but decreases along with the $SG_x$ frequency with increasing delay (right panel, cf. Fig 3c and 3d). In this state, the activity beings to move out of the frequency band used for the wave fit (alpha, 7–13 Hz) at approx. 30 ms delay (shaded area), therefore this range of delays is excluded in the fitting procedure below. B: Probabilites of FW and BW wave states as a function of delay, during stimulus On and Off periods. Simulations were run with an intermediate stimulus amplitude ($I$=1.2 nA) between the reversal threshold and the upper asymptote to ensure dynamic range in the FW-state. The results show that the wave direction is stable within the non-diverging range of delays (up to 20 ms) while for higher delays, the FW state probability steadily decreases.
(TIFF)

**S2 Fig. Analysis of wave propagation velocities.** A: Spatial frequencies of planar wave fits, normalized to the phase difference predicted between Oz and Fz electrodes. Top panel shows boxplots of averages across participants in the real EEG data [14], bottom panel shows histograms of trial-averages obtained across 50 model simulations, using the same stimulus duration as the EEG experiment. The model EEG projection reproduces the empirical values for the BW state (stimulus OFF), but slightly overestimates the velocity of FW waves during stimulus ON periods. This effect is likely explained by our non-realistic EEG projection, which includes stronger transient activity at lower areas during the FW state (largely driven by sources in Layer 4; cf. the source level activity shown in Fig 1b, panel 4), and does not affect the direction of the fit (FW/BW). B: Source-level propagation velocities for models with different numbers of areas, normalized to the absolute phase difference covered by the entire hierarchy ($Cx_1$ to $Cx_N$). Each model was run 50 times using the same stimulus configuration as in A, values shown are means +/- 1SD across model runs. The comparison between $SG_x$ and $IG_{IB}$ phase differences (respectively, during stimulus ON and OFF periods) shows that FW wave propagation is highly spatially consistent, while the spontaneous BW propagation velocity is more variable (cf. again Fig 1b, panels 4 and 5). Note that the linear relationship between number of modeled areas and total phase difference only holds under the simplifying assumption that all areas have the same neural time-constants.
(TIFF)
## Acknowledgments

The authors would like to thank Xiaoqi Xu and Martin Vinck for helpful comments on the model, as well as Leslie Marie-Louise for administrative support. Views and opinions expressed are those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council (ERC). Neither the European Union nor the granting authority can be held responsible for them.

## Author contributions

**Conceptualization:** Jakob C. B. Schwenk, Andrea Alamia.

**Formal analysis:** Jakob C. B. Schwenk.

**Funding acquisition:** Andrea Alamia.

**Investigation:** Jakob C. B. Schwenk, Andrea Alamia.

**Methodology:** Jakob C. B. Schwenk.

**Software:** Jakob C. B. Schwenk.

**Supervision:** Andrea Alamia.

**Validation:** Jakob C. B. Schwenk, Andrea Alamia.

**Visualization:** Jakob C. B. Schwenk.

**Writing – original draft:** Jakob C. B. Schwenk.

**Writing – review & editing:** Jakob C. B. Schwenk, Andrea Alamia.

## References

1. Buzsaki G. Rhythms of the Brain. Oxford University Press. 2006.

2. Jensen O, Mazaheri A. Shaping functional architecture by oscillatory alpha activity: gating by inhibition. Front Hum Neurosci. 2010;4:186. https://doi.org/10.3389/fnhum.2010.00186 PMID: 21119777

3. Klimesch W. α-band oscillations, attention, and controlled access to stored information. Trends Cogn Sci. 2012;16(12):606–17. https://doi.org/10.1016/j.tics.2012.10.007 PMID: 23141428

4. Mathewson KE, Lleras A, Beck DM, Fabiani M, Ro T, Gratton G. Pulsed out of awareness: EEG alpha oscillations represent a pulsed-inhibition of ongoing cortical processing. Front Psychol. 2011;2:99. https://doi.org/10.3389/fpsyg.2011.00099 PMID: 21779257

5. Sadaghiani S, Kleinschmidt A. Brain Networks and α-Oscillations: Structural and Functional Foundations of Cognitive Control. Trends Cogn Sci. 2016;20(11):805–17. https://doi.org/10.1016/j.tics.2016.09.004 PMID: 27707588

6. Busch NA, Dubois J, VanRullen R. The phase of ongoing EEG oscillations predicts visual perception. J Neurosci. 2009;29(24):7869–76. https://doi.org/10.1523/JNEUROSCI.0113-09.2009 PMID: 19535598

7. Dugué L, Marque P, VanRullen R. The phase of ongoing oscillations mediates the causal relation between brain excitation and visual perception. J Neurosci. 2011;31(33):11889–93. https://doi.org/10.1523/JNEUROSCI.1161-11.2011 PMID: 21849549

8. VanRullen R. Perceptual Cycles. Trends Cogn Sci. 2016;20(10):723–35. https://doi.org/10.1016/j.tics.2016.07.006 PMID: 27567317

9. Zhang H, Watrous AJ, Patel A, Jacobs J. Theta and Alpha Oscillations Are Traveling Waves in the Human Neocortex. Neuron. 2018;98(6):1269-1281.e4. https://doi.org/10.1016/j.neuron.2018.05.019 PMID: 29887341

10. Lozano-Soldevilla D, VanRullen R. The Hidden Spatial Dimension of Alpha: 10-Hz Perceptual Echoes Propagate as Periodic Traveling Waves in the Human Brain. Cell Rep. 2019;26(2):374-380.e4. https://doi.org/10.1016/j.celrep.2018.12.058 PMID: 30625320

11. Nunez PL, Wingeier BM, Silberstein RB. Spatial-temporal structures of human alpha rhythms: theory, microcurrent sources, multiscale measurements, and global binding of local networks. Hum Brain Mapp. 2001;13(3):125–64. https://doi.org/10.1002/hbm.1030 PMID: 11376500

12. Ito J, Nikolaev AR, van Leeuwen C. Spatial and temporal structure of phase synchronization of spontaneous alpha EEG activity. Biol Cybern. 2005;92(1):54–60. https://doi.org/10.1007/s00422-004-0533-z PMID: 15650899

13. Alamia A, Terral L, D'ambra MR, VanRullen R. Distinct roles of forward and backward alpha-band waves in spatial visual attention. Elife. 2023;12:e85035. https://doi.org/10.7554/eLife.85035 PMID: 36876909

14. Pang 庞兆阳 Z, Alamia A, VanRullen R. Turning the Stimulus On and Off Changes the Direction of α Traveling Waves. eNeuro. 2020;7(6): ENEURO.0218-20.2020. https://doi.org/10.1523/ENEURO.0218-20.2020 PMID: 33168617
15. Zhigalov A, Jensen O. Perceptual echoes as travelling waves may arise from two discrete neuronal sources. Neuroimage. 2023;272:120047. https://doi.org/10.1016/j.neuroimage.2023.120047 PMID: 37001836

16. Bahramisharif A, van Gerven MAJ, Aarnoutse EJ, Mercier MR, Schwartz TH, Foxe JJ, et al. Propagating neocortical gamma bursts are coordinated by traveling alpha waves. J Neurosci. 2013;33(48):18849-54. https://doi.org/10.1523/JNEUROSCI.2455-13.2013 PMID: 24285891

17. Bhattacharya S, Brincat SL, Lundqvist M, Miller EK. Traveling waves in the prefrontal cortex during working memory. PLoS Comput Biol. 2022;18(1):e1009827. https://doi.org/10.1371/journal.pcbi.1009827 PMID: 35089915

18. Halgren M, Ulbert I, Bastuji H, Fabó D, Erőss L, Rey M, et al. The generation and propagation of the human alpha rhythm. Proc Natl Acad Sci U S A. 2019;116(47):23772-82. https://doi.org/10.1073/pnas.1913092116 PMID: 31685634

19. Alamia A, VanRullen R. Alpha oscillations and traveling waves: Signatures of predictive coding?. PLoS Biol. 2019;17(10):e3000487. https://doi.org/10.1371/journal.pbio.3000487 PMID: 31581198

20. Rao RP, Ballard DH. Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects. Nat Neurosci. 1999;2(1):79-87. https://doi.org/10.1038/4580 PMID: 10195184

21. Shipp S. The functional logic of cortico-pulvinar connections. Philos Trans R Soc Lond B Biol Sci. 2003;358(1438):1605-24. https://doi.org/10.1098/rstb.2002.1213 PMID: 14561322

22. Fiebelkorn IC, Kastner S. The Puzzling Pulvinar. Neuron. 2019;101(2):201-3. https://doi.org/10.1016/j.neuron.2018.12.032 PMID: 30653933

23. Jaramillo J, Mejias JF, Wang X-J. Engagement of Pulvino-cortical Feedforward and Feedback Pathways in Cognitive Computations. Neuron. 2019;101(2):321-336.e9. https://doi.org/10.1016/j.neuron.2018.11.023 PMID: 30553546

24. Quax S, Jensen O, Tiesinga P. Top-down control of cortical gamma-band communication via pulvinar induced phase shifts in the alpha rhythm. PLoS Comput Biol. 2017;13(5):e1005519. https://doi.org/10.1371/journal.pcbi.1005519 PMID: 28472057

25. Saalmann YB, Pinsk MA, Wang L, Li X, Kastner S. The pulvinar regulates information transmission between cortical areas based on attention demands. Science. 2012;337(6095):753-6. https://doi.org/10.1126/science.1223082 PMID: 22879517

26. Zhou H, Schafer RJ, Desimone R. Pulvinar-Cortex Interactions in Vision and Attention. Neuron. 2016;89(1):209-20. https://doi.org/10.1016/j.neuron.2015.11.034 PMID: 26748092

27. Felleman DJ, Van Essen DC. Distributed hierarchical processing in the primate cerebral cortex. Cereb Cortex. 1991;1(1):1-47. https://doi.org/10.1093/cercor/1.1.1-a PMID: 1822724

28. Markov NT, Vezoli J, Chameau P, Falchier A, Quilodran R, Huissoud C, et al. Anatomy of hierarchy: feedforward and feedback pathways in macaque visual cortex. J Comp Neurol. 2014;522(1):225-59. https://doi.org/10.1002/cne.23458 PMID: 23983048

29. Vezoli J, Magrou L, Goebel R, Wang X-J, Knoblauch K, Vinck M, et al. Cortical hierarchy, dual counterstream architecture and the importance of top-down generative networks. Neuroimage. 2021;225:117479. https://doi.org/10.1016/j.neuroimage.2020.117479 PMID: 33099005

30. Shipp S. Structure and function of the cerebral cortex. Curr Biol. 2007;17(12):R443-9. https://doi.org/10.1016/j.cub.2007.03.044 PMID: 17580069

31. Thomson AM, Bannister AP. Interlaminar connections in the neocortex. Cereb Cortex. 2003;13(1):5-14. https://doi.org/10.1093/cercor/13.1.5 PMID: 12466210

32. Maier A, Adams GK, Aura C, Leopold DA. Distinct superficial and deep laminar domains of activity in the visual cortex during rest and stimulation. Front Syst Neurosci. 2010;4:31. https://doi.org/10.3389/fnsys.2010.00031 PMID: 20802856

33. Maier A, Aura CJ, Leopold DA. Infragranular sources of sustained local field potential responses in macaque primary visual cortex. J Neurosci. 2011;31(6):1971-80. https://doi.org/10.1523/JNEUROSCI.5300-09.2011 PMID: 21307235

34. Mendoza-Halliday D, Major AJ, Lee N, Lichtenfeld MJ, Carlson B, Mitchell B, et al. A ubiquitous spectrolaminar motif of local field potential power across the primate cortex. Nat Neurosci. 2024;27(3):547-60. https://doi.org/10.1038/s41593-023-01554-7 PMID: 38238431

35. Steriade M, Gloor P, Llinás RR, Lopes de Silva FH, Mesulam MM. Report of IFCN Committee on Basic Mechanisms. Basic mechanisms of cerebral rhythmic activities. Electroencephalogr Clin Neurophysiol. 1990;76(6):481-508. https://doi.org/10.1016/0013-4694(90)90001-z PMID: 1701118

36. Sun W, Dan Y. Layer-specific network oscillation and spatiotemporal receptive field in the visual cortex. Proc Natl Acad Sci U S A. 2009;106(42):17986-91. https://doi.org/10.1073/pnas.0903962106 PMID: 19805197

37. Mejias JF, Murray JD, Kennedy H, Wang X-J. Feedforward and feedback frequency-dependent interactions in a large-scale laminar network of the primate cortex. Sci Adv. 2016;2(11):e1601335. https://doi.org/10.1126/sciadv.1601335 PMID: 28138530

38. Vijayan S, Kopell NJ. Thalamic model of awake alpha oscillations and implications for stimulus processing. Proc Natl Acad Sci U S A. 2012;109(45):18553-8. https://doi.org/10.1073/pnas.1215385109 PMID: 23054840

39. Gauthier B, Eger E, Hesselmann G, Giraud A-L, Kleinschmidt A. Temporal tuning properties along the human ventral visual stream. J Neurosci. 2012;32(41):14433-41. https://doi.org/10.1523/JNEUROSCI.2467-12.2012 PMID: 23055513

40. Murray JD, Bernacchia A, Freedman DJ, Romo R, Wallis JD, Cai X, et al. A hierarchy of intrinsic timescales across primate cortex. Nat Neurosci. 2014;17(12):1661-3. https://doi.org/10.1038/nn.3862 PMID: 25383900

41. Alamia A, Gordillo D, Chkonia E, Roinishvili M, Cappe C, Herzog MH. Oscillatory traveling waves reveal predictive coding abnormalities in schizophrenia. Cold Spring Harbor Laboratory. 2023. https://doi.org/10.1101/2023.10.09.561490
42. Bastos AM, Usrey WM, Adams RA, Mangun GR, Fries P, Friston KJ. Canonical microcircuits for predictive coding. Neuron. 2012;76(4):695-711.
https://doi.org/10.1016/j.neuron.2012.10.038 PMID: 23177956

43. Saalmann YB, Kastner S. Cognitive and perceptual functions of the visual thalamus. Neuron. 2011;71(2):209-23. https://doi.org/10.1016/j.neuron.2011.06.027 PMID: 21791281

44. Sherman SM. Thalamus plays a central role in ongoing cortical functioning. Nat Neurosci. 2016;19(4):533-41. https://doi.org/10.1038/nn.4269
PMID: 27021938

45. Grieve KL, Acuña C, Cudeiro J. The primate pulvinar nuclei: vision and action. Trends Neurosci. 2000;23(1):35-9. https://doi.org/10.1016/S0166-2236(99)01482-4 PMID: 10631787

46. Kaas JH, Lyon DC. Pulvinar contributions to the dorsal and ventral streams of visual processing in primates. Brain Res Rev. 2007;55(2):285-96.
https://doi.org/10.1016/j.brainresrev.2007.02.008 PMID: 17433837

47. Sherman SM, Guillery RW. The role of the thalamus in the flow of information to the cortex. Philos Trans R Soc Lond B Biol Sci. 2002;357(1428):1695-708. https://doi.org/10.1098/rstb.2002.1161 PMID: 12626004

48. Cortes N, Abbas Farishta R, Ladret HJ, Casanova C. Corticothalamic Projections Gate Alpha Rhythms in the Pulvinar. Front Cell Neurosci. 2021;15:787170. https://doi.org/10.3389/fncel.2021.787170 PMID: 34938163

49. Händel BF, Haarmeier T, Jensen O. Alpha oscillations correlate with the successful inhibition of unattended stimuli. J Cogn Neurosci. 2011;23(9):2494-502. https://doi.org/10.1162/jocn.2010.21557 PMID: 20681750

50. Sauseng P, Klimesch W, Stadler W, Schabus M, Doppelmayr M, Hanslmayr S, et al. A shift of visual spatial attention is selectively associated with human EEG alpha activity. Eur J Neurosci. 2005;22(11):2917-26. https://doi.org/10.1111/j.1460-9568.2005.04482.x PMID: 16324126

51. Thut G, Nietzel A, Brandt SA, Pascual-Leone A. Alpha-band electroencephalographic activity over occipital cortex indexes visuospatial attention bias and predicts visual target detection. J Neurosci. 2006;26(37):9494-502. https://doi.org/10.1523/JNEUROSCI.0875-06.2006 PMID: 16971533

52. Shipp S. Neural Elements for Predictive Coding. Front Psychol. 2016;7:1792. https://doi.org/10.3389/fpsyg.2016.01792 PMID: 27917138

53. Alexander DM, Hermens DF, Keage HAD, Clark CR, Williams LM, Kohn MR, et al. Event-related wave activity in the EEG provides new marker of ADHD. Clin Neurophysiol. 2008;119(1):163-79. https://doi.org/10.1016/j.clinph.2007.09.119 PMID: 18054279

54. Alexander DM, Flynn GJ, Wong W, Whitford TJ, Harris AWF, Galletly CA, et al. Spatio-temporal EEG waves in first episode schizophrenia. Clin Neurophysiol. 2009;120(9):1667-82. https://doi.org/10.1016/j.clinph.2009.06.020 PMID: 19646922

55. Patten TM, Rennie CJ, Robinson PA, Gong P. Human cortical traveling waves: dynamical properties and correlations with responses. PLoS One.
2012;7(6):e38392. https://doi.org/10.1371/journal.pone.0038392 PMID: 22675555

56. Aggarwal A, Brennan C, Luo J, Chung H, Contreras D, Kelz MB, et al. Visual evoked feedforward-feedback traveling waves organize neural activity across the cortical hierarchy in mice. Nat Commun. 2022;13(1):4754. https://doi.org/10.1038/s41467-022-32378-x PMID: 35963850

57. Davis ZW, Muller L, Martinez-Trujillo J, Sejnowski T, Reynolds JH. Spontaneous travelling cortical waves gate perception in behaving primates. Nature. 2020;587(7834):432-6. https://doi.org/10.1038/s41586-020-2802-y PMID: 33029013

58. Zabeh E, Foley NC, Jacobs J, Gottlieb JP. Beta traveling waves in monkey frontal and parietal areas encode recent reward history. Nat Commun. 2023;14(1):5428. https://doi.org/10.1038/s41467-023-41125-9 PMID: 37669966

59. Hallez H, Vanrumste B, Grech R, Muscat J, De Clercq W, Vergult A, et al. Review on solving the forward problem in EEG source analysis. J Neuroeng Rehabil. 2007;4:46. https://doi.org/10.1186/1743-0003-4-46 PMID: 18053144

60. Hindriks R, van Putten MJAM, Deco G. Intra-cortical propagation of EEG alpha oscillations. Neuroimage. 2014;103:444-53. https://doi.org/10.1016/j.neuroimage.2014.08.027 PMID: 25168275

61. Muller L, Chavane F, Reynolds J, Sejnowski TJ. Cortical travelling waves: mechanisms and computational principles. Nat Rev Neurosci. 2018;19(5):255-68. https://doi.org/10.1038/nrn.2018.20 PMID: 29563572

62. Lopes Da Silva FH, Storm Van Leeuwen W. The cortical source of the alpha rhythm. Neurosci Lett. 1977;6(2-3):237-41. https://doi.org/10.1016/0304-3940(77)90024-6 PMID: 19605058

63. Bollimunta A, Chen Y, Schroeder CE, Ding M. Neuronal mechanisms of cortical alpha oscillations in awake-behaving macaques. J Neurosci. 2008;28(40):9976-88. https://doi.org/10.1523/JNEUROSCI.2699-08.2008 PMID: 18829955

64. Bollimunta A, Mo J, Schroeder CE, Ding M. Neuronal mechanisms and attentional modulation of corticothalamic α oscillations. J Neurosci. 2011;31(13):4935-43. https://doi.org/10.1523/JNEUROSCI.5580-10.2011 PMID: 21451032

65. van Kerkooerle T, Self MW, Dagnino B, Gariel-Mathis M-A, Poort J, van der Togt C, et al. Alpha and gamma oscillations characterize feedback and feedforward processing in monkey visual cortex. Proc Natl Acad Sci U S A. 2014;111(40):14332-41. https://doi.org/10.1073/pnas.1402773111
PMID: 25205811

66. David O, Friston KJ. A neural mass model for MEG/EEG: coupling and neuronal dynamics. Neuroimage. 2003;20(3):1743-55. https://doi.org/10.1016/j.neuroimage.2003.07.015 PMID: 14642484

67. Jansen BH, Rit VG. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol Cybern. 1995;73(4):357-66. https://doi.org/10.1007/BF00199471 PMID: 7578475
68. Liley DTJ, Cadusch PJ, Dafilis MP. A spatially continuous mean field theory of electrocortical activity. Network: Computation in Neural Systems. 2002;13(1):67–113. https://doi.org/10.1080/net.13.1.67.113

69. Raj A, Cai C, Xie X, Palacios E, Owen J, Mukherjee P, et al. Spectral graph theory of brain oscillations. Hum Brain Mapp. 2020;41(11):2980–98. https://doi.org/10.1002/hbm.24991 PMID: 32202027

70. Lopes da Silva FH, Hoeks A, Smits H, Zetterberg LH. Model of brain rhythmic activity. The alpha-rhythm of the thalamus. Kybernetik. 1974;15(1):27–37. https://doi.org/10.1007/BF00270757 PMID: 4853232

71. Robinson PA, Rennie CJ, Rowe DL. Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. Phys Rev E Stat Nonlin Soft Matter Phys. 2002;65(4 Pt 1):041924. https://doi.org/10.1103/PhysRevE.65.041924 PMID: 12005890

72. Bastos AM, Vezoli J, Bosman CA, Schoffelen J-M, Oostenveld R, Dowdall JR, et al. Visual areas exert feedforward and feedback influences through distinct frequency channels. Neuron. 2015;85(2):390–401. https://doi.org/10.1016/j.neuron.2014.12.018 PMID: 25556836

73. Jensen O, Bonnefond M, Marshall TR, Tiesinga P. Oscillatory mechanisms of feedforward and feedback visual processing. Trends Neurosci. 2015;38(4):192–4. https://doi.org/10.1016/j.tins.2015.02.006 PMID: 25765320

74. Michalareas G, Vezoli J, van Pelt S, Schoffelen J-M, Kennedy H, Fries P. Alpha-Beta and Gamma Rhythms Subserve Feedback and Feedforward Influences among Human Visual Cortical Areas. Neuron. 2016;89(2):384–97. https://doi.org/10.1016/j.neuron.2015.12.018 PMID: 26777277

75. Spaak E, Bonnefond M, Maier A, Leopold DA, Jensen O. Layer-Specific Entrainment of Gamma-Band Neural Activity by the Alpha Rhythm in Monkey Visual Cortex. 2012;2313–8.

76. Canolty RT, Knight RT. The functional role of cross-frequency coupling. Trends Cogn Sci. 2010;14(11):506–15. https://doi.org/10.1016/j.tics.2010.09.001 PMID: 20932795

77. Sotero RC. Topology, Cross-Frequency, and Same-Frequency Band Interactions Shape the Generation of Phase-Amplitude Coupling in a Neural Mass Model of a Cortical Column. PLoS Comput Biol. 2016;12(11):e1005180. https://doi.org/10.1371/journal.pcbi.1005180 PMID: 27802274

78. Sotero RC, Bortel A, Naaman S, Mocanu VM, Kropf P, Villeneuve MY, et al. Laminar Distribution of Phase-Amplitude Coupling of Spontaneous Current Sources and Sinks. Front Neurosci. 2015;9:454. https://doi.org/10.3389/fnins.2015.00454 PMID: 26733778

79. Dantzker JL, Callaway EM. Laminar sources of synaptic input to cortical inhibitory interneurons and pyramidal neurons. Nat Neurosci. 2000;3(7):701–7. https://doi.org/10.1038/76656 PMID: 10862703

80. Xu X, Callaway EM. Laminar specificity of functional input to distinct types of inhibitory cortical neurons. J Neurosci. 2009;29(1):70–85. https://doi.org/10.1523/JNEUROSCI.4104-08.2009 PMID: 19129386

81. Chapeton JI, Haque R, Wittig JH Jr, Inati SK, Zaghloul KA. Large-Scale Communication in the Human Brain Is Rhythmically Modulated through Alpha Coherence. Curr Biol. 2019;29(17):2801-2811.e5. https://doi.org/10.1016/j.cub.2019.07.014 PMID: 31422882

82. Bonnefond M, Kastner S, Jensen O. Communication between Brain Areas Based on Nested Oscillations. eNeuro. 2017;4(2):ENEURO.0153-16.2017. https://doi.org/10.1523/ENEURO.0153-16.2017 PMID: 28374013

83. Shipp S, Adams RA, Friston KJ. Reflections on agranular architecture: predictive coding in the motor cortex. Trends Neurosci. 2013;36(12):706–16. https://doi.org/10.1016/j.tins.2013.09.004 PMID: 24157198

84. Bastos G, Holmes JT, Ross JM, Rader AM, Gallimore CG, Wargo JA, et al. Top-down input modulates visual context processing through an interneuron-specific circuit. Cell Rep. 2023;42(9):113133. https://doi.org/10.1016/j.celrep.2023.113133 PMID: 37708021

85. Gallimore CG, Ricci DA, Hamm JP. Spatiotemporal dynamics across visual cortical laminae support a predictive coding framework for interpreting mismatch responses. Cereb Cortex. 2023;33(15):9417–28. https://doi.org/10.1093/cercor/bhad215 PMID: 37310190

86. Hamm JP, Shymkiv Y, Han S, Yang W, Yuste R. Cortical ensembles selective for context. Proc Natl Acad Sci U S A. 2021;118(14):e2026179118. https://doi.org/10.1073/pnas.2026179118 PMID: 33811144

87. Jordan R, Keller GB. Opposing Influence of Top-down and Bottom-up Input on Excitatory Layer 2/3 Neurons in Mouse Primary Visual Cortex. Neuron. 2020;108(6):1194-1206.e5. https://doi.org/10.1016/j.neuron.2020.09.024 PMID: 33091338

88. Thomas ER, Haarsma J, Nicholson J, Yon D, Kok P, Press C. Predictions and errors are distinctly represented across V1 layers. Curr Biol. 2024;34(10):2265-2271.e4. https://doi.org/10.1016/j.cub.2024.04.036 PMID: 38697110

89. Lee K, Pennartz CMA, Mejias JF. Cortical networks with multiple interneuron types generate oscillatory patterns during predictive coding. Cold Spring Harbor Laboratory. 2024. https://doi.org/10.1101/2024.10.27.620494

90. Bastos AM, Lundqvist M, Waite AS, Kopell N, Miller EK. Layer and rhythm specificity for predictive routing. Proc Natl Acad Sci U S A. 2020;117(49):31459–69. https://doi.org/10.1073/pnas.2014868117 PMID: 33229572

91. Nejat H, Sherfey J, Bastos AM. Predictive routing emerges from self-supervised stochastic neural plasticity. bioRxiv. 2025;2024.12.31.630823. https://doi.org/10.1101/2024.12.31.630823 39803482

92. Benigno GB, Budzinski RC, Davis ZW, Reynolds JH, Muller L. Waves traveling over a map of visual space can ignite short-term predictions of sensory input. Nat Commun. 2023;14(1):3409. https://doi.org/10.1038/s41467-023-39076-2 PMID: 37296131

93. Crick F. Function of the thalamic reticular complex: the searchlight hypothesis. Proc Natl Acad Sci U S A. 1984;81(14):4586–90. https://doi.org/10.1073/pnas.81.14.4586 PMID: 6589612
94. Sillito AM, Jones HE. The role of the thalamic reticular nucleus in visual processing. THL. 2007;4(01). https://doi.org/10.1017/s1472928807000295

95. Destexhe A, Contreras D, Steriade M. Mechanisms underlying the synchronizing action of corticothalamic feedback through inhibition of thalamic relay cells. J Neurophysiol. 1998;79(2):999-1016. https://doi.org/10.1152/jn.1998.79.2.999 PMID: 9463458

96. Abbas Farishta R, Boire D, Casanova C. Hierarchical Organization of Corticothalamic Projections to the Pulvinar. Cereb Cortex Commun. 2020;1(1):tgaa030. https://doi.org/10.1093/texcom/tgaa030 PMID: 34296104

97. Sherman SM. The thalamus is more than just a relay. Curr Opin Neurobiol. 2007;17(4):417-22. https://doi.org/10.1016/j.conb.2007.07.003 PMID: 17707635

98. Ogren MP, Hendrickson AE. The distribution of pulvinar terminals in visual areas 17 and 18 of the monkey. Brain Res. 1977;137(2):343-50. https://doi.org/10.1016/0006-8993(77)90344-4 PMID: 412565

99. Purushothaman G, Marion R, Li K, Casagrande VA. Gating and control of primary visual cortex by pulvinar. Nat Neurosci. 2012;15(6):905-12. https://doi.org/10.1038/nn.3106 PMID: 22561455

100. Cortes N, Ladret HJ, Abbas-Farishta R, Casanova C. The pulvinar as a hub of visual processing and cortical integration. Trends Neurosci. 2024;47(2):120-34. https://doi.org/10.1016/j.tins.2023.11.008 PMID: 38143202

101. O'Reilly RC, Russin JL, Zolfaghar M, Rohrlich J. Deep Predictive Learning in Neocortex and Pulvinar. J Cogn Neurosci. 2021;33(6):1158-96. https://doi.org/10.1162/jocn_a_01708 PMID: 34428793

102. Izhikevich EM. Simple model of spiking neurons. IEEE Trans Neural Netw. 2003;14(6):1569-72. https://doi.org/10.1109/TNN.2003.820440 PMID: 18244602

103. Oostenveld R, Fries P, Maris E, Schoffelen J-M. FieldTrip: Open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data. Comput Intell Neurosci. 2011;2011:156869. https://doi.org/10.1155/2011/156869 PMID: 21253357

104. Tzourio-Mazoyer N, Landeau B, Papathanassiou D, Crivello F, Etard O, Delcroix N, et al. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. Neuroimage. 2002;15(1):273-89. https://doi.org/10.1006/nimg.2001.0978 PMID: 11771995