It is well established that the brain’s cortical areas differ in their cyto- and myeloarchitectonic structure, local and long-range anatomical connectivity, activity and, by consequence, their function (Glasser et al., 2016; Huntenburg et al., 2017). Interestingly, many structural features that distinguish individual brain areas change gradually in an orderly manner across the cortex, leading to spatial feature gradients. The most prominent and best established gradients are evident along the posterior-anterior axis (Eickhoff et al., 2018; Felleman and Van Essen, 1991; Huntenburg et al., 2018). For instance, neuron density decreases and neuronal connectivity increases from posterior to anterior brain areas. These differences have been attributed to differences in neurogenesis for posterior compared to anterior brain areas (Hill et al., 2010; Huntenburg et al., 2018). A similar posterior-anterior gradient has been observed for myelin content, cortical thickness, and gene expression (Burt et al., 2018). Next to the posterior-anterior gradient, other global spatial organization principles have been proposed to explain the variation of microstructural features across the cortex. For instance, Huntenburg et al. suggest a sensorimotor to transmodal gradient as an important intrinsic organizing dimension of human cortex (Huntenburg et al., 2018) reflecting gradual changes in structural features from functionally unimodal (dedicated sensory or motor) areas to higher order, transmodal areas.

In addition to structural gradients as an organizing principle reflecting global cortical organization, it is well acknowledged that cortical areas are structurally connected into larger networks, which often display a hierarchical organization. Cortical hierarchies are typically established based on the degree of microstructural differentiation of the connected areas, and on the classification of the anatomical connections as feedforward or feedback using histological tract-tracing. Early sensory areas with predominantly feedforward outgoing connections are placed at the bottom of the hierarchy and higher order association areas with mostly feedback outgoing connections are placed at the top of the hierarchy (Felleman and Van Essen, 1991; Markov et al., 2014). A noninvasive, but indirect index of these hierarchies is cortical thickness, a macroscopic feature of the cortex, which can be estimated from MRI scans. It has been shown that cortical thickness mirrors global hierarchical organization of the cortex as well as local hierarchies in visual, auditory and somatosensory areas (Jasmin et al., 2019; Wagstyl et al., 2015), and, therefore, could be used as a basis for understanding hierarchy-gradient relationships in the cortex.

The presence of these anatomical gradients raises the question to what extent they are reflected in features of brain activity and brain function. Indeed, it has been shown that cortical areas follow a hierarchical ordering in their timescales of intrinsic fluctuations as for example measured in the autocorrelation of spiking activity (Murray et al., 2014). Sensory areas show faster fluctuations while frontal areas show slower fluctuations. Shorter timescales in sensory areas enables them to reflect dynamic changes in the environment, whereas the longer timescales in prefrontal areas allows for integration of information. Particularly, this gradient of ‘temporal receptive windows’ has been demonstrated in visual (Himberger et al., 2018) and auditory processing (Jasmin et al., 2019) and could be related to the frequency of spontaneous brain oscillations. Oscillations are a prominent feature of brain activity, and have been suggested to play a central role in coordinating neuronal activity (Fries, 2005; Wang, 2010). Similar to many anatomical features described above, the spectral activity patterns have been shown to be characteristic for each brain area (Keitel and Gross, 2016). This is consistent with the view that the individual anatomical structure of a brain area shapes its rhythmic neuronal activity, which led us to hypothesize the existence of a posterior-anterior gradient in the frequency of spontaneous brain rhythms.

Spontaneous rhythms have been studied in the past but typically by focusing on the power in specific frequency bands (Hillebrand et al., 2016; Keitel and Gross, 2016; Mellem et al., 2017). Overall, these MEG studies revealed strongest cortical generators for the dominant alpha rhythm (7–13 Hz) in occipito-parietal brain areas. The beta band (15–30 Hz) shows strongest activity in sensorimotor areas while delta (1–3 Hz) and theta (3–7 Hz) bands are associated with activity in wide-spread areas including frontal cortex.

Here, we adopt a different approach that is based on sophisticated identification of spectral peaks in the power spectra of source-localized resting-state MEG data and included modelling of the 1/f spectral background (Haller et al., 2018). This approach offers two distinct advantages. First, focusing on spectral peaks ensures that results are indeed based on brain oscillations. This is not necessarily the case when using the power in a pre-defined frequency band or using band-pass filtered data. Second, by explicitly modeling the 1/f spectral background across the entire cortex we can dissociate contributions due to aperiodic neuronal background activity from those originating from oscillatory activity.

We used this approach to specifically test the hypothesis of a posterior-anterior gradient in the frequency of spontaneous brain rhythms. We identified the frequencies of the dominant brain rhythm across the cortex in source-localized resting-state MEG data of 187 individuals.

As we describe below, we found spatial gradients of the dominant peak frequency across the cortex following the cortical hierarchy.

This study is the first comprehensive demonstration of frequency gradients across the human cortex using a large set of resting-state MEG recordings. We found that the dominant peak frequency in a brain area decreases significantly, gradually and robustly along the posterior-anterior axis, following the global cortical hierarchy from early sensory to higher order areas. This finding establishes a frequency gradient of resting-state brain rhythms that complements previous anatomical studies reporting a posterior-anterior gradient in microscale and macroscale anatomical features of animal and human cortex (Huntenburg et al., 2018).

Several MEG/EEG studies have demonstrated that alpha activity (~10 Hz) is strongest in occipito-parietal brain areas and theta activity (~5 Hz) strongest in more frontal brain areas (Chiang et al., 2011; Voytek et al., 2010). There is also evidence that the dominant frequencies differ between brain areas (Groppe et al., 2013; Hillebrand et al., 2016; Keitel and Gross, 2016). However, we provide the first comprehensive (in frequency and space) statistical model of frequency gradients in a large resting-state brain activity.

Our results are consistent with a recent invasive study showing a systematic decrease of peak frequency from posterior to anterior brain areas in ECoG recordings of epilepsy patients (Zhang et al., 2018).

A differentiating feature of our approach was that we used a large number of healthy participants (N = 187), reconstructed cortical activity from noninvasive MEG recordings, and considered further anatomical features (i.e. cortical thickness). Notably, estimating the power spectrum in finely parcellated ROIs allowed us an accurate and robust identification of peak frequencies and characterization of their spatial gradients across the entire cortical surface. Importantly, we focus on peaks in the power spectrum that indicate the presence of rhythmicity in the neuronal activity, instead of focusing on predefined frequency bands where these rhythms might be absent. As slope and offset of frequency gradient could dramatically vary across participants, averaging across participants may not yield a reliable representation of PF gradient. Instead, we used mixed effect modelling of PF along the cortical hierarchies, where the between-participant variability was taken into account as a random effect. For example, the peak frequency in occipital brain areas and the slope of the frequency gradient varies across individuals, and thus the intercept and slope of the corresponding least square fit change. These individual differences are specifically modelled by a mixed effect model as random effects. Importantly, a mixed effect model will only show a significant gradient if it is significant across cortical areas at the individual level, and consistent across participants at the group level.

Results of our analyses showed that, just as peak frequency significantly decreased along the posterior-anterior axis, CT significantly increased in the same direction, which resulted in a significant anticorrelation between PF and CT. The observed correlation holds after removing the effect of spatial location (x,y,z). This seems to indicate that PF and CT are more closely related to each other than can be explained by spatial location alone. Since cortical hierarchies do not strictly follow a single linear trajectory in space (e.g. posterior-anterior) our results are consistent with the idea that both PF and CT, follow cortical hierarchies. Indeed, such local spatial gradients have been reported in multiple features of cortex during auditory perception (Jasmin et al., 2019) and visual processing streams (Himberger et al., 2018). From a broader view, local gradients could mirror complex organization of gradients in human cortex and support the approach of global gradient along the sensory to transmodal areas (Huntenburg et al., 2018). On the other hand, posterior-anterior gradient of PF was significant after subtracting CT scores from PF values. This suggests a partial independence of both measures. Since PF is a measure derived from brain activity the reported gradient could be modulated dynamically depending on cognitive state or task demands. Further studies are needed to investigate this in more detail.

We further addressed the question if our results can be explained by the linear superposition of activity from an occipital alpha source and a frontal theta source. Along the posterior-anterior axis, differential superposition of both sources could lead to a frequency gradient, due to imperfect unmixing of the signals. However, our analysis revealed that a significant frequency gradient is already evident within 2–3 cm of V1 where the effect of a frontal theta source (which has on average a lower power compared to occipital alpha) is negligible. Furthermore, A more comprehensive analysis showed existence of this gradient for posterio-pariatal, central, and frontal areas, separately, acknowledging the globally and continuously decreasing nature of the posterior-anterior gradient.

Additional supporting evidence can be drawn from intracranial studies, where the data is directly recorded from cortex. Zhang et al., 2018 have shown that oscillations generally propagate in a posterior-to-anterior direction because they are coordinated by an overall decrease in intrinsic oscillation frequency from posterior to anterior regions (see Figures S6 and 7 of Zhang et al., 2018). Overall, this indicates the existence of a gradual decrease of PF along the posterior-anterior axis.

What is the potential functional role of this frequency gradient? Zhang et al. demonstrated the existence of travelling waves along the frequency gradient (Zhang et al., 2018). Interestingly, they found that local frequencies along the posterior-to-anterior direction are positively correlated with waves’ propagation speed and direction consistent with a proposed model of travelling waves based on weakly coupled oscillators (WCO) (Ermentrout and Kleinfeld, 2001). Travelling waves along the posterior-anterior axis have also been reported during visual stimulation (Alamia and VanRullen, 2019; Lozano-Soldevilla and VanRullen, 2019). These travelling waves might serve to drive neural communication along the cortical hierarchy possibly through nested gamma oscillations thereby linking travelling waves to the concept of pulsed inhibition (Bahramisharif et al., 2013). In addition, travelling waves have been associated with memory consolidation and learning (Muller et al., 2018). It is of interest to note that frequency gradients have been reported previously in the entorhinal cortex (Giocomo et al., 2011; Giocomo and Hasselmo, 2009). Here, a frequency decrease and corresponding travelling waves have been observed in the dorsal-ventral direction and have been related to a representational gradient of spatial scales from coarse to fine (Muller et al., 2018). Indeed, converging evidence across recording methods, species and cortical domains suggests that representations become more ‘integrated’ with decreasing ‘dominant’ frequency of the underlying neuronal population. A prime example is the auditory cortex where response latencies and complexity of processing increase along the posterior-anterior axis (Jasmin et al., 2019). This is also mirrored by an increase in cortical thickness and increased ratio of feedback to feedforward connections along this axis. Similar observations have been made across more widely distributed cortical areas where timescales of intrinsic fluctuations in spiking activity increase from posterior to anterior brain areas (Murray et al., 2014). Not surprisingly, these time scales are largely determined by the time constants of synaptic transmission (Duarte et al., 2017). But interestingly, in a computational model of activity in macaque cortex using anatomical connectivity a gradient of time scales also emerges with short, transient responses to input in sensory areas and slower, sustained responses in higher order areas (Chaudhuri et al., 2015) (see also Kiebel et al., 2008).

Our approach additionally revealed that cortical peak frequencies decrease systematically along the inferior-superior axis. As seen in Figure 1 this seems to result from the fact that higher order frontal areas with lower PF have higher z-coordinates compared to the early sensory areas with higher PF. A similar inferior-superior gradient was also observed for cortical thickness. This finding is supported by a recent study reporting an inferior-superior gradient organization of the CT in human and macaque monkeys (Valk et al., 2020), and has been attributed to the organization patterns expected based on the theory of dual origin (Goulas et al., 2018).

Our detailed analysis was based on the cortical ROIs’ spectral peak with strongest power (PF). However, we identified all peaks in the power spectrum of each ROI. Since spectral peaks indicate the presence of brain rhythms, this data represents a comprehensive overview of these rhythms across the cortex. The histogram of spectral peaks across ROIs and participants provided a data-driven definition of frequency bands. Interestingly, the histogram delineates the classical frequency bands with histogram peaks centering at 4–7.5 (theta), 8.5–13 (alpha), 15–25 (low-beta) 27.5–34 (high-beta) (see Figure 5). This is the first MEG study to our knowledge to identify frequency bands from peak frequencies in a large data set (see Groppe et al., 2013 for a similar approach in a smaller sample of ECoG data).

We further analyzed these specific frequency bands for gradients and found significant posterior-anterior frequency changes in the theta, alpha and beta frequency band. Results in the alpha band mirrored the previous results based on the overall strongest peak frequency. Interestingly, and in contrast to the alpha band, peak frequencies increased along the posterior-anterior direction in the theta and beta frequency band. In the model used by Zhang et al. this would correspond to travelling waves from anterior to posterior brain areas (Zhang et al., 2018) that might represent frequency channels for top-down effects (Michalareas et al., 2016; Wang, 2010).

An alternative explanation for the observed posterior-to-anterior or anterior-to-posterior changes of the band-specific PFs may come from the laminar organization of the cortex, where several layers exhibit distinct frequency profiles (Bastos et al., 2018) and thickness patterns (Wagstyl et al., 2020). In particular, similar to our results for the spatial gradients of band-specific PFs, Wagstyl and colleagues (Wagstyl et al., 2020) have identified both increasing and decreasing gradients of thickness along the posterior-anterior axis for cortical layers, in the somatosensory, auditory, and motor cortex (see Figure 6 of Wagstyl et al., 2020).

In summary, our findings show that peak frequencies of cortical areas form a spatial gradient, which follows the global posterior-anterior hierarchy as well as local anatomical hierarchies. Previous research also points to spatial gradients in multiple features of the human and animal cortex. Further research might explore implications of frequency gradients in different cognitive states, disease, and aging.

We used a publicly available dataset for this study (https://data.donders.ru.nl/collections/di/dccn/DSC_3011020.09_236?0).

1. 1. Alamia A
   2. VanRullen R

   (2019) Alpha oscillations and traveling waves: signatures of predictive coding?

   PLOS Biology 17:e3000487.

   https://doi.org/10.1371/journal.pbio.3000487

   • Google Scholar
2. 1. Bahramisharif A
   2. van Gerven MA
   3. Aarnoutse EJ
   4. Mercier MR
   5. Schwartz TH
   6. Foxe JJ
   7. Ramsey NF
   8. Jensen O

   (2013) Propagating neocortical gamma bursts are coordinated by traveling alpha waves

   Journal of Neuroscience 33:18849–18854.

   https://doi.org/10.1523/JNEUROSCI.2455-13.2013

   • PubMed
   • Google Scholar
3. 1. Bastos AM
   2. Loonis R
   3. Kornblith S
   4. Lundqvist M
   5. Miller EK

   (2018) Laminar recordings in frontal cortex suggest distinct layers for maintenance and control of working memory

   PNAS 115:1117–1122.

   https://doi.org/10.1073/pnas.1710323115

   • PubMed
   • Google Scholar
4. 1. Bell AJ
   2. Sejnowski TJ

   (1995) An information-maximization approach to blind separation and blind deconvolution

   Neural Computation 7:1129–1159.

   https://doi.org/10.1162/neco.1995.7.6.1129

   • PubMed
   • Google Scholar
5. Book

   1. Brodmann K

   (1909)

   Vergleichende Lokalisationslehre Der Grosshirnrinde in Ihren Prinzipien Dargestellt Auf Grund Des Zellenbaues

   Johann Ambrosius Barth.

   • Google Scholar
6. 1. Burt JB
   2. Demirtaş M
   3. Eckner WJ
   4. Navejar NM
   5. Ji JL
   6. Martin WJ
   7. Bernacchia A
   8. Anticevic A
   9. Murray JD

   (2018) Hierarchy of transcriptomic specialization across human cortex captured by structural neuroimaging topography

   Nature Neuroscience 21:1251–1259.

   https://doi.org/10.1038/s41593-018-0195-0

   • PubMed
   • Google Scholar
7. 1. Chaudhuri R
   2. Knoblauch K
   3. Gariel MA
   4. Kennedy H
   5. Wang XJ

   (2015) A Large-Scale circuit mechanism for hierarchical dynamical processing in the primate cortex

   Neuron 88:419–431.

   https://doi.org/10.1016/j.neuron.2015.09.008

   • PubMed
   • Google Scholar
8. 1. Chiang AK
   2. Rennie CJ
   3. Robinson PA
   4. van Albada SJ
   5. Kerr CC

   (2011) Age trends and sex differences of alpha rhythms including split alpha peaks

   Clinical Neurophysiology 122:1505–1517.

   https://doi.org/10.1016/j.clinph.2011.01.040

   • PubMed
   • Google Scholar
9. 1. Demirtaş M
   2. Burt JB
   3. Helmer M
   4. Ji JL
   5. Adkinson BD
   6. Glasser MF
   7. Van Essen DC
   8. Sotiropoulos SN
   9. Anticevic A
   10. Murray JD

   (2019) Hierarchical heterogeneity across human cortex shapes Large-Scale neural dynamics

   Neuron 101:1181–1194.

   https://doi.org/10.1016/j.neuron.2019.01.017

   • PubMed
   • Google Scholar
10. 1. Duarte R
    2. Seeholzer A
    3. Zilles K
    4. Morrison A

    (2017) Synaptic patterning and the timescales of cortical dynamics

    Current Opinion in Neurobiology 43:156–165.

    https://doi.org/10.1016/j.conb.2017.02.007

    • PubMed
    • Google Scholar
11. 1. Eickhoff SB
    2. Constable RT
    3. Yeo BTT

    (2018) Topographic organization of the cerebral cortex and brain cartography

    NeuroImage 170:332–347.

    https://doi.org/10.1016/j.neuroimage.2017.02.018

    • PubMed
    • Google Scholar
12. 1. Ermentrout GB
    2. Kleinfeld D

    (2001) Traveling electrical waves in cortex: insights from phase dynamics and speculation on a computational role

    Neuron 29:33–44.

    https://doi.org/10.1016/s0896-6273(01)00178-7

    • PubMed
    • Google Scholar
13. 1. Felleman DJ
    2. Van Essen DC

    (1991) Distributed hierarchical processing in the primate cerebral cortex

    Cerebral Cortex 1:1–47.

    https://doi.org/10.1093/cercor/1.1.1

    • PubMed
    • Google Scholar
14. 1. Fries P

    (2005) A mechanism for cognitive dynamics: neuronal communication through neuronal coherence

    Trends in Cognitive Sciences 9:474–480.

    https://doi.org/10.1016/j.tics.2005.08.011

    • PubMed
    • Google Scholar
15. 1. Giocomo LM
    2. Moser MB
    3. Moser EI

    (2011) Computational models of grid cells

    Neuron 71:589–603.

    https://doi.org/10.1016/j.neuron.2011.07.023

    • PubMed
    • Google Scholar
16. 1. Giocomo LM
    2. Hasselmo ME

    (2009) Knock-out of HCN1 subunit flattens dorsal-ventral frequency gradient of medial entorhinal neurons in adult mice

    Journal of Neuroscience 29:7625–7630.

    https://doi.org/10.1523/JNEUROSCI.0609-09.2009

    • PubMed
    • Google Scholar
17. 1. Glasser MF
    2. Coalson TS
    3. Robinson EC
    4. Hacker CD
    5. Harwell J
    6. Yacoub E
    7. Ugurbil K
    8. Andersson J
    9. Beckmann CF
    10. Jenkinson M
    11. Smith SM
    12. Van Essen DC

    (2016) A multi-modal parcellation of human cerebral cortex

    Nature 536:171–178.

    https://doi.org/10.1038/nature18933

    • PubMed
    • Google Scholar
18. 1. Goulas A
    2. Zilles K
    3. Hilgetag CC

    (2018) Cortical gradients and laminar projections in mammals

    Trends in Neurosciences 41:775–788.

    https://doi.org/10.1016/j.tins.2018.06.003

    • PubMed
    • Google Scholar
19. 1. Groppe DM
    2. Bickel S
    3. Keller CJ
    4. Jain SK
    5. Hwang ST
    6. Harden C
    7. Mehta AD

    (2013) Dominant frequencies of resting human brain activity as measured by the electrocorticogram

    NeuroImage 79:223–233.

    https://doi.org/10.1016/j.neuroimage.2013.04.044

    • PubMed
    • Google Scholar
20. Preprint

    1. Haller M
    2. Donoghue T
    3. Peterson E
    4. Varma P
    5. Sebastian P
    6. Gao R
    7. Noto T
    8. Knight RT
    9. Shestyuk A
    10. Voytek B

    (2018) Parameterizing neural power spectra

    bioRxiv.

    https://doi.org/10.1101/299859

    • Google Scholar
21. 1. Hill J
    2. Inder T
    3. Neil J
    4. Dierker D
    5. Harwell J
    6. Van Essen D

    (2010) Similar patterns of cortical expansion during human development and evolution

    PNAS 107:13135–13140.

    https://doi.org/10.1073/pnas.1001229107

    • PubMed
    • Google Scholar
22. 1. Hillebrand A
    2. Tewarie P
    3. van Dellen E
    4. Yu M
    5. Carbo EW
    6. Douw L
    7. Gouw AA
    8. van Straaten EC
    9. Stam CJ

    (2016) Direction of information flow in large-scale resting-state networks is frequency-dependent

    PNAS 113:3867–3872.

    https://doi.org/10.1073/pnas.1515657113

    • PubMed
    • Google Scholar
23. 1. Himberger KD
    2. Chien HY
    3. Honey CJ

    (2018) Principles of temporal processing across the cortical hierarchy

    Neuroscience 389:161–174.

    https://doi.org/10.1016/j.neuroscience.2018.04.030

    • PubMed
    • Google Scholar
24. 1. Huntenburg JM
    2. Bazin PL
    3. Goulas A
    4. Tardif CL
    5. Villringer A
    6. Margulies DS

    (2017) A systematic relationship between functional connectivity and intracortical myelin in the human cerebral cortex

    Cerebral Cortex 27:981–997.

    https://doi.org/10.1093/cercor/bhx030

    • PubMed
    • Google Scholar
25. 1. Huntenburg JM
    2. Bazin PL
    3. Margulies DS

    (2018) Large-Scale gradients in human cortical organization

    Trends in Cognitive Sciences 22:21–31.

    https://doi.org/10.1016/j.tics.2017.11.002

    • PubMed
    • Google Scholar
26. 1. Ito T
    2. Kulkarni KR
    3. Schultz DH
    4. Mill RD
    5. Chen RH
    6. Solomyak LI
    7. Cole MW

    (2017) Cognitive task information is transferred between brain regions via resting-state network topology

    Nature Communications 8:1027.

    https://doi.org/10.1038/s41467-017-01000-w

    • PubMed
    • Google Scholar
27. 1. Jasmin K
    2. Lima CF
    3. Scott SK

    (2019) Understanding rostral-caudal auditory cortex contributions to auditory perception

    Nature Reviews Neuroscience 20:425–434.

    https://doi.org/10.1038/s41583-019-0160-2

    • PubMed
    • Google Scholar
28. 1. Keitel A
    2. Gross J

    (2016) Individual human brain Areas can be identified from their characteristic spectral activation fingerprints

    PLOS Biology 14:e1002498.

    https://doi.org/10.1371/journal.pbio.1002498

    • PubMed
    • Google Scholar
29. 1. Kiebel SJ
    2. Daunizeau J
    3. Friston KJ

    (2008) A hierarchy of time-scales and the brain

    PLOS Computational Biology 4:e1000209.

    https://doi.org/10.1371/journal.pcbi.1000209

    • PubMed
    • Google Scholar
30. 1. Lozano-Soldevilla D
    2. VanRullen R

    (2019) The hidden spatial dimension of alpha: 10-hz perceptual echoes propagate as periodic traveling waves in the human brain

    Cell Reports 26:374–380.

    https://doi.org/10.1016/j.celrep.2018.12.058

    • PubMed
    • Google Scholar
31. 1. Markov NT
    2. Vezoli J
    3. Chameau P
    4. Falchier A
    5. Quilodran R
    6. Huissoud C
    7. Lamy C
    8. Misery P
    9. Giroud P
    10. Ullman S
    11. Barone P
    12. Dehay C
    13. Knoblauch K
    14. Kennedy H

    (2014) Anatomy of hierarchy: feedforward and feedback pathways in macaque visual cortex

    Journal of Comparative Neurology 522:225–259.

    https://doi.org/10.1002/cne.23458

    • PubMed
    • Google Scholar
32. 1. Mellem MS
    2. Wohltjen S
    3. Gotts SJ
    4. Ghuman AS
    5. Martin A

    (2017) Intrinsic frequency biases and profiles across human cortex

    Journal of Neurophysiology 118:2853–2864.

    https://doi.org/10.1152/jn.00061.2017

    • PubMed
    • Google Scholar
33. 1. Michalareas G
    2. Vezoli J
    3. van Pelt S
    4. Schoffelen JM
    5. Kennedy H
    6. Fries P

    (2016) Alpha-Beta and gamma rhythms subserve feedback and feedforward influences among human visual cortical Areas

    Neuron 89:384–397.

    https://doi.org/10.1016/j.neuron.2015.12.018

    • PubMed
    • Google Scholar
34. 1. Muller L
    2. Chavane F
    3. Reynolds J
    4. Sejnowski TJ

    (2018) Cortical travelling waves: mechanisms and computational principles

    Nature Reviews Neuroscience 19:255–268.

    https://doi.org/10.1038/nrn.2018.20

    • PubMed
    • Google Scholar
35. 1. Murray JD
    2. Bernacchia A
    3. Freedman DJ
    4. Romo R
    5. Wallis JD
    6. Cai X
    7. Padoa-Schioppa C
    8. Pasternak T
    9. Seo H
    10. Lee D
    11. Wang XJ

    (2014) A hierarchy of intrinsic timescales across primate cortex

    Nature Neuroscience 17:1661–1663.

    https://doi.org/10.1038/nn.3862

    • PubMed
    • Google Scholar
36. 1. Oostenveld R
    2. Fries P
    3. Maris E
    4. Schoffelen JM

    (2011) FieldTrip: open source software for advanced analysis of MEG, EEG, and invasive electrophysiological data

    Computational Intelligence and Neuroscience 2011:1–9.

    https://doi.org/10.1155/2011/156869

    • PubMed
    • Google Scholar
37. Book

    1. Penny WD
    2. Friston KJ
    3. Ashburner JT
    4. Kiebel SJ
    5. Nichols TE

    (2011)

    Statistical Parametric Mapping: The Analysis of Functional Brain Images

    Elsevier Academic Press.

    • Google Scholar
38. 1. Pernet CR
    2. Wilcox R
    3. Rousselet GA

    (2012) Robust correlation analyses: false positive and power validation using a new open source matlab toolbox

    Frontiers in Psychology 3:606.

    https://doi.org/10.3389/fpsyg.2012.00606

    • PubMed
    • Google Scholar
39. 1. Schoffelen JM
    2. Hultén A
    3. Lam N
    4. Marquand AF
    5. Uddén J
    6. Hagoort P

    (2017) Frequency-specific directed interactions in the human brain network for language

    PNAS 114:8083–8088.

    https://doi.org/10.1073/pnas.1703155114

    • PubMed
    • Google Scholar
40. 1. Schoffelen JM
    2. Oostenveld R
    3. Lam NHL
    4. Uddén J
    5. Hultén A
    6. Hagoort P

    (2019) A 204-subject multimodal neuroimaging dataset to study language processing

    Scientific Data 6:17.

    https://doi.org/10.1038/s41597-019-0020-y

    • PubMed
    • Google Scholar
41. Preprint

    1. Valk SL
    2. Xu T
    3. Margulies DS
    4. Kahrabian Masouleh S
    5. Paquola C
    6. Goulas A
    7. Kochunov P
    8. Smallwood J
    9. Yeo BTT
    10. Bernhardt BC
    11. Eickhoff SB

    (2020) Shaping brain structure: genetic and phylogenetic axes of macro scale organization of cortical thickness

    bioRxiv.

    https://doi.org/10.1101/2020.02.10.939561

    • Google Scholar
42. 1. Van Essen DC
    2. Glasser MF
    3. Dierker DL
    4. Harwell J
    5. Coalson T

    (2012) Parcellations and hemispheric asymmetries of human cerebral cortex analyzed on surface-based atlases

    Cerebral Cortex 22:2241–2262.

    https://doi.org/10.1093/cercor/bhr291

    • PubMed
    • Google Scholar
43. 1. Van Veen BD
    2. van Drongelen W
    3. Yuchtman M
    4. Suzuki A

    (1997) Localization of brain electrical activity via linearly constrained minimum variance spatial filtering

    IEEE Transactions on Biomedical Engineering 44:867–880.

    https://doi.org/10.1109/10.623056

    • PubMed
    • Google Scholar
44. 1. Voytek B
    2. Canolty RT
    3. Shestyuk A
    4. Crone NE
    5. Parvizi J
    6. Knight RT

    (2010) Shifts in gamma phase-amplitude coupling frequency from theta to alpha over posterior cortex during visual tasks

    Frontiers in Human Neuroscience 4:191.

    https://doi.org/10.3389/fnhum.2010.00191

    • PubMed
    • Google Scholar
45. 1. Wagstyl K
    2. Ronan L
    3. Goodyer IM
    4. Fletcher PC

    (2015) Cortical thickness gradients in structural hierarchies

    NeuroImage 111:241–250.

    https://doi.org/10.1016/j.neuroimage.2015.02.036

    • PubMed
    • Google Scholar
46. 1. Wagstyl K
    2. Larocque S
    3. Cucurull G
    4. Lepage C
    5. Cohen JP
    6. Bludau S
    7. Palomero-Gallagher N
    8. Lewis LB
    9. Funck T
    10. Spitzer H
    11. Dickscheid T
    12. Fletcher PC
    13. Romero A
    14. Zilles K
    15. Amunts K
    16. Bengio Y
    17. Evans AC

    (2020) BigBrain 3D atlas of cortical layers: cortical and laminar thickness gradients diverge in sensory and motor cortices

    PLOS Biology 18:e3000678.

    https://doi.org/10.1371/journal.pbio.3000678

    • PubMed
    • Google Scholar
47. 1. Wang XJ

    (2010) Neurophysiological and computational principles of cortical rhythms in cognition

    Physiological Reviews 90:1195–1268.

    https://doi.org/10.1152/physrev.00035.2008

    • PubMed
    • Google Scholar
48. 1. Zhang H
    2. Watrous AJ
    3. Patel A
    4. Jacobs J

    (2018) Theta and alpha oscillations are traveling waves in the human neocortex

    Neuron 98:1269–1281.

    https://doi.org/10.1016/j.neuron.2018.05.019

    • PubMed
    • Google Scholar

Acceptance summary:

The manuscript makes an important contribution to the current literature showing how the frequency gradients of human resting-state neuronal oscillations are organized in accordance with cortical hierarchy.

Decision letter after peer review:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Thank you for submitting your work entitled "The frequency gradient of human resting-state brain oscillations follows cortical hierarchies" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Laura Dugué as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by a Senior Editor.

Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife.

This manuscript presents results from a large dataset of human resting-state MEG recordings. The authors investigated the timely topic of the spatial distribution of low frequency brain oscillations, and found an anterior-posterior gradient in the frequency of brain oscillations. They further show a gradient in cortical thickness (CT), and suggest that CT correlates with oscillations frequency. Although investigating the spatial distribution of brain oscillations goes beyond simply studying their temporal dynamics, other studies have shown before that oscillation frequency varies across the brain. Reviewers have also noticed that such spatial gradients are in fact more specific to the alpha frequency. Some methodological aspects have additionally been raised including addressing the question of spatial leakage, and the possible confound of inter individual differences. Based on these comments and others described below, and on following discussions, the three reviewers concluded that this manuscript will not be considered further for publication in eLife.

Reviewer #1:

Mahjoory et al. present new results on a publicly available, resting-state MEG recording dataset from a large pool of participants. This well-written manuscript provides a systematic analysis of the spatial distribution of high-power, low-frequency, oscillations, based on non-invasive recordings of healthy, human participants. The manuscript is interesting in that the authors go beyond the mere temporal characterization of neural activity, using advanced source reconstruction approach and linear mixed effect modeling for statistics. Based on their results, the authors report a posterior-anterior spatial gradient of the dominant peak frequency, and further our understanding of the link between structure and function.

Numbered summary of any substantive concerns

1) My main concern regards the authors' statement that there is a decrease of the intrinsic resonance frequency from early sensory to higher-order area.

i) In their analysis, the authors identify PF on the 1/f-free power spectra and then select the PF that has the highest power based on the original power spectra. Isn't this approach necessarily increasing the probability to detect a pic at low (around alpha) frequencies? In other words, such decrease in oscillations' frequency is actually a decrease of, specifically, alpha oscillations. The band-specific PF analysis actually suggests such conclusion.

ii) It has been shown before (e.g. Rosanova et al., 2009) that different lobes of the brain have different "intrinsic resonance frequencies," e.g. beta frequency appears in anterior regions, and the proposed analysis does not seem sensitive to it. The Band-specific PF analysis, however, shows that indeed beta increases from posterior to anterior regions.

iii) Based on the results reported in Figure 4A, the authors suggest that there is a change in PF and CT as a function of the region's function (i.e. sensory vs. associative regions). Yet, the data seem to suggest that the effect is driven by the location of the different regions: we see high PF in VIS and DAN, which are mainly located in the posterior part of the brain, and low PF in the other networks. And doesn't this interpretation fit better the fact that there is a link between CT and PF? In other words, change in CT are probably not sparse across the head as some of these networks seem to be, right?

2) Could the authors give an intuition of the variance in ROIs' size for a given participant and across participant, as well as how the size affects the power spectra?

3) The authors argue that there is a posterior-anterior gradient of CT, which strongly correlates with the PF gradient. From the results presented in Figure 3, it seems that the CT gradient is more strongly explained by the Z-axis (thicker in ventral than dorsal regions).

4) Could the authors give an intuition of how stable (within participant) these PF are over time (see for instance Haegens et al., 2014)?

5) Could the authors comment on the absence of gamma oscillations?

6) Could the authors comment on the generalizability of their results given the fact that they analyzed resting-state data?

Reviewer #2:

This study investigated a novel and interesting topic on whether there is anterior-posterior gradient in the frequency of brain oscillations using source-localized MEG data. It was found that the spatial gradient of strongest peak-frequency decreases gradually and robustly along the posterior-anterior axis following the hierarchy of early sensory and higher-order areas. The manuscript is themed around a robust spatial gradient of the frequency of brain oscillations while the main finding was the spatial gradient in the frequency of alpha oscillation but not across oscillatory frequency hierarchy. I would suggest reframing the text to be better in line the findings. Further, the main analysis are mostly sound but methodological details are missing and some parts are difficult to evaluate. Also the contribution of spatial leakage as confounding factor to the results has not been sufficiently well tested and should be better controlled for.

1) The strongest peak in the original power spectrum was selected. Taken the robust alpha oscillations in rest, these peaks were in the alpha-band range (8-12 Hz) as visible in Figure 1A. The results hence describe the spatial gradient of alpha oscillations in narrow frequency range of 7 to 12 Hz. Thus, I find the title of the manuscript "The frequency gradient of human resting-state brain oscillations follows cortical hierarchies" as well as structure of text somewhat misleading and suggestive of the presence of spatial oscillatory frequency hierarchy, which seems not to be the case. Also the spatial gradients of peaks in theta and beta frequencies were estimated, but these data are now only shown in the supplementary material, not thoroughly analysed and only discussed very superficially. It looks like these data were added only for the reviewers. I think the authors should decide do they want to present a frequency gradient or an alpha-frequency gradient and modify the manuscript accordingly.

2) The correlation of the frequency and cortical thickness was estimated to investigate if the frequency is correlated rather with cortical hierarchy than spatial location. A correlation of -0.14 was found (Figure 3B). Albeit very significant, the correlation itself is very weak as also shown by the large scattered peak frequencies. Compared to correlation of -0.84 of the peak frequency with posterior-anterior gradient, this result indicates that alpha peak frequency is not strongly related to cortical hierarchy. The Results and Discussion should be modified to emphasize this point. Moreover, the authors could perform a partial correlation analysis among the peak frequency, posterior-anterior direction and cortical thickness – does the correlation with thickness persist when the posterior-anterior gradient is factored out?

3) Related to the above, in Materials and methods the authors write that to obtain correlation between PF and CT scores, the robust correlation was performed (Pernet et al., 2012). Pernet et al. released a new open source Matlab toolbox for analysis of robust correlations using multiple statistical tests. Which of tests would the authors have used for the correlation analysis?

4) The observed spatial gradient of peak frequencies can reflect a true continuous frequency variation along an anatomical axis or, alternatively, it could be caused by source leakage of distributed alpha sources at different frequencies and different cortical areas. The authors argue that if the spatial gradient is due to leakage then there should not be frequency change in alpha oscillations in V1. As there was a significant correlation of PF of V1 and its nearby 0.5-1.5 cm sources, the authors concluded that there is no source-leakage from the frontal source that could explain the spatial gradient. However, the alpha sources that could become mixed and cause confounding the results do not need to be localized to V1 and frontal cortex and hence estimating the source-leakage of V1 alpha source appears inadequate in my opinion. Classical MEG alpha literature suggests that the strongest alpha sources could actually be located higher in the visual cortical hierarchy, in the parieto-occipital sulcus or in areas near precuneus which could well cause significant leakage affecting the analyses. Second, the spatial accuracy of source-localized MEG is > 1-2 cm and hence the area for which the source-leakage was estimated could have been too small to test whether there is leakage or not.

5) The authors have used a Beamforming approach for source reconstruction. Several details of the source reconstruction approach are missing or unclear. First, if Beamforming was used for the entire brain volume in voxel space, how was data transformed to cortical parcellations and how were the values for each parcel obtained? The sensor covariance matrix was computed for 2 sec trials. Which trials were these?

Reviewer #3:

Mahjoory et al. examined a large dataset of human MEG-recorded neuronal oscillations to probe the relation between oscillation properties across space and cortical thickness. Their main findings are that the brain shows anterior-posterior gradients in the frequency of neuronal oscillations and in the slope of 1/f nonoscillatory MEG patterns. They also find a matching A-P gradient in cortical thickness, and they suggest that cortical thickness correlates closely with oscillation frequency in individual subjects, even after accounting for mean anatomical variations in each pattern.

My enthusiasm for publishing this paper at eLife is limited because most of their findings are not especially novel. Specifically, it is already known that oscillation frequency varies across the brain (see work by Voytek et al., Groppe et al., Zhang et al., and others) and merely replicating this finding in a large open dataset is not extremely innovative. Similarly, Figure 5 on oscillatory peaks is not novel and neither is the result showing an overall A-P gradient in mean cortical thickness.

The paper's novel claim is showing that oscillation frequency correlates with cortical thickness even after accounting for mean anatomical patterns. However, this result is not adequate to justify publishing the entire paper, especially because the data underlying this result were not examined in much detail and there was no compelling and specific proposed mechanism to directly link these two phenomena. Further, I am concerned that this apparent correlation could actually be a reflection of intersubject differences in mean thickness and oscillation properties, rather than by a detailed region-by-region correspondence between these variables within-subject, as the authors suggest. Failing to rule out this possibility is a substantial weakness in this analysis and the authors could do more to demonstrate this effect at the within-subject level.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "The frequency gradient of human resting-state brain oscillations follows cortical hierarchies" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Laura Colgin as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Satu Palva (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.

Summary:

The authors have revised the manuscript thoroughly based on the previous comments. The manuscript is a very nice addition to the current literature showing how the frequency gradients of human resting-state neuronal oscillations are organized in accordance with cortical hierarchy. Revisions are still required to clarify in the main text, the Abstract and the figures that the patterns are present within subject, as well as within frequencies in contrast to across frequencies.

Revisions:

1) A key part of the paper is showing frequency gradients both within and across subjects. The rebuttal letter does a good job explaining that the authors' statistical framework identifies this pattern robustly both within and across subjects. But the text related to this is still unclear. The authors should revise the text of the results to more clearly explain how their statistical framework identifies gradients both within and across subjects.

2) The authors should also consider revising their figures to show clear examples of within and across subject gradients. The scatter plots and brain plots in Figure 1, for example, are hard to understand because they seem to combine data both across subjects and regions. It would be very informative if the figures followed the statistical results.

3) In response to one of the reviewers, the revised paper now includes an analysis of analyses within specific bands. This new analysis is a bit hard to follow in the context of the paper because it is unclear how it relates to the paper's primary analyses. Is the idea that there are multiple oscillatory patterns at different frequencies that all show gradients simultaneously? Additional clarity here would be helpful.

[Editors’ note: The authors appealed the original decision. What follows is the authors’ response to the first round of review.]

This manuscript presents results from a large dataset of human resting-state MEG recordings. The authors investigated the timely topic of the spatial distribution of low frequency brain oscillations, and found an anterior-posterior gradient in the frequency of brain oscillations. They further show a gradient in cortical thickness (CT), and suggest that CT correlates with oscillations frequency. Although investigating the spatial distribution of brain oscillations goes beyond simply studying their temporal dynamics, other studies have shown before that oscillation frequency varies across the brain.

This is correct and we acknowledge that in our manuscript. However, this is not the main result of our study. The novel contributions of our study are the following: (1) We provide the first comprehensive (in space and frequency) statistical model of frequency gradients in a large resting-state data set. The main result is not that different brain areas feature oscillations of different frequencies but that the frequency changes systematically and globally along spatial (and hierarchical) gradients. The only study that had previously tested this hypothesis is Zhang et al., 2018. However, their result was based on a smaller sample, data was recorded from epilepsy patients with inherent pathological brain activity and the ECoG data was constrained by limited electrode coverage. In addition, only anterior-posterior changes were studied. We provide the first full 3D statistical model at the level of individual brain areas (showing for example a significant gradient in inferior to superior direction). (2) We provide the first full 3D statistical model of cortical thickness in a large data set and show that cortical thickness changes systematically in space and is correlated with peak frequency. (3) We show for the first time that frequency gradients follow cortical thickness (as a proxy of hierarchical level) more closely than can be explained purely by spatial location.

Reviewers have also noticed that such spatial gradients are in fact more specific to the alpha frequency.

We agree with the reviewers in that our framework for PF identification is more specific to low frequencies. In our approach, we determine peaks from the 1/f-free spectrum and define the PF as the peak with the largest amplitude in the original spectrum, which leads to the preference of low frequencies. We also had tested another approach, where a PF was defined as a peak with the largest amplitude in the 1/f-free spectrum, thereby, the resulting PFs were in a wide range covering both alpha and beta rhythms. However, we found the PFs obtained from the latter approach extremely dependent on goodness of the 1/f fit, lacking a systematic pattern, and heavily inconsistent across participants. To address reviewers’ perfectly valid comment, we reframed the manuscript and included band-specific results in the main text as a complementary analysis to mirror exactly the analysis for low-frequency peaks. Accordingly, we modified the Results and Discussion sections. In the Results section; we now present the sub-section “Spatial Gradients are frequency specific” in more depth and with new plots.

Some methodological aspects have additionally been raised including addressing the question of spatial leakage, and the possible confound of inter individual differences.

On this point, we believe that we have addressed the reviewers’ comments by clarifying the points in more details and providing further evidences and analyses in the corresponding sections.

Reviewer #1:

Mahjoory et al. present new results on a publicly available, resting-state MEG recording dataset from a large pool of participants. This well-written manuscript provides a systematic analysis of the spatial distribution of high-power, low-frequency, oscillations, based on non-invasive recordings of healthy, human participants. The manuscript is interesting in that the authors go beyond the mere temporal characterization of neural activity, using advanced source reconstruction approach and linear mixed effect modeling for statistics. Based on their results, the authors report a posterior-anterior spatial gradient of the dominant peak frequency, and further our understanding of the link between structure and function.

Numbered summary of any substantive concerns

1) My main concern regards the authors' statement that there is a decrease of the intrinsic resonance frequency from early sensory to higher-order area.

i) In their analysis, the authors identify PF on the 1/f-free power spectra and then select the PF that has the highest power based on the original power spectra. Isn't this approach necessarily increasing the probability to detect a pic at low (around alpha) frequencies? In other words, such decrease in oscillations' frequency is actually a decrease of, specifically, alpha oscillations. The band-specific PF analysis actually suggests such conclusion.

Yes, we agree with the reviewer in that our results mainly shows the spatial gradients of low frequency (theta and alpha), which originates from the fact that we determine peaks from 1/f-free spectrum and define the PF as the peak with the largest amplitude in the original spectrum, which leads to the preference for low frequencies. In the revised manuscript, we present the results for both PF and band-specific PFs.

We have further analyzed (also in response to other comments) spatial gradients of the bandspecific PFs (theta, alpha, and beta) and their correlation with parameters of the 1/f component as well as cortical thickness. These analyses together with a new figure (Figure 6) are now included in the main text to show the results in more detail. Interestingly, in agreement with the main PF gradient, we found similar significant gradients (positive for theta and beta, negative for beta) of the band-specific PFs along the posterior-anterior axis. Overall these analyses support the notion that the posterior-anterior axis constitutes a principal dimension for frequency changes (both PF and band-specific PF) in the cortex.

ii) It has been shown before (e.g. Rosanova et al., 2009) that different lobes of the brain have different "intrinsic resonance frequencies," e.g. beta frequency appears in anterior regions, and the proposed analysis does not seem sensitive to it. The Band-specific PF analysis, however, shows that indeed beta increases from posterior to anterior regions.

Indeed, Rosanova (and others) have shown that different brain areas have different resonance frequencies. Similarly, in a previous study from our group we have demonstrated that different brain areas have characteristic spectral signatures (Keitel and Gross, 2016). Notably, the approach and aim of these previous studies is very different to ours. These previous studies focused on oscillatory power and treated each brain area separately, here, we focus on oscillatory frequency and study systematic changes in all spatial dimensions. So, we are indeed sensitive to different aspects of brain oscillations. This is also the case because Rosanova et al., 2009 used transcranial magnetic stimulation (TMS) to perturb directly different regions of the human cortex and then record the TMS-evoked brain response. We do not use a perturbation approach but rather study the resting-state brain activity.

We also agree with the reviewer that – irrespective of the different approaches – results of our band-specific PF analysis are consistent with findings from Rosanova and others. Also, we acknowledge that the term “intrinsic resonance frequency” can be confusing. To address the reviewer’s comment and avoid confusions for readers, we replaced “intrinsic resonance frequency” with “dominant frequency” in the revised manuscript.

iii) Based on the results reported in Figure 4A, the authors suggest that there is a change in PF and CT as a function of the region's function (i.e. sensory vs. associative regions). Yet, the data seem to suggest that the effect is driven by the location of the different regions: we see high PF in VIS and DAN, which are mainly located in the posterior part of the brain, and low PF in the other networks. And doesn't this interpretation fit better the fact that there is a link between CT and PF? In other words, change in CT are probably not sparse across the head as some of these networks seem to be, right?

We agree with the reviewer. There is a change of PF and CT associated with location and function (sensory versus associative). Both PF and CT show a posterior-anterior gradient. At the same time higher order associative brain areas are located more anterior compared to sensory brain areas (most prominently in the visual domain). Networks that are mostly localized in posterior brain areas (such as VIS) are expected to show different PF and CT compared to more anterior networks (such as DAN) merely by virtue of their location. Therefore, we agree with the reviewer that changes in PF and CT are not sparse as some of the networks are. Still, we believe that the “network representation” of our results is important to emphasize the point that these quantities (PF and CT) differ between networks. This is especially important for studies focusing a-priori on these anatomically defined networks and where these differences (in PF and CT) may confound the results.

However, we have also carried out further extensive analysis to test if PF and CT differ between networks irrespective of their location. This was done by removing from both measures changes that can be explained by linear dependencies on x,y,z. Interestingly, after regressing out the effect of spatial location, a significant difference between networks remained for PF (Figure 4—figure supplement 1A) and CT (Figure 4—figure supplement 1B). This was the case when looking at the standard resting-state networks and also when testing sensory against association areas. These results indicate that, beyond a global effect of location, networks still differ significantly in PF and CT after removing linear effects of location.

These results are now described in the main text and presented in detail in Figure 4— figure supplement 1 and 2.

2) Could the authors give an intuition of the variance in ROIs' size for a given participant and across participant, as well as how the size affects the power spectra?

In our analyses, we down-sampled the individual cortical surfaces to 8196 vertices and 16384 (triangular) faces and co-registered to a finer version of the Conte69 atlas. The atlas provides anatomical labels for each vertex but not for the faces. To obtain the surface area of a given cortical area, we first assigned a label to each face based on the label of the nearest vertex (Euclidean distance was computed between centroid of the given face and all vertices) and summed across the area values of the identically labeled faces. The obtained area values are a rough approximation but fine in this case because the spatial resolution of MEG of resting state reconstructed activity and the potential co-registration inaccuracy does not really warrant extremely fine-grained spatial details. Figure 1—figure supplement 2A shows the histogram of ROI’s size for a given participant. Figure 1—figure supplement 2B represents the mean and standard deviation of the ROI size across all participants. To assess the impact of ROI size on PF, we used again a linear mixed effect model (LMEM) where PF was defined as a function of the individual ROI area, defining the participant’s ID as a random effect to control for variability between individuals. Indeed, this model assesses the impact of ROIs’ size on PF gradient at the individual level and tests its consistency across participants. We found no significant impact of ROI size on PF (t = -1, p > 0.05). Next, we asked the question whether ROIs’ size could affect the spatial gradients of PF across the cortex irrespective of the PF itself. We used LMEM according to Equation 1 of the manuscript where we included the ROIs’ size as an additional fixed effect variable. We again found no significant impact of ROIs’ size on spatial gradients of PF along the Z and Y axes (t = -1.1, p > 0.05) (Figure 1—figure supplement 2C).

3) The authors argue that there is a posterior-anterior gradient of CT, which strongly correlates with the PF gradient. From the results presented in Figure 3, it seems that the CT gradient is more strongly explained by the Z-axis (thicker in ventral than dorsal regions).

Thanks for this constructive comment. Our focus so far has been on the posterior-anterior gradient. However, there is indeed strong evidence in our data for an inferior-superior gradient for both, CT and PF. Interestingly, following our preprint, another preprint has just confirmed these two gradients for CT (Valk et al., 2020). To address the reviewer's comment, we updated the manuscript and added z axis to the result and Discussion sections as follows:

“Notably, similar to the spatial gradients of PF, LMEM uncovered a secondary inferior-superior axis of CT. Having established the organizational axes of CT, we then tested for a significant relationship between CT and PF.”

4) Could the authors give an intuition of how stable (within participant) these PF are over time (see for instance Haegens et al, 2014)?

Thanks for raising this point. There is indeed increasing evidence that peak frequency is not constant but increases over tasks with increasing cognitive demand and decreases with increasing time on task (Benwell et al., 2019). Here, we did not expect to observe this type of change because participants did not engage in a task. But to fully address this question, we reconstructed a cortical activity time course similar to our main analysis approach. To explore consistency of our results across time, we halved source time courses of all individuals and ROIs, created two groups of data, and analyzed each group separately, as described in the Materials and methods section of the manuscript. Figure 1—figure supplement 1A illustrates the distribution of trimmed-mean-PF along the cortex for the first and second halves of the data. To statistically test the significance of the observed gradients, and their consistency across individuals, we applied LMEM on PF values of each half as a function of the ROIs’ centroid coordinates and their interactions (Figure 1—figure supplement 1B). To test whether the impact of time is significant or not, we defined the time as a categorical variable (“1^st half”, “2^nd half”) and added to our LMEM, where we defined PF as a function of centroid coordinates, their interactions and time. The statistical model showed a non-significant effect of time on our results (Figure 1—figure supplement 1C). We therefore conclude that PF and its gradient are stable over time.

5) Could the authors comment on the absence of gamma oscillations?

We initially had included gamma oscillations in our analysis. But after analyzing the data, we found gamma oscillations absent in most areas of the cortex across individuals, resulting in unreliable estimates of global gradients, which hindered us in deriving a conclusive result. This is likely at least partly due to the fact that the analysis of gamma oscillations would require a specifically optimized spectral analysis that is different from the analysis of oscillations with lower frequencies (e.g. the use of multi-taper spectral analysis might be required). Since we did not want to introduce a bias into our analysis (by having spectra with different properties), we focused on low-frequency oscillations.

However, we agree with the reviewer that changes of low frequencies across the cortex, may likely be coupled with occurrence of local gamma activities. For instance, (Voytek et al., 2010) have shown that high gamma amplitude is modulated by both the theta and alpha phase, and noted that there is a topographic separation of preferred coupling such that theta PAC is larger than alpha PAC over anterior brain regions, regardless of behavioral task, using ECOG recordings of 2 participants. We would like to keep these analyses for a future study where we account for cross-frequency couplings between low- and high-frequency oscillations and explore existence of systematic patterns across the cortex.

6) Could the authors comment on the generalizability of their results given the fact that they analyzed resting-state data?

Thank you for pointing this out. We are indeed interested in exploring implications of the observed frequency gradients in different cognitive states, disease, and aging. So far, studies have demonstrated that PF shifts as a function of cognitive load (Haegens et al., 2014), experiment duration (Benwell et al., 2019), and aging (Scally et al., 2018). But since we are the first to report a full 3-D statistical model of PF gradients it is currently unknown how they change as the brain starts to engage in a task. This will be the focus of a future study.

Reviewer #2:

This study investigated a novel and interesting topic on whether there is anterior-posterior gradient in the frequency of brain oscillations using source-localized MEG data. It was found that the spatial gradient of strongest peak-frequency decreases gradually and robustly along the posterior-anterior axis following the hierarchy of early sensory and higher-order areas. The manuscript is themed around a robust spatial gradient of the frequency of brain oscillations while the main finding was the spatial gradient in the frequency of alpha oscillation but not across oscillatory frequency hierarchy. I would suggest reframing the text to be better in line the findings. Further, the main analysis are mostly sound but methodological details are missing and some parts are difficult to evaluate. Also the contribution of spatial leakage as confounding factor to the results has not been sufficiently well tested and should be better controlled for.

We appreciate that the reviewer found our study ‘novel and interesting’. We have

performed further extensive analysis to address these important points in full. Analysis and results are described below in detail.

1) The strongest peak in the original power spectrum was selected. Taken the robust alpha oscillations in rest, these peaks were in the alpha-band range (8-12 Hz) as visible in Figure 1A. The results hence describe the spatial gradient of alpha oscillations in narrow frequency range of 7 to 12 Hz. Thus, I find the title of the manuscript "The frequency gradient of human resting-state brain oscillations follows cortical hierarchies" as well as structure of text somewhat misleading and suggestive of the presence of spatial oscillatory frequency hierarchy, which seems not to be the case. Also the spatial gradients of peaks in theta and beta frequencies were estimated, but these data are now only shown in the supplementary material, not thoroughly analysed and only discussed very superficially. It looks like these data were added only for the reviewers. I think the authors should decide do they want to present a frequency gradient or an alpha-frequency gradient and modify the manuscript accordingly.

We agree with the reviewer that spatial gradients are most pronounced at low frequencies. We now state that more clearly in the revised manuscript. But, as the reviewer correctly noted, we also show significant gradients in other frequency bands. In response to this comment we have now (1) included band-specific results (including results of new analysis) as a new subsection in the main text, (2) explained the details of analysis in the Materials and methods section, and (3) discussed the results in the Discussion section. In the Results section we now present the sub-section “Spatial Gradients are frequency specific” in the main text together with a new figure (Figure 6). Moreover, we thoroughly analyzed the spatial gradients of the band-specific PFs (Figure 6) and their correlation with parameters of the 1/f signal (Figure 6—figure supplement 1) and with cortical thickness (Figure 6—figure supplement 2). Interestingly, our results show that, just as for the PF, the spatial gradients of the band-specific PFs follow the posterior-anterior axis. However, this gradient is increasing from posterior to anterior for theta and beta range oscillations but decreasing for alpha band.

2) The correlation of the frequency and cortical thickness was estimated to investigate if the frequency is correlated rather with cortical hierarchy than spatial location. A correlation of -0.14 was found (Figure 3B). Albeit very significant, the correlation itself is very weak as also shown by the large scattered peak frequencies. Compared to correlation of -0.84 of the peak frequency with posterior-anterior gradient, this result indicates that alpha peak frequency is not strongly related to cortical hierarchy. The Results and Discussion should be modified to emphasize this point. Moreover, the authors could perform a partial correlation analysis among the peak frequency, posterior-anterior direction and cortical thickness – does the correlation with thickness persist when the posterior-anterior gradient is factored out?

The low correlation actually reflects the fact that correlation analysis is not optimal in this case. It does not take into consideration the variability between participants (e.g. the fact that some participants have an overall higher occipital alpha frequency compared to others). This is also one of the reasons for the large scattering of values in Figure 3—figure supplement 2. Therefore, LMEM is, again, the preferred statistical tool because it models specifically the interindividual differences and, as a result, leads to more robust and highly significant results. Since both are a bit redundant and LMEM is more appropriate we now moved the correlation analysis to a supplement figure (Figure 3—figure supplement 2). We still included Figures 3—figure supplement 2, and the correlation results because they are intuitive and easily accessible to the readers.

We had already addressed the reviewer’s second point. Factoring out the dependency on location (x,y,z) from PF and CT still preserves the significant association between cortical thickness and peak frequency (Figures 3—figure supplement 1). This is indeed an important result because it suggests that PF and CT are related in a way that cannot be explained by a common linear dependence on location.

3) Related to the above, in Materials and methods the authors write that to obtain correlation between PF and CT scores, the robust correlation was performed (Pernet et al., 2012). Pernet et al. released a new open source Matlab toolbox for analysis of robust correlations using multiple statistical tests. Which of tests would the authors have used for the correlation analysis?

After plotting the data (scatter plot) and inspecting all correlation alternatives offered by Robust Correlation Toolbox, we chose Skipped Pearson Correlation for our analysis because it captured most of the variation in the data. However, all alternative techniques also led to similar results. To address this comment, we updated the manuscript and added further details for the applied correlation analysis.

Materials and methods section:

“To obtain the correlation between PF and CT scores, we initially computed the 10% trimmed mean across participants for each ROI and performed a robust correlation (skipped Pearson correlation, Robust Correlation Toolbox) (Pernet et al., 2012)”

4) The observed spatial gradient of peak frequencies can reflect a true continuous frequency variation along an anatomical axis or, alternatively, it could be caused by source leakage of distributed alpha sources at different frequencies and different cortical areas. The authors argue that if the spatial gradient is due to leakage then there should not be frequency change in alpha oscillations in V1. As there was a significant correlation of PF of V1 and its nearby 0.5-1.5 cm sources, the authors concluded that there is no source-leakage from the frontal source that could explain the spatial gradient. However, the alpha sources that could become mixed and cause confounding the results do not need to be localized to V1 and frontal cortex and hence estimating the source-leakage of V1 alpha source appears inadequate in my opinion. Classical MEG alpha literature suggests that the strongest alpha sources could actually be located higher in the visual cortical hierarchy, in the parieto-occipital sulcus or in areas near precuneus which could well cause significant leakage affecting the analyses. Second, the spatial accuracy of source-localized MEG is > 1-2 cm and hence the area for which the source-leakage was estimated could have been too small to test whether there is leakage or not.

We agree with the reviewer in that alpha generators could be distributed at different frequencies and different cortical areas. However, we believe that our analysis is relatively insensitive to the exact reference location for the following reason: The simplest scenario for a gradient caused by leakage consists of two sources of different frequencies (for example an occipito-parietal alpha source and a frontal theta source). Differential leakage of occipital and frontal sources at different points on the posterior-anterior line could result in a gradient. To address this point, according to the reviewer’s valid comment, we computed the geodesic distance between the reference ROI, V1, and all areas located 2–3cm away from V1, and applied linear mixed effect modelling of PF as a function of the distance values. We found a highly significant negative dependence between PF and distance (t = -21.11, p << 0.001).

To further account for the potential confounding effect of spatial leakage we performed a new comprehensive analysis, where we computed the geodesic distance between centroid of all ROIs and the centroid of V1 and used this as a new axis because it well explains the posterior-anterior axis in the cortex (Figure 1D bottom panel). We applied an LMEM between the PF as a response variable and the geodesic distance values as an independent variable to assess the distribution of PF along the distance values while considering inter-individual variability. We found a highly significant gradient of PF along the specified geodesic distance (t = -18.9, p << 0.001). Furthermore, to answer the question whether the spatial gradient of PF constantly exists in different areas of the cortex, we split the cortex to 3 equal, consecutive and non-overlapping windows (about 4 cm) based on the Y axis, applied LMEM for each window modelling PF as a function of geodesic distance, and found a significant gradient (window 1: t = -8.1, window 2: t = -10.9, window 3: t = -5.2, all p < 0.001) (Figure 1E). Indeed, our analysis demonstrates a significant organization of PF along the posterior-anterior direction for all windows indicating that this axis constitutes a systematic and constant gradient of PF. We updated the manuscript as follows:

“To address this question, we computed the geodesic distance between the reference ROI, V1, and all areas located 2–3 cm away from V1, and applied linear mixed effect modelling of PF as a function of the distance values. […] Indeed, our analysis demonstrates a significant organization of PF along the posterior-anterior direction for all windows indicating that this axis constitutes a systematic and constant gradient of PF.”

5) The authors have used a Beamforming approach for source reconstruction. Several details of the source reconstruction approach are missing or unclear. First, if Beamforming was used for the entire brain volume in voxel space, how was data transformed to cortical parcellations and how were the values for each parcel obtained? The sensor covariance matrix was computed for 2 sec trials. Which trials were these?

In the revised manuscript we have added further details. Beamforming was used not in the entire brain volume but directly on vertices of the cortical surface. For each anatomical ROI we performed dimensionality reduction using singular value decomposition (svd) on all vertex timeseries. We retained the required number of components to account for 95% of the variance for each ROI (typically 15 components). The final spectrum was computed as the 10% trimmed mean across spectra. The 2-sec trials refer to data segments obtained from cutting the continuous data into 2-sec long segments. These details have now been further clarified in the revised manuscript. The relevant parts read:

“[…] Using the Caret software, the mid-thickness cortical surface (halfway between the pial and white matter surfaces) was generated. The cortex surface was parceled into 384 ROIs (192 per hemisphere) according to Schoffelen et al. (Schoffelen et al., 2017). 648 vertices located in the medial wall (sub-cortical areas) were excluded from further analysis. […] The sensor covariance matrix was computed between all MEG-sensor pairs, as the average covariance matrix after cutting the continuous data into 2-second data segments. […] For each anatomical ROI we performed dimensionality reduction using singular value decomposition (SVD) on all vertex timeseries. We retained the required number of components to account for 95% of the variance for each ROI (typically 15 components). Component time courses were segmented to twosecond epochs, from which power spectra were computed using a multi-tapered Fast Fourier transform, using discrete prolate spheroidal sequences (dpss) as windowing function, with 2 Hz spectral smoothing. To obtain a single spectrum for each ROI (ROI spectrum), we pooled spectra of epochs across components and computed the 10% trimmed mean across them. Averaging after leaving out 10% of data from left and right tails of the spectra distribution offers a more robust estimate.”

Reviewer #3:

Mahjoory et al. examined a large dataset of human MEG-recorded neuronal oscillations to probe the relation between oscillation properties across space and cortical thickness. Their main findings are that the brain shows anterior-posterior gradients in the frequency of neuronal oscillations and in the slope of 1/f nonoscillatory MEG patterns. They also find a matching A-P gradient in cortical thickness, and they suggest that cortical thickness correlates closely with oscillation frequency in individual subjects, even after accounting for mean anatomical variations in each pattern.

My enthusiasm for publishing this paper at eLife is limited because most of their findings are not especially novel. Specifically, it is already known that oscillation frequency varies across the brain (see work by Voytek et al., Groppe et al., Zhang et al., and others) and merely replicating this finding in a large open dataset is not extremely innovative. Similarly, Figure 5 on oscillatory peaks is not novel and neither is the result showing an overall A-P gradient in mean cortical thickness.

The reviewers state that it is already known that the frequency of oscillations changes across the brain. This is correct and we acknowledge that in our manuscript. However, this is not the main result of our study. The novel contributions of our study are the following:

(1) We provide the first comprehensive (in space and frequency) statistical model of frequency gradients in a large resting-state data set. The main result is not that frequency changes across the cortex but that it changes systematically and globally along spatial (and hierarchical) gradients. To the best of our knowledge, the only study that had previously tested this hypothesis is Zhang et al., 2018. However, their result was based on a smaller sample, data was recorded from epilepsy patients with inherent pathological brain activity and the ECoG data was constrained by limited electrode coverage. In addition, only anterior-posterior changes were studied. We provide the first full 3D statistical model at the level of individual brain areas (showing for example a significant gradient in inferior to superior direction).

(2) We provide the first full 3D statistical model of cortical thickness in a large data set and show that cortical thickness changes systematically in space and is correlated with peak frequency.

(3) We show for the first time that frequency gradients follow cortical thickness (as a proxy of hierarchical level) more closely than can be explained purely by spatial location.

We believe that the relevance of these three novel findings based on state-of-the-art statistical models in a large data set meet the criteria for this prestigious journal.

The paper's novel claim is showing that oscillation frequency correlates with cortical thickness even after accounting for mean anatomical patterns. However, this result is not adequate to justify publishing the entire paper, especially because the data underlying this result were not examined in much detail and there was no compelling and specific proposed mechanism to directly link these two phenomena. Further, I am concerned that this apparent correlation could actually be a reflection of intersubject differences in mean thickness and oscillation properties, rather than by a detailed region-by-region correspondence between these variables within-subject, as the authors suggest. Failing to rule out this possibility is a substantial weakness in this analysis and the authors could do more to demonstrate this effect at the within-subject level.

The correlation of oscillations frequency and cortical thickness is not the only novelty but only one of three main novel results (as listed in our response to the previous comment). In addition, we can rule out the main concern of the reviewer that the correlation reflects intersubject variability. This is an important concern but it is already addressed in our analysis. First, we made sure that our analysis does not introduce unnecessary variability across participants. Our source analysis and subsequent parcellation was based on a standard anatomical atlas and the procedure ensured that the number of vertices for a given labelled area is the same across participants. Second, and more importantly, our use of linear-mixed-effect models (LMEM) ensured that individual differences are properly accounted for and that significant gradients are consistently present in individual participants. For example, participants have different alpha peak frequency in occipital brain areas and the slope of the frequency gradient is different. These individual differences are specifically modelled by LMEM as random effects. Importantly, LMEM applies two-level statistics, and therefore, will only show a significant gradient if it is significant across cortical areas at the individual level, and consistent across participants at the group level. We have further corroborated the existence of spatial gradients by computing linear correlations for each individual and then analyzing individual correlations across the group. At the individual level, we computed the Spearman correlation between Y coordinates and PF values. For group level analysis we first computed the inverse hyperbolic tangent of the obtained correlation values, and applied one sample t-test across them. This approach tests the hypothesis that, across the group, correlations deviate significantly from zero. The results showed a highly significant negative correlation (t-value = -15.52, p < 0.001) between PF and Y coordinates. As expected, the results are consistent with the LMEM results. Both approaches show significant gradients at the individual level. However, we use LMEM throughout the manuscript because it is the statistically superior and more versatile approach.

[Editors’ note: what follows is the authors’ response to the second round of review.]

Revisions:

1) A key part of the paper is showing frequency gradients both within and across subjects. The rebuttal letter does a good job explaining that the authors' statistical framework identifies this pattern robustly both within and across subjects. But the text related to this is still unclear. The authors should revise the text of the results to more clearly explain how their statistical framework identifies gradients both within and across subjects.

To address the mentioned points we updated the following sections of the text:

• We updated the first subsection of the result section, adding more details about the significance of the posterior-to-anterior gradient of PF at individual level across ROIs and its consistency across individuals.

• In the result section, we now describe further details about the LMEM and how it deals with the intra- and inter-individual variability.

• We added new plots to Figure 1, which illustrate the correlation between PF and Y-coordinates at the individual level and their distribution across subjects.

• We updated the Discussion section further discussing the inter-individual variability.

2) The authors should also consider revising their figures to show clear examples of within and across subject gradients. The scatter plots and brain plots in Figure 1, for example, are hard to understand because they seem to combine data both across subjects and regions. It would be very informative if the figures followed the statistical results.

Thanks for suggesting this important point, which indeed requires further clarification. To illustrate within- and across-individual gradient we added new plots to Figures 1B and 1D, top-right panel, to display the individual correlation values across ROIs and to show the distribution of correlation values across all participants.

3) In response to one of the reviewers, the revised paper now includes an analysis of analyses within specific bands. This new analysis is a bit hard to follow in the context of the paper because it is unclear how it relates to the paper's primary analyses. Is the idea that there are multiple oscillatory patterns at different frequencies that all show gradients simultaneously? Additional clarity here would be helpful.

In our analysis pipeline, we obtained power spectrum for each region of interest (ROI), from which we determined all spectral peaks at the range of 3–45 Hz. This analysis led to identification of multiple peaks for (some) ROIs. Figure 5 represents a histogram plot of all identified peaks pooled across all ROIs and individuals. The plot nicely delineates canonical oscillatory frequency bands. We then extracted peaks, separately, within each band for each participant and ROI and repeated the LMEM analysis. Indeed, this revealed existence of band-specific gradients in addition to the dominant peak frequency gradient. This means that peak frequencies change consistently across space in these canonical frequency bands. We have now clarified this point in the Results and Discussion section. In addition, we now describe more clearly the term “peak frequency” which also helps to clarify the distinction between gradient of peak frequency and gradients of band-specific peaks.

References

Benwell, C.S.Y., London, R.E., Tagliabue, C.F., Veniero, D., Gross, J., Keitel, C., Thut, G., 2019. Frequency and power of human alpha oscillations drift systematically with time-on-task. Neuroimage 192, 101–114. doi:10.1016/j.neuroimage.2019.02.067

Haegens, S., Cousijn, H., Wallis, G., Harrison, P.J., Nobre, A.C., 2014. Inter- and intra-individual variability in alpha peak frequency. Neuroimage 92, 46–55. doi:10.1016/j.neuroimage.2014.01.049

Scally, B., Burke, M.R., Bunce, D., Delvenne, J.-F., 2018. Resting-state EEG power and connectivity are associated with alpha peak frequency slowing in healthy aging. Neurobiol. Aging 71, 149–155. doi:10.1016/j.neurobiolaging.2018.07.004

Author details

1. Keyvan Mahjoory

   Institute for Biomagnetism and Biosignalanalysis (IBB), University of Muenster, Muenster, Germany

   Contribution

   Resources, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

   For correspondence

   kmahjoory@gmail.com

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-0386-1135

2. Jan-Mathijs Schoffelen

   Radboud University Nijmegen, Donders Institute for Brain, Cognition and Behaviour, Nijmegen, Netherlands

   Contribution

   Data curation, Software, Formal analysis, Methodology, Writing - review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-0923-6610

3. Anne Keitel

   Psychology, University of Dundee, Scrymgeour Building, Dundee, United Kingdom

   Contribution

   Writing - review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-4498-0146

4. Joachim Gross

   1. Institute for Biomagnetism and Biosignalanalysis (IBB), University of Muenster, Muenster, Germany
   2. Centre for Cognitive Neuroimaging (CCNi), University of Glasgow, Glasgow, United Kingdom
   3. Otto-Creutzfeldt-Center for Cognitive and Behavioral Neuroscience, University of Muenster, Muenster, Germany

   Contribution

   Supervision, Funding acquisition, Validation, Methodology, Writing - review and editing

   For correspondence

   joachim.gross@uni-muenster.de

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-3994-1006

• 5,535

  Page views

• 698

  Downloads

• 38

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.