# Neuronal timescales are functionally dynamic and shaped by cortical microarchitecture

1. RichardGao r.dg.gao@gmail.com 1
2. RuudLvandenBrink 2
3. ThomasPfeffer 3
4. BradleyVoytek 1 4 5 6

1. Department of Cognitive Science, University of California, San Diego
2. Section Computational Cognitive Neuroscience, Department of Neurophysiology and Pathophysiology, University Medical Center Hamburg-Eppendorf
3. Center for Brain and Cognition, Computational Neuroscience Group, Universitat Pompeu Fabra
4. Halıcıoğlu Data Science Institute, University of California, San Diego
5. Neurosciences Graduate Program, University of California, San Diego
6. Kavli Institute for Brain and Mind, University of California, San Diego

2020-11-23

- Research Article

- Computational and Systems Biology
- Neuroscience

- neuronal timescales
- cortical gradients
- functional specialization
- transcriptomics
- spectral analysis
- Human
- Rhesus macaque

- publisher-id61277
- doi10.7554/eLife.61277
- elocation-ide61277

## Abstract

Complex cognitive functions such as working memory and decision-making require information maintenance over seconds to years, from transient sensory stimuli to long-term contextual cues. While theoretical accounts predict the emergence of a corresponding hierarchy of neuronal timescales, direct electrophysiological evidence across the human cortex is lacking. Here, we infer neuronal timescales from invasive intracranial recordings. Timescales increase along the principal sensorimotor-to-association axis across the entire human cortex, and scale with single-unit timescales within macaques. Cortex-wide transcriptomic analysis shows direct alignment between timescales and expression of excitation- and inhibition-related genes, as well as genes specific to voltage-gated transmembrane ion transporters. Finally, neuronal timescales are functionally dynamic: prefrontal cortex timescales expand during working memory maintenance and predict individual performance, while cortex-wide timescales compress with aging. Thus, neuronal timescales follow cytoarchitectonic gradients across the human cortex and are relevant for cognition in both short and long terms, bridging microcircuit physiology with macroscale dynamics and behavior.

    # set up behind the scenes
    %matplotlib inline
    import matplotlib.pyplot as plt
    from matplotlib.gridspec import GridSpec
    plt.style.use('./matplotlibrc_notebook')
    import numpy as np
    from scipy import stats
    import pandas as pd
    from fooof import FOOOF, FOOOFGroup
    from neurodsp import sim, spectral
    from statsmodels.tsa.stattools import acf
    import altair as alt

    # For compatability with Stencila output Altair plots
    # using MIME type renderer
    alt.renderers.enable('mimetype')

    import warnings
    warnings.filterwarnings('ignore')
    C_ORD = plt.rcParams['axes.prop_cycle'].by_key()['color']

    import nibabel as ni
    from surfer import Brain

    def compute_perm_corr(x, y, y_nulls, corr_method='spearman'):
        corr_func = stats.spearmanr if corr_method is 'spearman' else stats.pearsonr
        rho, pv = corr_func(x,y)
        rho_null = np.array([corr_func(x, n_)[0] for n_ in y_nulls])
        pv_perm = (abs(rho)<abs(rho_null)).sum()/y_nulls.shape[0]
        return rho, pv, pv_perm, rho_null

    def convert_knee_val(knee, exponent=2.):
        """
        Convert knee parameter to frequency and time-constant value.
        Can operate on array or float.

        Default exponent value of 2 means take the square-root, but simulation shows
        taking the exp-th root returns a more accurate drop-off frequency estimate
        when the PSD is actually Lorentzian.
        """
        knee_freq = knee**(1./exponent)
        knee_tau = 1./(2*np.pi*knee_freq)
        return knee_freq, knee_tau

    def sig_str(rho, pv, pv_thres=[0.05, 0.01, 0.005, 0.001], form='*', corr_letter=r'$\rho$'):
        """Generates the string to print rho and p-value.

        Parameters
        ----------
        rho : float
        pv : float
        pv_thres : list
            P-value thresholds to for successive # of stars to print.
        form : str
            '*' to print stars after rho, otherwise print p-value on separate line.

        Returns
        -------
        str
        """
        if form == '*':
            s = corr_letter+' = %.2f '%rho + np.sum(pv<=np.array(pv_thres))*'*'
        else:
            if pv<pv_thres[-1]:
                s = corr_letter+' = %.2f'%rho+ '\np < %.3f'%pv_thres[-1]
            else:
                s = corr_letter+' = %.2f'%rho+ '\np = %.3f'%pv
        return s

    # load data and helper funcs
    # helper function for plotting the parcellation results
    def plot_MMP(data, minmax=None, thresh=None, cmap='inferno', alpha=1, add_border=False, bp=0, title=None):
        """
        Plots arbitrary array of data onto MMP parcellation
        """
        from mpl_toolkits.axes_grid1.inset_locator import inset_axes
        ## Following this tutorial
        # https://github.com/nipy/PySurfer/blob/master/examples/plot_parc_values.py
        mmp_labels, ctab, names = ni.freesurfer.read_annot('./data/lh.HCPMMP1.annot')

        if minmax is None:
            minmax = [np.min(data), np.max(data)]

        if len(names)-len(data)==1:
            # short one label cuz 0 is unassigned in the data, fill with big negative number
            data_app = np.hstack((minmax[0],data))
            vtx_data = data_app[mmp_labels]
        else:
            vtx_data = data[mmp_labels]
            vtx_data[mmp_labels < 1] = minmax[0]

        # plot brain
        brain = Brain('fsaverage', 'lh', 'inflated', subjects_dir='./data/',
                      cortex=None, alpha=1, background='white', size=800, show_toolbar=False, offscreen=True)

        if add_border:
            brain.add_annotation((mmp_labels, np.array([[0,0,0, c[3], c[4]] for c in ctab])))

        # add data
        if thresh is None:
            thresh = minmax[0] - abs(minmax[0]*0.05)
        brain.add_data(vtx_data, minmax[0], minmax[1], colormap=cmap, alpha=alpha, colorbar=False, thresh=thresh)

        # plot brain views
        brainviews = brain.save_imageset(None, ['lat', 'med'])

        # merge brainviews and plot horizontally
        plt.imshow(np.concatenate(brainviews, axis=1), cmap=cmap)
        plt.gca().spines['bottom'].set_visible(False)
        plt.gca().spines['left'].set_visible(False)

        plt.xticks([])
        plt.yticks([])
        cbaxes = inset_axes(plt.gca(), width="50%", height="4%", loc=8, borderpad=bp)
        plt.colorbar(cax=cbaxes, orientation='horizontal')
        plt.clim(minmax[0], minmax[1])
        if title is not None:
            plt.title(title)

    <>:27: SyntaxWarning: "is" with a literal. Did you mean "=="?
    <>:27: SyntaxWarning: "is" with a literal. Did you mean "=="?

## Introduction

Human brain regions are broadly specialized for different aspects of behavior and cognition, and the temporal dynamics of neuronal populations across the cortex are thought to be an intrinsic property (i.e., neuronal timescale) that enables the representation of information over multiple durations in a hierarchically embedded environment 57Kiebel et al.2008. For example, primary sensory neurons are tightly coupled to changes in the environment, firing rapidly to the onset and removal of a stimulus, and showing characteristically short intrinsic timescales 70Ogawa and Komatsu2010 76Runyan et al.2017. In contrast, neurons in cortical association (or transmodal) regions, such as the prefrontal cortex (PFC), can sustain their activity for many seconds when a person is engaged in working memory 105Zylberberg and Strowbridge2017, decision-making 40Gold and Shadlen2007, and hierarchical reasoning 77Sarafyazd and Jazayeri2019. This persistent activity in the absence of immediate sensory stimuli reflects longer neuronal timescales, which is thought to result from neural attractor states 94Wang2002 102Wimmer et al.2014 shaped by N-methyl-D-aspartate receptor (NMDA)-mediated recurrent excitation and fast feedback inhibition 95Wang2008 93Wang1999, with contributions from other synaptic and cell-intrinsic properties 23Duarte and Morrison2019 37Gjorgjieva et al.2016. How connectivity and various cellular properties combine to shape neuronal dynamics across the cortex remains an open question.

Anatomical connectivity measures based on tract tracing data, such as laminar feedforward vs. feedback projection patterns, have classically defined a hierarchical organization of the cortex 26Felleman and Van Essen1991 46Hilgetag and Goulas2020 85Vezoli et al.2020. Recent studies have also shown that variations in many microarchitectural features follow continuous and coinciding gradients along a sensory-to-association axis across the cortex, including cortical thickness, cell density, and distribution of excitatory and inhibitory neurons 49Huntenburg et al.2018 98Wang2020. In particular, gray matter myelination 39Glasser and Van Essen2011—a noninvasive proxy of anatomical hierarchy consistent with laminar projection data—varies with the expression of genes related to microcircuit function in the human brain, such as NMDA receptor and inhibitory cell-type marker genes 10Burt et al.2018. Functionally, specialization of the human cortex, as well as structural and functional connectivity 64Margulies et al.2016, also follow similar macroscopic gradients. Moreover, in addition to the broad differentiation between sensory and association cortices, there is evidence for an even finer hierarchical organization within the frontal cortex 77Sarafyazd and Jazayeri2019. For example, the anterior-most parts of the PFC are responsible for long timescale goal-planning behavior 2Badre and D'Esposito2009 88Voytek et al.2015, while healthy aging is associated with a shift in these gradients such that older adults become more reliant on higher-level association regions to compensate for altered lower-level cortical functioning 17Davis et al.2008.

Despite convergent observations of cortical gradients in structural features and cognitive specialization, there is no direct evidence for a similar gradient of neuronal timescales across the human cortex. Such a gradient of neuronal dynamics is predicted to be a natural consequence of macroscopic variations in synaptic connectivity and microarchitectural features 13Chaudhuri et al.2015 22Duarte et al.2017 48Huang and Doiron2017 49Huntenburg et al.2018 98Wang2020, and would be a primary candidate for how functional specialization emerges as a result of hierarchical temporal processing 57Kiebel et al.2008. Single-unit recordings in rodents and non-human primates demonstrated a hierarchy of timescales that increase, or lengthen, progressively along a posterior-to-anterior axis 21Dotson et al.2018 68Murray et al.2014 76Runyan et al.2017 99Wasmuht et al.2018, while intracranial recordings and functional neuroimaging data collected during perceptual and cognitive tasks suggest likewise in humans 3Baldassano et al.2017 47Honey et al.2012 61Lerner et al.2011 100Watanabe et al.2019. However, these data are either sparsely sampled across the cortex or do not measure neuronal activity at the cellular and synaptic level directly, prohibiting the full construction of an electrophysiological timescale gradient across the human cortex. As a result, while whole-cortex data of transcriptomic and anatomical variations exist, we cannot take advantage of them to dissect the contributions of synaptic, cellular, and circuit connectivity in shaping fast neuronal timescales, nor ask whether regional timescales are dynamic and relevant for human cognition.

Here we combine several publicly available datasets to infer neuronal timescales from invasive human electrocorticography (ECoG) recordings and relate them to whole-cortex transcriptomic and anatomical data, as well as probe their functional relevance during behavior (Figure 1A for schematic of study; Tables 1 and 2 for dataset information). Unless otherwise specified, (neuronal) timescale in the following sections refers to ECoG-derived timescales, which are more reflective of fast synaptic and transmembrane current timescales than single-unit or population spiking timescales (Figure 1A, left box), though we demonstrate in macaques a close correspondence between the two. In humans, neuronal timescales increase along the principal sensorimotor-to-association axis across the cortex and align with macroscopic gradients of gray matter myelination (T1w/T2w ratio) and synaptic receptor and ion channel gene expression. Finally, we find that human PFC timescales expand during working memory maintenance and predict individual performance, while cortex-wide timescales compress with aging. Thus, neuronal timescales follow cytoarchitectonic gradients across the human cortex and are relevant for cognition in both short and long terms, bridging microcircuit physiology with macroscale dynamics and behavior.

Table 1.

### Summary of open-access datasets used.

Data | Ref. | Specific source/format used | Participant info | Relevant figures
--- | --- | --- | --- | ---
MNI Open iEEG Atlas | 29Frauscher et al.2018; 30Frauscher et al.2018 |  | N = 105 (48 females) Ages: 13–65, 33.4 ± 10.6 | Figure 2A–D, Figure 3, Figure 4E and F
T1w/T2w and cortical thickness maps from Human Connectome Project | 38Glasser et al.2016; 39Glasser and Van Essen2011 | Release S1200, March 1, 2017 | N = 1096 (596 females) Age: 22–36+ (details restricted due to identifiability) | Figure 2C and D, Figure 3D–F
Neurotycho macaque ECoG | 69Nagasaka et al.2011; 103Yanagawa et al.2013 | Eyes-open state from anesthesia datasets (propofol and ketamine) | Two animals (Chibi and George) four sessions each | Figure 2E–G
Macaque single-unit timescales | 68Murray et al.2014 | Figure 1 of reference |  | Figure 2E–G
Whole-cortex interpolated Allen Brain Atlas human gene expression | 42Gryglewski et al.2018; 43Hawrylycz et al.2012 | Interpolated maps downloadable from http://www.meduniwien.ac.at/neuroimaging/mRNA.html | N = 6 (one female) Age: 24, 31, 39, 49, 55, 57 (42.5 ± 12.2) | Figure 3
Single-cell timescale-related genes | 5Bomkamp et al.2019; 81Tripathy et al.2017 | Table S3 from 81Tripathy et al.2017, Online Table 1 from 5Bomkamp et al.2019 | N = 170 (81Tripathy et al.2017) and 4168 (5Bomkamp et al.2019) genes | Figure 3C and D
Human working memory ECoG | 54Johnson2019; 53Johnson2018; 51Johnson et al.2018, 52Johnson et al.2018 | CRCNS fcx-2 and fcx-3 | N = 14 (five females) Age: 22–50, 30.9 ± 7.8 | Figure 4A–D

Table 2.

### Reproducing figures from code repository.

All IPython notebooks 35Gao2020: https://github.com/rdgao/field-echos/tree/master/notebooks

Notebook | Results
--- | ---
1_sim_method_schematic.ipynb | Simulations: Figure 1B–E
2_viz_NeuroTycho-SU.ipynb | Macaque timescales: Figure 2E–G, Figure 2—figure supplement 4
3_viz_human_structural.ipynb | Human timescales vs. T1w/T2w and gene expression: Figure 2A–D, Figure 2—figure supplements 1 and 3, Figure 3, Figure 3—figure supplements 1 and 2, Supplementary file 1–, Supplementary file 2, Supplementary file 3.
4b_viz_human_wm.ipynb | Human working memory: Figure 4A–D, Figure 4—figure supplement 1
4a_viz_human_aging.ipynb | Human aging: Figure 4E and F, Figure 4—figure supplement 2
supp_spatialautocorr.ipynb | Spatial autocorrelation-preserving nulls:
supp_spatialautocorr.ipynb | Spatial autocorrelation-preserving nulls: Figure 2—figure supplement 2

Figure 1A.

### Schematic of study and timescale inference technique.

(A) In this study, we infer neuronal timescales from intracranial field potential recordings, which reflect integrated synaptic and transmembrane current fluctuations over large neural populations 12Buzsáki et al.2012. Combining multiple open-access datasets (Table 1), we link timescales to known human anatomical hierarchy, dissect its cellular and physiological basis via transcriptomic analysis, and demonstrate its functional modulation during behavior and through aging.

Figure 1.

### Schematic of study and timescale inference technique.

(A) In this study, we infer neuronal timescales from intracranial field potential recordings, which reflect integrated synaptic and transmembrane current fluctuations over large neural populations 12Buzsáki et al.2012. Combining multiple open-access datasets (Table 1), we link timescales to known human anatomical hierarchy, dissect its cellular and physiological basis via transcriptomic analysis, and demonstrate its functional modulation during behavior and through aging. (B) Simulated time series and their (C) autocorrelation functions (ACFs), with increasing (longer) decay time constant, τ (which neuronal timescale is defined to be). (D) Example human electrocorticography (ECoG) power spectral density (PSD) showing the aperiodic component fit (red dashed), and the ‘knee frequency’ at which power drops off (f k, red circle; insets: time series and ACF). (E) Estimation of timescale from PSDs of simulated time series in (B), where the knee frequency, f k, is converted to timescale, τ, via the embedded equation (inset: correlation between ground truth and estimated timescale values).

    # simulate noise
    T = 240
    fs = 2000.
    t_ds = np.arange(0.005,0.08,0.01)
    f_to_plot=100
    noise, ac = [], []
    for t_d in t_ds:
        noise.append(sim.sim_synaptic_current(T, fs, tau_d = t_d))
        ac.append(acf(noise[-1], nlags=int(fs), fft=True))

    noise = np.vstack(noise)
    ac = np.vstack(ac).T
    f_axis, PSD = spectral.compute_spectrum(noise,fs)

    # FOOOF PSDs without knee
    fg = FOOOFGroup(aperiodic_mode='knee', max_n_peaks=0, verbose=False)
    fg.fit(freqs=f_axis, power_spectra=PSD, freq_range=(2,200))
    fit_knee = fg.get_params('aperiodic_params', 'knee')
    fit_exp = fg.get_params('aperiodic_params', 'exponent')
    knee_freq, taus = convert_knee_val(fit_knee, fit_exp)
    P_knee = [PSD[i,np.argmin(np.abs((f_axis[:f_to_plot]-(knee_freq[i]))))] for i in range(len(t_ds))]

    ### simulation data and analysis ###
    color = plt.cm.inferno(np.linspace(0,1,len(t_ds)))
    plt.rcParams['axes.prop_cycle'] = plt.cycler('color', color)
    i_plot = [0,3,7]
    t = np.arange(0,T,1/fs)
    plt.figure(figsize=(12,3))
    for p in range(3):
        plt.subplot(3,1,p+1)
        plt.plot(t[:int(fs)], noise.T[:int(fs),i_plot[p]], alpha=0.9, color=np.array(plt.cycler('color', plt.cm.inferno(np.linspace(0,0.85,3))))[p]['color'])
        plt.xticks([]); plt.yticks([])
        plt.xlim([0,1]);

    plt.xlabel('time (s)', labelpad=-10)
    plt.xticks([0,1], fontsize=15);
    plt.subplot(3,1,1)
    plt.title('simulated data')
    plt.tight_layout(pad=0)

    plt.figure(figsize=(12,4))
    plt.rcParams['axes.prop_cycle'] = plt.cycler('color', plt.cm.inferno(np.linspace(0,0.85,3)))
    plt.subplot(1,3,1)
    plt.plot(t[:500]*1000, ac[:500,i_plot]/ac[0,i_plot])
    plt.xlim([-5,150]);
    plt.yticks([0,np.exp(-1), 1], ['0','e','1'])
    plt.xlabel('lag time (ms)'); plt.ylabel('autocorrelation')
    plt.legend(['%i ms'%np.round(tt) for tt in t_ds[i_plot]*1000], frameon=False, loc='upper right', title='decay time constant')

    plt.rcParams['axes.prop_cycle'] = plt.cycler('color',color)
    plt.subplot(1,3,2)
    plt.loglog(f_axis[:f_to_plot], PSD[:,:f_to_plot].T);
    plt.scatter(knee_freq,P_knee, c=color, s=50, edgecolor='k', zorder=100)
    plt.yticks([]); plt.tick_params('y', which='minor', left=False, labelleft=False)
    plt.xlabel('frequency (Hz)'); plt.ylabel(r'power ($V^2/Hz$)')
    plt.xticks([1, 10, 100], ['1','10','100']);
    plt.xlim([1,200]);

    plt.subplot(1,3,3)
    plt.scatter(t_ds*1000,taus*1000, c=color, s=50, edgecolor='k', zorder=100)
    plt.xlim([0,t_ds.max()*1200]);plt.ylim([0,t_ds.max()*1200])
    plt.plot(plt.xlim(), plt.xlim(), 'k--', alpha=0.8);
    plt.xlabel('simulated time constant (ms)'); plt.ylabel('fit time constant (ms)')
    plt.annotate('r = %.2f'%stats.spearmanr(t_ds,taus)[0]**2, xy=(0.75, 0.25), xycoords='axes fraction')
    plt.tight_layout()

    findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.
    findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.

    # defining the simulation, analysis, and plotting function
    def sim_timescale_schematic(tau_sim=0.025, osc_freq=10.5, rel_osc_amp=0.2):
        """
        simulates the neural signal with a given decay timescale and oscillation

        """
        # set fit frequency range
        fit_range=[1,100]
        plt_inds = np.arange(fit_range[0],fit_range[1]+1)
        # simulation time and sampling frequency
        T, fs = 120, 2000
        # PSD parameters
        nperseg = int(fs)
        noverlap = int(0.5*nperseg)

        # simulate signal with given timescale and oscillation
        # Define the components of the combined signal to simulate
        components = {'sim_synaptic_current' : {'tau_d' : tau_sim},
                  'sim_bursty_oscillation' : {'freq' : osc_freq}}
        component_variances = [1, rel_osc_amp]
        x = sim.sim_combined(T, fs, components, component_variances)
        t = np.arange(0,T, 1/fs)

        # compute autocorrelation function and PSD
        autocor = acf(x, nlags=int(fs), fft=True)
        f_axis, psd = spectral.compute_spectrum(x, fs, nperseg=nperseg, noverlap=noverlap)

        # fit with spectral parametrization and get fit values
        # and compute timescale
        ff = FOOOF(max_n_peaks=2, aperiodic_mode='knee', verbose=False)
        ff.fit(f_axis, psd, fit_range)
        offset, knee, exp = ff.get_params('aperiodic_params')
        knee_freq = knee**(1./exp)
        tau_fit = 1./(2*np.pi*knee_freq)
        knee_power = 10**offset/(knee+f_axis**exp)[np.where(f_axis==np.round(knee_freq))[0]]
        fit_spectrum = 10**offset/(knee+f_axis**exp)[plt_inds]

        ### plotting ###
        fig = plt.figure(figsize=(12,8))
        gs = GridSpec(3,2, figure=fig)

        # plot time series
        ax1 = fig.add_subplot(gs[0,:])
        ax1.plot(t[:1000],x[:1000])
        ax1.set_xlabel('time (s)'); ax1.set_ylabel('voltage (au)');
        ax1.set_xlim([0,t[1000]])
        if tau_sim==0.042: print('you found easter egg! nice.')

        # plot autocorrelation
        ax2 = fig.add_subplot(gs[1:,0])
        ax2.plot(t[:1001], autocor[:1001], label='data autocorrelation', lw=2, alpha=0.8)
        ax2.axvline(tau_sim, ls='--', lw=2, label='true tau: %.2f ms'%(tau_sim*1000))
        ax2.axvline(tau_fit, color='r', lw=4, alpha=0.5, label='fit tau: %.2f ms'%(tau_fit*1000))
        ax2.set_xticks([0,tau_sim, 0.1,0.2]); ax2.set_xlim([0,0.2])
        ax2.set_ylim([autocor.min(),1])
        ax2.set_xlabel('lag time (s)'); ax2.set_ylabel('acf');
        ax2.legend()

        # plot spectrum
        ax3 = fig.add_subplot(gs[1:,1])
        ax3.loglog(f_axis[plt_inds], psd[plt_inds], lw=2, label='data PSD')
        ax3.loglog(f_axis[plt_inds], 10**ff.fooofed_spectrum_, 'r-', alpha=0.4, lw=5, label='fit')
        ax3.plot(knee_freq, knee_power, 'ro', ms=20, mec='k', alpha=0.8, label='knee frequency')
        ax3.set_xlim([1,None])
        ax3.set_xlabel('frequency (Hz)'); ax3.set_ylabel('power (V^2/Hz)');
        ax3.legend()

        plt.tight_layout()

    # compute and plotting function is defined in the cell above
    # here, try different timescale values for the simulation by changing `tau_sim`
    # there is stochasticity in the simulation so run it multiple times
    # try to stay within 0.005 to 0.15 seconds

    # you can also change the frequency and relative power of the bursty oscillation
    # to see how the autocorrelation is corrupted by the oscillatory component
    sim_timescale_schematic(tau_sim = 0.02, osc_freq=22.5, rel_osc_amp=0.1)

    # also see here for great autocorrelation gif: https://twitter.com/saydnay/status/1355228493089361921

    # load example macaque data
    data_load = np.load('./data/fig1D_data.npz')
    data, fs = data_load['data'], data_load['fs']
    data_load.close()

    # compute autocorrelation function
    max_lag=250
    t_ac = np.arange(0,max_lag+1)/fs
    ac = acf(data, nlags=max_lag, fft=True)
    # compute psd
    fit_range=[1,80]
    plt_inds = np.arange(fit_range[0],fit_range[1]+1)
    faxis, psd = spectral.compute_spectrum(data,fs,avg_type='median')
    # fit fooof with knee
    fok = FOOOF(max_n_peaks=2, aperiodic_mode='knee', verbose=False)
    fok.fit(faxis, psd, fit_range)
    offset, knee, exp = fok.get_params('aperiodic_params')
    kfreq, tau = convert_knee_val(knee,exp)

    # plot time series
    plt.figure(figsize=(8,2))
    plt.plot(np.arange(0,fs)/fs, data[int(fs*1):int(fs*2)])
    plt.xlim([0,1]);
    plt.xlabel('time (s)');plt.ylabel('voltage (au)');
    plt.title('example macaque ECoG')
    plt.tight_layout()

    plt.figure(figsize=(8,4))
    # plot acf
    plt.subplot(1,2,1)
    plt.plot(t_ac*1000, ac, 'k', lw=2)
    plt.axhline(0, lw=1, ls='--', color='k')
    plt.xlim([0,250])
    plt.xlabel('lag time (ms)'); plt.ylabel('autocorrelation')
    # plot psd
    plt.subplot(1,2,2)
    plt.loglog(faxis[plt_inds], psd[plt_inds], 'k', lw=5, alpha=0.3)
    plt.loglog(faxis[plt_inds], 10**fok.fooofed_spectrum_, 'k-')
    plt.loglog(faxis[plt_inds], 10**offset/(knee+faxis**exp)[plt_inds], '--r', lw=2)
    plt.xlim([1,100])
    plt.xlabel('frequency (Hz)');plt.ylabel('power');
    plt.xticks([1, 10, 100], ['1','10', '100']);plt.yticks([]);
    plt.legend(['data', 'full fit', 'aperiodic fit'], frameon=False)
    plt.yticks([]); plt.tick_params('y', which='minor', left=False, labelleft=False)
    plt.xlabel('frequency (Hz)'); plt.ylabel(r'power ($V^2/Hz$)')
    # plot dot at knee frequency
    plt.plot(kfreq, 10**offset/(knee+faxis**exp)[np.where(faxis==np.round(kfreq))[0]], 'ro', ms=10, alpha=0.8)
    plt.tight_layout()
    # print('knee frequency: %.3fHz, time constant: %.3fms, exponent: %.3f'%(kfreq, 1000*tau, exp))

## Results

### Neuronal timescale can be inferred from the frequency domain

Neural time series often exhibit time-lagged correlation (i.e., autocorrelation), where future values are partially predictable from past values, and predictability decreases with increasing time lags. For demonstration, we simulate the aperiodic (non-rhythmic) component of ECoG recordings by convolving Poisson population spikes with exponentially decaying synaptic kernels with varying decay constant (Figure 1B). Empirically, the degree of self-similarity is characterized by the autocorrelation function (ACF), and ‘timescale’ is defined as the time constant (τ) of an exponential decay function (e−t τ) fit to the ACF, i.e., the time it takes for the autocorrelation to decrease by a factor of e (Figure 1C).

Equivalently, we can estimate timescale in the frequency domain from the power spectral density (PSD). PSDs of neural time series often follow a Lorentzian function of the form 1 f k 2+f 2, where power is approximately constant until the ‘knee frequency’ (fk, Figure 1D), then decays following a power law. This approach is similar to the one presented in 14Chaudhuri et al.2017, but here we further allow the power law exponent (fixed at two in the equation above) to be a free parameter representing variable scale-free activity 45He et al.2010 65Miller et al.2009 75Podvalny et al.2015 89Voytek et al.2015. We also simultaneously parameterize oscillatory components as Gaussians peaks, allowing us to remove their effect on the power spectrum, providing more accurate estimates of the knee frequency. From the knee frequency of the aperiodic component, neural timescale (decay constant) can then be computed exactly as τ=1 2 π f k.

Compared to fitting exponential decay functions in the time domain (e.g., 68Murray et al.2014)—which can be biased even without the presence of additional components 104Zeraati et al.2020—the frequency domain approach is advantageous when a variable power law exponent and strong oscillatory components are present, as is often the case for neural signals (example of real data in Figure 1D). While the oscillatory component can corrupt naive measurement of τ as time for the ACF to reach 1/e (Figure 1D, inset), it can be more easily accounted for and removed in the frequency domain as Gaussian-like peaks. This is especially important considering neural oscillations with non-stationary frequencies. For example, a broad peak in the power spectrum (e.g., ~10 Hz in bandwidth in Figure 1D) represents drifts in the oscillation frequency over time, which is easily accounted for with a single Gaussian, but requires multiple cosine terms to capture well in the autocorrelation. Therefore, in this study, we apply spectral parameterization to extract timescales from intracranial recordings 20Donoghue et al.2020. We validate this approach on PSDs computed from simulated neural time series and show that the extracted timescales closely match their ground-truth values (Figure 1E).

### Timescales follow anatomical hierarchy and are ~10 times faster than spiking timescales

Applying this technique, we infer a continuous gradient of neuronal timescales across the human cortex by analyzing a large dataset of human intracranial (ECoG) recordings of task-free brain activity 29Frauscher et al.2018. The MNI-iEEG dataset contains 1 min of resting state data across 1772 channels from 106 patients (13–62 years old, 48 females) with variable coverages, recorded using either surface strip/grid or stereoEEG electrodes, and cleaned of visible artifacts. Figure 2A shows example data traces along the cortical hierarchy with increasing timescales estimated from their PSDs (Figure 2B; circles denote fitted knee frequency). Timescales from individual channels were extracted and projected from MNI coordinates onto the left hemisphere of HCP-MMP1.0 surface parcellation 38Glasser et al.2016 for each patient using a Gaussian-weighted mask centered on each electrode. While coverage is sparse and idiosyncratic in individual patients, it does not vary as a function of age, and when pooling across the entire population, 178 of 180 parcels have at least one patient with an electrode within 4 mm (Figure 2—figure supplement 1A–F).

Figure 2.

#### Timescale increases along the anatomical hierarchy in humans and macaques.

(A) Example time series from five electrodes along the human cortical hierarchy (M1: primary motor cortex; SMC: supplementary motor cortex; OFC: orbitofrontal cortex; ACC: anterior cingulate cortex; MTL: medial temporal lobe), and (B) their corresponding power spectral densities (PSDs) computed over 1 min. Circle and dashed line indicate the knee frequency for each PSD, derived from the aperiodic component fits (inset). Data: MNI-iEEG database, N = 106 participants. (C) Human cortical timescale gradient (left) falls predominantly along the rostrocaudal axis, similar to T1w/T2w ratio (right; z-scored, in units of standard deviation). Colored dots show electrode locations of example data. (D) Neuronal timescales are negatively correlated with cortical T1w/T2w, thus increasing along the anatomical hierarchy from sensory to association regions (Spearman correlation; p-value corrected for spatial autocorrelation, Figure 2—figure supplement 2A–C). (E) Example PSDs from macaque ECoG recordings, similar to (B) (LIP: lateral intraparietal cortex; LPFC: lateral prefrontal cortex; S1 and S2: primary and secondary somatosensory cortex). PSDs are averaged over electrodes within each region (inset of [F]). Data: Neurotycho, N = 8 sessions from two animals. (F) Macaque ECoG timescales track published single-unit spiking timescales 68Murray et al.2014 in corresponding regions (error bars represent mean ± s.e.m). Inset: ECoG electrode map of one animal and selected electrodes for comparison. (G) ECoG-derived timescales are consistently correlated with (left), and ~10 times faster than (right), single-unit timescales across individual sessions. Hollow markers: individual sessions; shapes: animals; solid circles: grand average from (F).

    # load data
    # timescale data
    df_tau = pd.read_csv('./data/df_tau.csv', index_col=0)
    df_tau.columns=['timescale (ms)', 'log10 timescale (ms)']
    # timescale SAP-surrogates
    msr_nulls = pd.read_csv('./data/df_tau_shuffles.csv', index_col=0).values.T

    # structural data
    df_struct = pd.read_csv('./data/df_structural_avg.csv', index_col=0)
    df_struct.columns = df_struct.columns.str.upper()

    # load macroparcel data
    df_macro = pd.read_csv('./data/df_human_features_macro.csv')
    df_macro.columns = df_macro.columns.str.upper()
    df_macro.columns = ['index', 'timescale (ms)'] + list(df_macro.columns[2:])
    df_macro['timescale (ms)'] = df_macro['timescale (ms)']*1000

    ## load example ECoG data
    data_load = np.load('./data/fig2AB_data.npz')
    fs, data, psds, f_axis = data_load['fs'], data_load['data'], data_load['psds'], data_load['f_axis']
    data_load.close()

    # plot example time series
    labels = ['M1', 'SMC', 'OFC', 'ACC', 'MTL']
    plt.figure(figsize=(8,4))

    # these are the indices of the channels in the MNI dataset
    # but only these 5 are included here, you can find the whole .mat file online
    plt_inds = [125, 220, 1123, 1666, 573]
    c_ord = [3, 2, 0, 4, 5]

    plt.figure(figsize=(13.5,3))
    for i, i_p in enumerate(plt_inds):
        plt.plot(stats.zscore(data[:, i])-4*i, color=C_ORD[c_ord[i]], label=labels[i])

    plt.xticks([0, fs*5], ['0', '5']); plt.yticks([]);
    plt.xlabel('time (s)');plt.ylabel('voltage (au)');
    plt.xlim([0, fs*5])
    plt.legend(loc='lower left', bbox_to_anchor= (1, 0), ncol=1, frameon=False, handletextpad=0.5)
    plt.tight_layout()

    # plot example PSDs and fit
    plt.figure(figsize=(12,4))
    for i, i_p in enumerate(plt_inds):
        fit_range=[1,70]
        fok = FOOOF(max_n_peaks=3, aperiodic_mode='knee', verbose=False)
        fok.fit(f_axis, psds[:,i], fit_range)
        offset, knee, exp = fok.get_params('aperiodic_params')
        kfreq, tau = convert_knee_val(knee,exp)
        ap_spectrum = (10**offset/(knee+f_axis**exp))

        color = C_ORD[c_ord[i]]

        plt.subplot(1,3,1)
        plt.loglog(f_axis[1:100], psds[1:100,i]/psds[2,i], lw=2, color=color)
        plt.axvline(kfreq, ls='--', color=color, lw=2, alpha=0.3)
        plt.plot(kfreq, 23, 'o', color=color, ms=10, label=labels[i])
        plt.xticks([1, 10, 100], ['1', '10','100']); plt.yticks([]); plt.xlabel('frequency (Hz)'); plt.ylabel(r'power ($V^2/Hz$)')
        plt.xlim([1,100]); plt.ylim([None,30])

        plt.subplot(1,3,2)
        plt.loglog(f_axis[2:100], ap_spectrum[2:100]/ap_spectrum[1], '-', color=color, lw=4, alpha=0.8, label=labels[i])
        plt.xticks([1, 10, 100], ['1', '10','100']); plt.yticks([]); plt.xlabel('frequency (Hz)'); plt.ylabel(r'power ($V^2/Hz$)')
        plt.xlim([2,100]);
        plt.legend(loc='lower left', bbox_to_anchor= (0, 0.01), ncol=1, frameon=False, handletextpad=1)

    x=stats.zscore(df_struct['T1T2'])
    y=df_tau['log10 timescale (ms)']
    rho, pv, pv_perm, rho_null = compute_perm_corr(x,y.values,msr_nulls)
    m,b,_,_,_ = stats.linregress(x,y)
    plt.subplot(1,3,3)
    plt.plot(x, y, 'o', color=C_ORD[3], alpha=0.5, ms=5)
    XL= np.array(plt.xlim())
    plt.plot(XL,XL*m+b, '--', lw=2, color=C_ORD[3], alpha=0.8)
    plt.xlabel(r'z-scored T1w/T2w ($\sigma$)');  plt.ylabel('timescale (ms)');
    plt.tick_params('y', which='minor', left=False, labelleft=False)
    plt.yticks(np.log10(np.arange(10,60,10)), (np.arange(10, 60, 10)).astype(int))
    s = sig_str(rho, pv_perm, form='text')
    plt.annotate(s, xy=(0.55, 0.75), xycoords='axes fraction')
    plt.tight_layout()

    <Figure size 576x288 with 0 Axes>

    ## Not currently working, the plot_MMP function is throwing a KeyError

    # brain_fig_size = (6,6)
    # struct_query = 'T1T2' # try 'T1T2' or 'THICKNESS' or any of the genes in Fig 3. supplement 2
    # # plot on brains
    # plt.figure(figsize=brain_fig_size)
    # plot_MMP(y, minmax=[1,np.log10(50)], title='log10 timescale (ms)', cmap='inferno');
    # plt.show();

    # plt.figure(figsize=brain_fig_size)
    # plot_MMP(x, title=struct_query, cmap='pink', minmax=[-2, 4]);
    # plt.show();

    # get mean population time constants, values grabbed from Murray et al., 2014
    cell_ts = {'MT':[77.,64.], 'LIP':[138., 91.], 'LPFC':[184.,180.,195.,162.], 'OFC':[176.,188.], 'ACC':[313.,340.,257.], 'S1':[65.], 'S2':[149.]}
    cell_ts_avg = {k: np.array([np.mean(np.array(v)), stats.sem(np.array(v))]) for k,v in cell_ts.items()}

    # electrode indices for each of the corresponding areas in each monkey
    loc_inds_chibi = {'MT':[3,4,109], 'LIP':[10,11], 'LPFC':[14,15,25,26], 'OFC':[23,34,45], 'ACC':[52,57,58,59], 'S1':[9,19], 'S2':[95,108]}
    loc_inds_george = {'MT':[4,13,22], 'LIP':[10,11,20,21], 'LPFC':[15,24,25,26], 'OFC':[45,66], 'ACC':[52,57,58,59], 'S1':[18,19,30], 'S2':[1,2,9,108]}
    loc_inds = {'Chibi': loc_inds_chibi, 'George': loc_inds_george}

    area_ord = [3,1,2,0,4,6,5] # color order to match Murray figure

    ##### plot example PSDs ######
    data_load = np.load('./data/fig2E_data.npz')
    psds, f_axis = data_load['psds'], data_load['f_axis']
    data_load.close()
    plt.figure(figsize=(8,4))
    for i_r, (reg, inds) in enumerate(loc_inds_chibi.items()):

        psds_reg = np.mean(np.log10(psds[np.array(inds)-1]),0)
        psds_reg = psds_reg-psds_reg[40]
        fit_range=[1,70]
        plt_inds = np.arange(fit_range[0],fit_range[1]+1)
        fok = FOOOF(max_n_peaks=3, aperiodic_mode='knee', verbose=False)
        fok.fit(f_axis, 10**psds_reg, fit_range)
        offset, knee, exp = fok.get_params('aperiodic_params')
        kfreq, tau = convert_knee_val(knee,exp)
        ap_spectrum = np.log10((10**offset/(knee+f_axis**exp)))

        plt.subplot(1,2,1)
        plt.semilogx(f_axis[3:100], psds_reg[3:100] , color=C_ORD[area_ord[i_r]], label=reg, lw=2)
        plt.axvline(kfreq, ls='--', color=C_ORD[area_ord[i_r]], lw=2, alpha=0.3)
        plt.plot(kfreq, 1.9, 'o', color=C_ORD[area_ord[i_r]], ms=10)
        plt.xticks([10, 100], ['10','100']); plt.yticks([]); plt.xlabel('frequency (Hz)'); plt.ylabel(r'power ($V^2/Hz$)')
        plt.xlim([3,100]); plt.ylim([None,2])

        plt.subplot(1,2,2)
        plt.semilogx(f_axis[3:100], ap_spectrum[3:100], '-', color=C_ORD[area_ord[i_r]], lw=4, alpha=0.8)
        plt.xticks([10, 100], ['10','100']); plt.yticks([]); plt.xlabel('frequency (Hz)'); plt.ylabel(r'power ($V^2/Hz$)')
        plt.xlim([3,100]); plt.ylim([None,2])

    ##### plot stats and results ######
    # load ecog results dataframe
    df_combined = pd.read_csv('./data/df_macaque.csv', index_col=0)

    # define querying condition
    feature = 'tau'
    cond_query = 'EyesOpen'
    df_cond = df_combined[df_combined['cond']==cond_query]

    # plot
    plt.figure(figsize=(5,4))
    ecog_ts_avg = {}

    # plot per session average across electrodes
    for i, k in enumerate(cell_ts.keys()):
        sesh_mkrs = [] # hack to save the session marker for next plot
        ecog_ts_avg[k] = []
        for s in df_cond['session_id'].unique():
            df_sesh = df_cond[df_cond['session_id']==s]
            patient = df_sesh['patient'].iloc[0]

            # loc_inds has the ecog electrode indices that fall into each area
            region_inds = loc_inds[patient][k]
            marker = 's' if patient == 'George' else '^'
            sesh_mkrs.append(marker)
            ecog_ts_sess_avg = df_sesh.loc[df_sesh['chan'].isin(region_inds)].mean()[feature]*1e3 # use ms
            ecog_ts_avg[k].append(ecog_ts_sess_avg)

    # plot grand average
    for i,k in enumerate(cell_ts.keys()):
        plt.errorbar(cell_ts_avg[k][0], np.mean(ecog_ts_avg[k]), xerr=cell_ts_avg[k][1], yerr=stats.sem(ecog_ts_avg[k]), fmt='o', color=C_ORD[area_ord[i]], ms=10, label=k)

    # fit & plot line
    ts_mat = np.array([(cell_ts_avg[k][0], np.mean(ecog_ts_avg[k])) for k in cell_ts.keys()])
    m,b,r,pv,stderr = stats.linregress(ts_mat)
    XL = np.array(plt.xlim())
    plt.plot(XL, m*XL+b, 'k--', lw=1)

    rho, pv = stats.spearmanr(ts_mat[:,0], ts_mat[:,1])
    s = sig_str(rho, pv, form='text')
    plt.annotate(s, xy=(0.05, 0.8), xycoords='axes fraction')

    plt.tick_params('x', which='minor', bottom=False, labelbottom=False)
    plt.tick_params('y', which='minor', left=False, labelleft=False)
    plt.xlim([-10,400]);plt.ylim([8,50]);

    plt.legend(loc='lower left', bbox_to_anchor= (0.9, 0), ncol=1, frameon=False, handletextpad=0.05)
    plt.xlabel(r'$\tau_{spiking}$ (ms)', fontsize=16);
    plt.ylabel(r'$\tau_{ECoG}$ (ms)', fontsize=16); #plt.title(cond_query)
    plt.tight_layout();

    ### compute and plot across-trial stats ###
    ecog_ts_mat = np.array([ecog_ts_avg[k] for k in cell_ts.keys()])
    cell_ts_mat = np.array([cell_ts_avg[k][0] for k in cell_ts.keys()])
    session_stats = np.array([stats.linregress(cell_ts_mat, ecog_ts_mat[:,i]) for i in range(ecog_ts_mat.shape[1])])
    session_rhos = np.array([stats.spearmanr(cell_ts_mat, ecog_ts_mat[:,i]) for i in range(ecog_ts_mat.shape[1])])

    plt.figure(figsize=(4,4))
    plt.subplot(1,2,1)
    for s_i, s in enumerate(sesh_mkrs): plt.plot(np.random.randn(1)/10, session_rhos[s_i,0],mec='k', mfc='w', ms=6, marker=s, alpha=0.9)
    plt.plot(0, rho,'ko', ms=10, alpha=0.6)
    plt.xlim([-0.5,1]); plt.ylim([0,1]);
    plt.xticks([]); plt.ylabel(r'Spearman $\rho$')

    plt.subplot(1,2,2)
    for s_i, s in enumerate(sesh_mkrs): plt.plot(np.random.randn(1)/10, session_stats[s_i,0]*100, mec='k', mfc='w', ms=6, marker=s)
    plt.plot(0, m*100,'ko', ms=12, alpha=0.6)
    plt.xlim([-0.5,1]);
    plt.yticks([0, 10, 20], ['0%', '10%', '20%'])
    plt.xticks([]); plt.ylabel(r'$\tau_{ECoG} : \tau_{spiking}$ (%)', fontsize=16)
    plt.tight_layout();

    # print('open to explore data!')

    struct_query = 'PVALB'

    # make temp dataframe because loading the whole thing in altair is slow af
    df_query = pd.concat([df_tau['timescale (ms)'], df_struct[struct_query]], axis=1)
    rho, pv = stats.spearmanr(df_tau['timescale (ms)'], df_struct[struct_query])

    alt.Chart(df_query.reset_index()).mark_circle(size=200).encode(
        alt.X(struct_query, scale=alt.Scale(zero=False)),
        alt.Y('timescale (ms)', scale=alt.Scale(zero=False, type='log', base=2, domain=(10, 64))),
        color=alt.value('black'),
        tooltip=['index']
    ).properties(
        width=600,
        height=600
    ).configure_axis(
        labelFontSize=20,
        titleFontSize=20
    ).configure_title(fontSize=24).interactive()

    # print('open to explore data!')
    struct_query = 'PVALB'

    # make temp dataframe because loading the whole thing in altair is slow af
    df_query = pd.concat([df_macro['index'], df_macro['timescale (ms)'], df_macro[struct_query]], axis=1)
    rho, pv = stats.spearmanr(df_macro['timescale (ms)'], df_macro[struct_query])

    alt.Chart(df_query).mark_circle(size=200).encode(
        alt.X(struct_query, scale=alt.Scale(zero=False)),
        alt.Y('timescale (ms)', scale=alt.Scale(zero=False, domain=(10, 40))),
        color=alt.value('black'),
        tooltip=['index']
    ).properties(
        width=600,
        height=600
    ).configure_axis(
        labelFontSize=20,
        titleFontSize=20
    ).configure_title(fontSize=24).interactive()

Figure 2—figure supplement 1.

#### MNI-iEEG dataset electrode coverage.

(A) Per-parcel Gaussian-weighted mask values showing how close the nearest electrode was to a given HCP-MMP1.0 parcel for each participant. Brighter means closer, 0.5 corresponds to the nearest electrode being 4 mm away. (B) Maximum mask weight for each parcel across all participants. Most parcels have electrodes very close by in at least one participant across the entire participant pool. (C) The number of valid HCP-MMP parcels each participant has above the confidence threshold of 0.5 is uncorrelated with age. (D) Cortical map of the number of participants with confidence above threshold at each parcel. Sensorimotor, frontal, and lateral temporal regions have the highest coverage. (E) Cortical map of the average age of participants with confidence above threshold at each parcel. (F) Age distribution of participants with confidence above threshold at each parcel. Average age per parcel (red line) is relatively stable while age distribution varies from parcel to parcel (each subject is a black dot). (G) Average neuronal timescale when further aggregating the 180 Glasser parcels into 21 macro-regions (mean ± s.e.m across parcels within the macro-region).

Figure 2—figure supplement 2.

#### Comparison of spatial autocorrelation-preserving null map generation methods.

(A) Distributions of Spearman correlation values between empirical T1w/T2w map and 1000 spatial-autocorrelation preserving null timescale maps generated using Moran Spectral Randomization (MSR), spatial variogram fitting (VF), and spin permutation. Red dashed line denotes correlation between empirical timescale and T1w/T2w maps, p-values indicate two-tailed significance, i.e., proportion of distribution with values more extreme than empirical correlation. (B) Spatial variogram for empirical timescale map (black) and 10 null maps (blue) generated using MSR (left) and VF (right). Inset shows distribution of distances between pairs of HCP-MMP parcels. (C) Distribution of Spearman correlations between empirical and 1000 null timescale maps generated using MSR (green) and VF (red), showing similar levels of correlation between empirical and null maps for both methods.

Figure 2—figure supplement 3.

#### Cortical thickness.

Cortical thickness from the HCP dataset is positively correlated with neuronal timescale (left) and negatively correlated with T1w/T2w, i.e., thicker brain regions have longer (slower) timescales and less gray matter myelination, corresponding to higher-order association areas.

    # comparing cortical thickness to timescale
    x = df_struct['THICKNESS'].values
    y = df_tau['log10 timescale (ms)']

    # --- group average correlation ----
    rho, pv, pv_perm, rho_null = compute_perm_corr(x,y.values,msr_nulls)
    m,b,_,_,_ = stats.linregress(x,y.values)

    plt.figure(figsize=(8,4))
    plt.subplot(1,2,1)
    plt.plot(x, y, 'o', color=C_ORD[4], alpha=0.5, ms=5);
    XL= np.array(plt.xlim())
    plt.plot(XL,XL*m+b, '--', lw=2, color=C_ORD[4], alpha=0.8)
    plt.xlabel('cortical thickness (mm)');  plt.ylabel('timescale (ms)');
    plt.tick_params('y', which='minor', left=False, labelleft=False)
    plt.yticks(np.log10(np.arange(10,60,10)), (np.arange(10, 60, 10)).astype(int))
    s = sig_str(rho, pv_perm, form='text')
    plt.annotate(s, xy=(0.65, 0.85), xycoords='axes fraction')
    plt.tight_layout()

Figure 2—figure supplement 4.

#### Macaque ECoG and single-unit coverage.

(A) Locations of 180-electrode ECoG grid from two animals in the Neurotycho dataset; colors correspond to locations used for comparison with single-unit timescales. (B) Electrode indices of the sampled areas from the two animals, corresponding to those colored in (A).

Across the human cortex, timescales of fast electrophysiological dynamics (~10–50 ms) predominantly follow a rostrocaudal gradient (Figure 2C, circles denote location of example data from 2A). Consistent with numerous accounts of a principal cortical axis spanning from primary sensory to association regions 46Hilgetag and Goulas2020 64Margulies et al.2016 98Wang2020, timescales are shorter in sensorimotor and early visual areas, and longer in association regions, especially cingulate, ventral/medial frontal, and medial temporal lobe (MTL) regions (Figure 2—figure supplement 1G shows further pooling into 21 labeled macro-regions). We then compare the timescale gradient to the average T1w/T2w map from the Human Connectome Project, which captures gray matter myelination and indexes the proportion of feedforward vs. feedback connections between cortical regions, thus acting as a noninvasive proxy of connectivity-based anatomical hierarchy 10Burt et al.2018 39Glasser and Van Essen2011. We find that neuronal timescales are negatively correlated with T1w/T2w across the entire cortex (Figure 2D, ρ = −0.47, p<0.001; corrected for spatial autocorrelation [SA], see Materials and methods and Figure 2—figure supplement 2A–C for a comparison of correction methods), such that timescales are shorter in more heavily myelinated (i.e., lower-level, sensory) regions. Timescales are also positively correlated with cortical thickness (Figure 2—figure supplement 3, ρ = 0.37, p=0.035)—another index of cortical hierarchy that is itself anti-correlated with T1w/T2w. Thus, we observe that neuronal timescales lengthen along the human cortical hierarchy, from sensorimotor to association regions.

While surface ECoG recordings offer much broader spatial coverage than extracellular single-unit recordings, they are fundamentally different signals: ECoG and field potentials largely reflect integrated synaptic and other transmembrane currents across many neuronal and glial cells, rather than putative action potentials from single neurons (12Buzsáki et al.2012; Figure 1A, yellow box). Considering this, we ask whether timescales measured from ECoG in this study (τ ECoG) are related to single-unit spiking timescales along the cortical hierarchy (τ spiking). To test this, we extract neuronal timescales from task-free ECoG recordings in macaques 69Nagasaka et al.2011 and compare them to a separate dataset of single-unit spiking timescales from a different group of macaques 68Murray et al.2014 (see Figure 2—figure supplement 4 for electrode locations). Consistent with τ spiking estimates 68Murray et al.2014 99Wasmuht et al.2018, τ ECoG also increase along the macaque cortical hierarchy. While there is a strong correspondence between spiking and ECoG timescales (Figure 2F; ρ = 0.96, p<0.001)—measured from independent datasets—across the macaque cortex, τ ECoG are ~10 times faster than τ spiking and are conserved across individual sessions (Figure 2G). This suggests that neuronal spiking and transmembrane currents have distinct but related timescales of fluctuations, and that both are hierarchically organized along the primate cortex.

### Synaptic and ion channel genes shape timescales of neuronal dynamics

Next, we identify potential cellular and synaptic mechanisms underlying timescale variations across the human cortex. Theoretical accounts posit that NMDA-mediated recurrent excitation coupled with fast inhibition 13Chaudhuri et al.2015 95Wang2008 93Wang1999, as well as cell-intrinsic properties 23Duarte and Morrison2019 37Gjorgjieva et al.2016 59Koch et al.1996, are crucial for shaping neuronal timescales. While in vitro and in vivo studies in model organisms 83van Vugt et al.2020 97Wang et al.2013 can test these hypotheses at the single-neuron level, causal manipulations and large-scale recordings of neuronal networks embedded in the human brain are severely limited. Here, we apply an approach analogous to multimodal single-cell profiling 5Bomkamp et al.2019 and examine the transcriptomic basis of neuronal dynamics at the macroscale.

Leveraging whole-cortex interpolation of the Allen Human Brain Atlas bulk mRNA expression 42Gryglewski et al.2018 43Hawrylycz et al.2012, we project voxel-wise expression maps onto the HCP-MMP1.0 surface parcellation, and find that the neuronal timescale gradient overlaps with the dominant axis of gene expression (i.e., first principal component of 2429 brain-related genes) across the human cortex (Figure 3A, ρ = −0.60, p<0.001; see Figure 3—figure supplement 1 for similar results with all 18,114 genes). Consistent with theoretical predictions (Figure 3B), timescales significantly correlate with the expression of genes encoding for NMDA (GRIN2B) and GABA-A (GABRA3) receptor subunits, voltage-gated sodium (SCN1A) and potassium (KCNA3) ion channel subunits, as well as inhibitory cell-type markers (parvalbumin, PVALB), and genes previously identified to be associated with single-neuron membrane time constants (PRR5) 5Bomkamp et al.2019 (all Spearman correlations corrected for SA in gradients).

Figure 3.

#### Timescale gradient is linked to expression of genes related to synaptic receptors and transmembrane ion channels across the human cortex.

(A) Timescale gradient follows the dominant axis of gene expression variation across the cortex (z-scored PC1 of 2429 brain-specific genes, arbitrary direction). (B) Timescale gradient is significantly correlated with expression of genes known to alter synaptic and neuronal membrane time constants, as well as inhibitory cell-type markers, but (C) members within a gene family (e.g., NMDA receptor subunits) can be both positively and negatively associated with timescales, consistent with predictions from in vitro studies. (D) Macroscale timescale-transcriptomic correlation captures association between RNA-sequenced expression of the same genes and single-cell timescale properties fit to patch clamp data from two studies, and the correspondence improves for genes (separated by quintiles) that are more strongly correlated with timescale (solid: N = 170 81Tripathy et al.2017, dashed: N = 4168 genes 5Bomkamp et al.2019; horizontal lines: correlation across all genes from the two studies, ρ = 0.36 and 0.25, p<_0.001 for both). (E) T1w/T2w gradient is regressed out from timescale and gene expression gradients, and a partial least squares (PLS) model is fit to the residual maps. Genes with significant PLS weights (filled blue boxes) compared to spatial autocorrelation (SA)-preserved null distributions are submitted for gene ontology enrichment analysis (GOEA), returning a set of significant GO terms that represent functional gene clusters (filled green boxes). (F) Enriched genes are primarily linked to potassium and chloride transmembrane transporters, and GABA-ergic synapses; genes specifically with strong negative associations further over-represent transmembrane ion exchange mechanisms, especially voltage-gated potassium and cation transporters. Branches indicate GO items that share higher-level (parent) items, e.g., _voltage-gated cation channel activity is a child of cation channel activity in the molecular functions (MF) ontology, and both are significantly associated with timescale. Color of lines indicate curated ontology (BP—biological process, CC—cellular components, or MF). Dotted, dashed, and solid lines correspond to analysis performed using all genes or only those with positive or negative PLS weights. Spatial correlation p-values in (A–C) are corrected for SA (see Materials and methods; asterisks in (B,D) indicate p<0.05, 0.01, 0.005, and 0.001 respectively; filled markers in (C,D) indicate p<0.05).

    # perform PCA on brain genes
    from sklearn.decomposition import PCA
    brain_gene_symbols = pd.read_csv('./data/df_brain_genes_symbols.csv', index_col=0).values[:,1].tolist()
    df_genes = df_struct[brain_gene_symbols]

    n_pcs = 10
    gene_pca = PCA(n_pcs)
    gene_grad = gene_pca.fit_transform(df_genes.values)
    df_pc_weights = pd.DataFrame(gene_pca.components_.T, index=df_genes.columns, columns=['pc%i'%i for i in range(1, n_pcs+1)])

    # compute correlation
    x = stats.zscore(gene_grad[:,0])
    y = df_tau['log10 timescale (ms)']
    rho, pv, pv_perm, rho_null = compute_perm_corr(x,y.values,msr_nulls)
    m,b,_,_,_ = stats.linregress(x,y.values)

    plt.figure(figsize=(12,4))

    plt.subplot(1,3,1)
    plt.plot(x, y, 'o', color=C_ORD[0], alpha=0.5, ms=5)
    plt.xlim([-2.75, 2.75]);
    XL= np.array([-2.5, 2.5])
    plt.plot(XL,XL*m+b, '--', lw=2, color=C_ORD[0], alpha=0.8)
    plt.xlabel(r'gene expression PC1 ($\sigma$)'); plt.ylabel('timescale (ms)');
    plt.yticks(np.log10(np.arange(10,60,10)), (np.arange(10, 60, 10)).astype(int))
    plt.tick_params('y', which='minor', left=False, labelleft=False)
    s = sig_str(rho, pv_perm, form='text')
    plt.annotate(s, xy=(0.05, 0.05), xycoords='axes fraction')

    plt.subplot(1,3,2)
    plt.bar(range(1,n_pcs+1), gene_pca.explained_variance_ratio_, fc='k', alpha=0.5)
    plt.xticks([1, 10], ['1', '10']);
    plt.xlabel('principal component'); plt.ylabel('proportion of var. explained');

    plt.subplot(1,3,3)
    all_pc_rhos = np.array([compute_perm_corr(x, y, msr_nulls)[:3] for x in gene_grad.T])
    plt.bar(range(1,n_pcs+1), np.abs(all_pc_rhos[:,0]), fc='k', alpha=0.5)
    for i in range(1,11):
        plt.annotate('*'*sum(all_pc_rhos[i-1,2]<[0.05, 0.01, 0.005, 0.001]), (i, abs(all_pc_rhos[i-1,0])+0.005), horizontalalignment='center')
    plt.xticks([1, 10], ['1', '10']);
    plt.xlabel('principal component'); plt.ylabel(r'absolute $\rho$');
    plt.tight_layout()

    ## Not currently working, the plot_MMP function is throwing a ValueError:
    # all the input arrays must have same number of dimensions, but the array at index 0 has 1 dimension(s) and the array at index 1 has 2 dimension(s)

    # plt.figure(figsize=(6,6))
    # plot_MMP(x, bp=0, minmax=[-2.5,2.5], title='PC1', cmap='viridis')

    xlb = ['GRIN2B', 'GABRA3','SCN1A', 'KCNA3','PVALB', 'PRR5']

    plt.figure(figsize=(4,6))
    y = df_tau['log10 timescale (ms)']
    xs = [df_struct[g] for g in xlb]

    for i_x, x in enumerate(xs):
        rho, pv, pv_perm, rho_null = compute_perm_corr(x,y.values,msr_nulls)
        m,b,_,_,_ = stats.linregress(x,y.values)

        plt.subplot(3,2,i_x+1)
        plt.plot(x, y, '.', color=C_ORD[0], alpha=0.5, ms=5);
        XL= np.array([-2.5, 2.5])
        plt.xlim([-2.75, 2.75]);
        plt.plot(XL,XL*m+b, '--', lw=2, color=C_ORD[0], alpha=0.8)
        plt.gca().set_xticklabels([])
        plt.gca().set_yticklabels([])
        plt.annotate(sig_str(rho, pv_perm, form='*'), xy=(0.02,0.05), xycoords='axes fraction', fontsize=14);
        plt.xlabel(xlb[i_x], fontsize=14, labelpad=0)

    plt.tight_layout()

    findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.

    geneset_nmda = ['GRIN1', 'GRINA', 'GRIN2A','GRIN2B','GRIN2C','GRIN2D','GRIN3A', 'GRIN3B'] # nmda receptor
    geneset_gabra = ['GABRA1','GABRA2','GABRA3','GABRA4','GABRA5','GABRA6'] # GABA-A alpha subchannels
    geneset_sodium = ['SCN1A',  'SCN2A',  'SCN3A', 'SCN4A',  'SCN5A', 'SCN7A', 'SCN8A', 'SCN9A', 'SCN10A']#,'SCN1B','SCN2B','SCN3B',SCN4B'] # sodium ion channels
    geneset_potassium = ['KCNA1','KCNA2','KCNA3','KCNA4','KCNA5','KCNA6'] # GABA-A alpha subchannels
    geneset_inh = ['CORT', 'CALB1', 'CALB2', 'SST', 'PVALB', 'CCK', 'NPY', 'PNOC', 'VIP'] # inhibitory markers
    geneset_sctau = ['CELF6', 'PRR5', 'FAM81A', 'LRRC4C','OXTR', 'CTXN1', 'ENC1', 'AKAIN1'] # single-cell membrane time constant
    gene_families = [geneset_nmda, geneset_gabra, geneset_sodium, geneset_potassium, geneset_inh, geneset_sctau]

    df_tau_gene_corrfam = pd.DataFrame([], index=sum(gene_families, []), columns=['rho', 'pv', 'pv_adj'])
    for i_g, g in df_tau_gene_corrfam.iterrows():
        df_tau_gene_corrfam.loc[i_g] = compute_perm_corr(df_struct[i_g].values, y.values, msr_nulls)[0:3]

    # plot all correlations
    plt.figure(figsize=(3,4))
    plt.axvline(0, color='k', lw=1)
    for i_g, gs in enumerate(gene_families):
        set_x = df_tau_gene_corrfam.loc[gs]['rho']
        set_y = np.random.randn(len(gs))/15+(5-i_g)

        # color based on permutation significance
        set_color = [C_ORD[0] if p<0.05 else 'w' for p in df_tau_gene_corrfam.loc[gs]['pv_adj']]
        plt.axhline(5-i_g, lw=0.2, color='k')
        plt.scatter(set_x, set_y, alpha=0.5, s=50, ec = C_ORD[0], c = set_color)

    plt.xlabel(r'$\rho$', labelpad=-15); plt.xticks([-0.5, 0.5])
    plt.yticks(np.arange(len(gene_families)), ['sc-tau','inhibitory','K+ chan.','Na+ chan.','GABAA', 'NMDA']);
    plt.tight_layout()

    # another view, plot families of correlations as hbar
    color = plt.cm.RdBu_r(np.linspace(0,1,200))
    plt.figure(figsize=(18,4))
    for i_g, g_subset in enumerate(gene_families):
        plt.subplot(1,len(gene_families),i_g+1)
        geneset_plot = df_tau_gene_corrfam.loc[g_subset]
        for i_p in range(len(geneset_plot)):
            rho, pv_perm = df_tau_gene_corrfam.loc[geneset_plot.iloc[i_p].name][['rho', 'pv_adj']]
            plt.barh(i_p, rho, ec='k', fc=color[int((1+rho*1.5)*100)])
            s = np.sum(pv_perm<=np.array([0.05, 0.01, 0.005, 0.001]))*'*'
            plt.text(rho, i_p-0.4, s, fontsize=18, horizontalalignment='left' if rho>0 else 'right')

        plt.plot([0,0], plt.ylim(), 'k')
        plt.ylim([-0.5,len(geneset_plot)-0.5])
        plt.yticks(range(len(geneset_plot)), geneset_plot.index.values, rotation=0, ha='right', va='center', rotation_mode='anchor', fontsize=14)
        plt.tick_params(axis='y', which=u'both',length=0)
        plt.xlim([-1,1]); plt.xlabel(r'$\rho$', labelpad=-15); plt.xticks([-1,1])

    plt.tight_layout()

    findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.

    def collect_micro_macro_corr(df_micro_corr, df_macro_corr, col_names):
        sc_feats, corr_metric,  prop_col, gene_col = col_names
        micro_macro_corr = []
        for i_f, feat in enumerate(sc_feats):
            match_genes = []
            for i_g, g in df_micro_corr[df_micro_corr[prop_col]==feat].iterrows():
                if g[gene_col].upper() in df_macro_corr.index:
                    match_genes.append([i_g, g[gene_col].upper()])

            match_genes = np.array(match_genes, dtype='object')
            micro_macro_corr.append([match_genes[:,1], df_micro_corr.loc[match_genes[:,0]][corr_metric].values, df_macro_corr.loc[match_genes[:,1]]['rho'].values])

        return np.hstack(micro_macro_corr), micro_macro_corr

    # collect micro/macro correlations of relevant genes
    df_bomkamp = pd.read_csv('./data/bomkamp_online_table1.txt', index_col=0)
    df_tripathy = pd.read_csv('./data/tripathy_tableS3.csv', index_col=0)

    df_genes = df_struct[df_struct.columns[2:]]
    df_tau_gene_corrall = pd.DataFrame([stats.spearmanr(g.values, y.values) for g_i, g in df_genes.iteritems()], columns=['rho','pv'], index=df_genes.columns)
    df_tau_gene_macro = df_tau_gene_corrall[df_tau_gene_corrall['pv']<0.05]
    #df_tau_gene_macro = df_tau_gene_corrall_rmvt1t2[df_tau_gene_corrall_rmvt1t2['pv']<0.05]

    col_names = [['tau', 'ri', 'cap'], 'beta_gene' , 'property', 'gene_symbol']
    micro_macro_bomkamp, mmc_bomkamp_split = collect_micro_macro_corr(df_bomkamp, df_tau_gene_macro, col_names)

    col_names = [['Tau', 'Rin', 'Cm'], 'DiscCorr', 'EphysProp', 'GeneSymbol']
    micro_macro_tripathy, mmc_tripathy_split = collect_micro_macro_corr(df_tripathy, df_tau_gene_macro, col_names)

    plt.figure(figsize=(4,4))
    pct_bins = np.arange(0,100,20)
    lss = ['o-', '^--']
    labels=['Tripathy17', 'Bomkamp18']

    for i_m, micro_macro in enumerate([micro_macro_tripathy, micro_macro_bomkamp]):
        meta_corr = []
        micro_macro_abscorr = abs(micro_macro[2].astype(float))
        bins = np.percentile(micro_macro_abscorr, pct_bins)
        corr_quant_inds = np.digitize(micro_macro_abscorr, bins)
        for i in np.unique(corr_quant_inds):
            meta_corr.append(stats.pearsonr(micro_macro[1,corr_quant_inds==i],micro_macro[2,corr_quant_inds==i]))

        meta_corr=np.array(meta_corr)
        r, pv = stats.pearsonr(micro_macro[1],micro_macro[2])
        plt.plot(pct_bins+pct_bins[1], meta_corr[:,0], lss[i_m], color=C_ORD[0], ms=8, mfc='w', alpha=0.8, label=labels[i_m])
        plt.plot(pct_bins[meta_corr[:,1]<0.05]+pct_bins[1], meta_corr[meta_corr[:,1]<0.05,0], lss[i_m][0], color=C_ORD[0], ms=8, alpha=0.8)
        plt.plot(pct_bins[[0,-1]]+pct_bins[1], [r, r], ls=lss[i_m][1:], color='k', alpha=0.7, zorder=-20)

    plt.xticks([20,100],['1/5', '5/5']);
    plt.yticks(np.arange(-0.25, 1, 0.25))
    plt.legend(frameon=False, fontsize=14, loc=[0,0.8])
    plt.xlabel(r'macroscale absolute $\rho$ quintile'); plt.ylabel('single-cell vs. macroscale\nassociation correlation (r)')
    plt.tight_layout()

Figure 3—figure supplement 1.

#### Transcriptomic principal component analysis results.

(A) Proportion of variance explained by the top 10 principal components (PCs) of brain-specific genes (top) and all AHBA genes (bottom). (B) Absolute Spearman correlation between timescale map and top 10 PCs from brain-specific or full gene dataset. Asterisks indicate resampled significance while accounting for spatial autocorrelation; **** indicate p_<_0.001. Top PCs explain similar amounts of variance, while only PC1 in both cases is significantly correlated with timescale.

Figure 3—figure supplement 2.

#### Individual timescale-gene correlations magnitudes.

Correlation between timescale and expression of genes from Figure 3C, with gene symbols labeled and grouped into functional families for ease of interpretation.

More specifically, in vitro electrophysiological studies have shown that, for example, increased expression of receptor subunit 2B extends the NMDA current time course 27Flint et al.1997, while 2A expression shortens it 66Monyer et al.1994. Similarly, the GABA-A receptor time constant lengthens with increasing a3:a1 subunit ratio 24Eyre et al.2012. We show that these relationships are recapitulated at the macroscale, where neuronal timescales positively correlate with GRIN2B and GABRA3 expression and negatively correlate with GRIN2A and GABRA1 (Figure 3C). These results demonstrate that timescales of neural dynamics depend on specific receptor subunit combinations with different (de)activation timescales 22Duarte et al.2017 37Gjorgjieva et al.2016, in addition to broad excitation–inhibition interactions 34Gao et al.2017 98Wang2020 94Wang2002. Notably, almost all genes related to voltage-gated sodium and potassium ion channel alpha-subunits—the main functional subunits—are correlated with timescale, while all inhibitory cell-type markers except parvalbumin have strong positive associations with timescale (Figure 3C and Figure 3—figure supplement 2).

We further test whether single-cell timescale-transcriptomic associations are captured at the macroscale as follows: for a given gene, we can measure how strongly its expression correlates with membrane time constant parameters at the single-cell level using patch-clamp and RNA sequencing (scRNASeq) data 5Bomkamp et al.2019 81Tripathy et al.2017. Analogously, we can measure its macroscopic transcriptomic-timescale correlation using the cortical gradients above. If the association between the expression of this gene and neuronal timescale is preserved at both levels, then the correlation across cells at the microscale should be similar to the correlation across cortical regions at the macroscale. Comparing across these two scales for all previously identified timescale-related genes in two studies (N = 170 81Tripathy et al.2017 and 4168 5Bomkamp et al.2019 genes), we find a significant correlation between the strength of association at the single-cell and macroscale levels (Figure 3D, horizontal black lines; ρ = 0.36 and 0.25 for the two datasets, p<0.001 for both). Furthermore, genes with stronger associations to timescale tend to conserve this relationship across single-cell and macroscale levels (Figure 3D, separated by macroscale correlation magnitude). Thus, the association between cellular variations in gene expression and cell-intrinsic temporal dynamics is captured at the macroscale, even though scRNAseq and microarray data represent entirely different measurements of gene expression.

While we have shown associations between cortical timescales and genes suspected to influence neuronal dynamics, these data present an opportunity to discover additional novel genes that are functionally related to timescales through a data-driven approach. However, since transcriptomic variation and anatomical hierarchy overlap along a shared macroscopic gradient 10Burt et al.2018 49Huntenburg et al.2018 64Margulies et al.2016, we cannot specify the role certain genes play based on their level of association with timescale alone: gene expression differences across the cortex first result in cell-type and connectivity differences, sculpting the hierarchical organization of cortical anatomy. Consequently, anatomy and cell-intrinsic properties jointly shape neuronal dynamics through connectivity differences 13Chaudhuri et al.2015 19Demirtaş et al.2019 and expression of ion transport and receptor proteins with variable activation timescales, respectively. Therefore, we ask whether variation in gene expression still accounts for variation in timescale beyond the principal structural gradient, and if associated genes have known functional roles in biological processes (BP) (schematic in Figure 3E). To do this, we first remove the contribution of anatomical hierarchy by linearly regressing out the T1w/T2w gradient from both timescale and individual gene expression gradients. We then fit partial least squares (PLS) models to simultaneously estimate regression weights for all genes 101Whitaker et al.2016, submitting those with significant associations for gene ontology enrichment analysis (GOEA) 58Klopfenstein et al.2018.

We find that genes highly associated with neuronal timescales are preferentially related to transmembrane ion transporter complexes, as well as GABAergic synapses and chloride channels (see Figure 3F and Supplementary file 1 for GOEA results with brain genes only, and Supplementary file 2 for all genes). When restricted to positively associated genes only (expression increases with timescales), one functional group related to phosphatidate phosphatase activity is uncovered, including the gene PLPPR1, which has been linked to neuronal plasticity 78Savaskan et al.2004—a much slower timescale physiological process. Conversely, genes that are negatively associated with timescale are related to numerous groups involved in the construction and functioning of transmembrane transporters and voltage-gated ion channels, especially potassium and other inorganic cation transporters. To further ensure that these genes specifically relate to neuronal timescale, we perform the same enrichment analysis with T1w/T2w vs. gene maps as a control. The control analysis yields no significant GO terms when restricted to brain-specific genes (in contrast to Figure 3F), while repeating the analysis with all genes does yield significant GO terms related to ion channels and synapses, but are much less specific to those (see Supplementary file 3), including a variety of other gene clusters associated with general metabolic processes, signaling pathways, and cellular components (CC). This further strengthens the point that removing the contribution of T1w/T2w aids in identifying genes that are more specifically associated with neurodynamics, suggesting that inhibition 80Teleńczuk et al.2017—mediated by GABA and chloride channels—and voltage-gated potassium channels have prominent roles in shaping neuronal timescale dynamics at the macroscale level, beyond what is expected based on the anatomical hierarchy alone.

### Timescales lengthen in working memory and shorten in aging

Finally, having shown that neuronal timescales are associated with stable anatomical and gene expression gradients across the human cortex, we turn to the final question of the study: are cortical timescales relatively static, or are they functionally dynamic and relevant for human cognition? While previous studies have shown hierarchical segregation of task-relevant information corresponding to intrinsic timescales of different cortical regions 3Baldassano et al.2017 15Chien and Honey2020 47Honey et al.2012 76Runyan et al.2017 77Sarafyazd and Jazayeri2019 99Wasmuht et al.2018, as well as optimal adaptation of behavioral timescales to match the environment 33Ganupuru et al.2019 71Ossmy et al.2013, evidence for functionally relevant changes in regional neuronal timescales is lacking. Here, we examine whether timescales undergo short- and long-term shifts during working memory maintenance and aging, respectively.

We first analyze human ECoG recordings from parietal, frontal (PFC and orbitofrontal cortex [OFC]), and medial temporal lobe (MTL) regions of patients (N = 14) performing a visuospatial working memory task that requires a delayed cued response (Figure 4A; 51Johnson et al.2018). Neuronal timescales were extracted for pre-stimulus baseline and memory maintenance delay periods (900 ms, both stimulus-free). Replicating our previous result in Figure 2—figure supplement 1G, we observe that baseline neuronal timescales follow a hierarchical progression across association regions, where parietal cortex (PC), PFC, OFC, and MTL have gradually longer timescales (pairwise Mann–Whitney U-test, Figure 4B). If neuronal timescales track the temporal persistence of information in a functional manner, then they should expand during delay periods. Consistent with our prediction, timescales in all regions are ~20% longer during delay periods (Figure 4C; Wilcoxon rank-sum test). Moreover, only timescale changes in the PFC are significantly correlated with behavior across participants, where longer delay-period timescales relative to baseline are associated with better working memory performance (Figure 4D, ρ = 0.75, p=0.003). No other spectral features in the recorded brain regions show consistent changes from baseline to delay periods while also significantly correlating with individual performance, including the 1/f-like spectral exponent, narrowband theta (3–8 Hz), and high-frequency (high gamma; 70–100 Hz) activity power (Figure 4—figure supplement 1).

Figure 4.

#### Timescales expand during working memory maintenance while tracking performance, and task-free average timescales compress in older adults.

(A) Fourteen participants with overlapping intracranial coverage performed a visuospatial working memory task, with 900 ms of baseline (pre-stimulus) and delay period data analyzed (PC: parietal, PFC: prefrontal, OFC: orbitofrontal, MTL: medial temporal lobe; n denotes number of subjects with electrodes in that region). (B) Baseline timescales follow hierarchical organization within association regions (*: p<0.05, Mann–Whitney U-test; small dots represent individual participants, large dots and error bar for mean ± s.e.m. across participants). (C) All regions show significant timescale increase during delay period compared to baseline (asterisks represent p<0.05, 0.01, 0.005, 0.001, Wilcoxon signed-rank test). (D) PFC timescale expansion during delay periods predicts average working memory accuracy across participants (dot represents individual participants, mean ± s.e.m. across PFC electrodes within participant); inset: correlation between working memory accuracy and timescale change for all regions. (E) In the MNI-iEEG dataset, participant-average cortical timescales decrease (become faster) with age (n = 71 participants with at least 10 valid parcels, see Figure 4—figure supplement 2B). (F) Most cortical parcels show a negative relationship between timescales and age, with the exception being parts of the visual cortex and the temporal poles (one-sample t-test, t = −7.04, p<0.001; n = 114 parcels where at least six participants have data, see Figure 4—figure supplement 2C).

    df_patient_info = pd.read_csv('./data/fig4D_patientinfo.csv', index_col=0)
    plt.figure(figsize=(16,4))
    #### baseline timescale
    df_mean = pd.read_csv('./data/fig4B.csv', index_col=0)
    region_labels = ['Parietal', 'PFC', 'OFC', 'MTL']
    rl_short = ['PC', 'PFC', 'OFC', 'MTL']
    plt.subplot(1,4,1)
    for i_r, reg in enumerate(region_labels):
        plt.plot([i_r]*len(df_mean)+np.random.randn(len(df_mean))/10, df_mean[reg].values, '.', ms=5, alpha=0.7, color=C_ORD[i_r])
        plt.errorbar(i_r, df_mean[reg].mean(), df_mean[reg].sem(), color=C_ORD[i_r], fmt='o', ms=10, alpha=1)

    plt.xticks(range(len(rl_short)), rl_short)
    plt.yticks(np.arange(0.02,0.07, 0.01), (np.round(np.arange(0.02,0.07, 0.01)*1000)).astype(int))
    plt.ylabel('baseline timescale (ms)')
    plt.xlim([-0.5, 3.5])
    plt.tight_layout()

    # print('-----prestim-----')
    # print('PC-OFC: ', stats.mannwhitneyu(df_mean['Parietal'], df_mean['OFC'], alternative='two-sided'))
    # print('PC-PFC: ', stats.mannwhitneyu(df_mean['Parietal'], df_mean['PFC'], alternative='two-sided'))
    # print('PC-MTL: ', stats.mannwhitneyu(df_mean['Parietal'], df_mean['MTL'], alternative='two-sided'))
    # print('PFC-OFC: ', stats.mannwhitneyu(df_mean['PFC'], df_mean['OFC'], alternative='two-sided'))
    # print('PFC-MTL: ', stats.mannwhitneyu(df_mean['PFC'], df_mean['MTL'], alternative='two-sided'))
    # print('OFC-MTL: ', stats.mannwhitneyu(df_mean['OFC'], df_mean['MTL'], alternative='two-sided'))

    #### pre-post change in timescale
    # print('\n-----pre:post change-----')
    df_mean = pd.read_csv('./data/fig4C_mean.csv', index_col=0)
    df_sem = pd.read_csv('./data/fig4C_sem.csv', index_col=0)
    plt.subplot(1,4,2)
    for i_r, reg in enumerate(region_labels):
        plt.plot([i_r]*len(df_mean)+np.random.randn(len(df_mean))/10, df_mean[reg].values, '.', ms=5, alpha=0.7, color=C_ORD[i_r])
        plt.errorbar(i_r, df_mean[reg].mean(), df_mean[reg].sem(), color=C_ORD[i_r], fmt='o', ms=10, alpha=1)
        pv = stats.wilcoxon(df_mean[reg][~np.isnan(df_mean[reg])])[1]
        # print(reg, pv)
        s = sig_str(0, pv)
        plt.annotate(s.split(' ')[-1], xy=(i_r, 0.145), horizontalalignment='center', fontsize=20)
        for j_r, reg2 in enumerate(region_labels):
            pass # print(reg, '-', reg2, stats.mannwhitneyu(df_mean[reg][~np.isnan(df_mean[reg])], df_mean[reg2][~np.isnan(df_mean[reg2])], alternative='two-sided'))

    plt.xticks(range(len(rl_short)), rl_short)
    # change axis to show percent change
    yt = np.arange(1,1.41,0.2)
    plt.yticks(np.log10(yt), ['%i%%'%int(yt_*100) for yt_ in yt]); plt.ylim(np.log10([yt[0], yt[-1]]))
    plt.ylabel('delay period timescale\n(as % of baseline)')
    plt.xlim([-0.5, 3.5])
    plt.tight_layout()

    #### pre-post change per region
    reg='PFC'
    df_mean = pd.read_csv('./data/fig4C_mean.csv', index_col=0)
    df_sem = pd.read_csv('./data/fig4C_sem.csv', index_col=0)
    y = df_mean[reg][~np.isnan(df_mean[reg])]
    y_err = df_sem[reg][~np.isnan(df_sem[reg])]
    x = df_patient_info[['acc_identity','acc_spatial','acc_temporal']].mean(1)
    rho, pv = stats.spearmanr(x[y.index].values,y.values)

    plt.subplot(1,4,3)
    plt.errorbar(x[y.index].values, y.values, y_err.values, fmt='o', ms=8, alpha=0.8, color=C_ORD[1])
    s = sig_str(rho, pv, form='text')
    plt.annotate(s, xy=(0.7, 0.05), xycoords='axes fraction');
    plt.yticks(np.log10(yt), ['%i%%'%int(yt_*100) for yt_ in yt]); plt.ylim(np.log10([yt[0], yt[-1]]))
    plt.xticks()
    plt.xlim(left=0.68);plt.ylim(bottom=-0.01);
    plt.xlabel('working memory accuracy'); plt.ylabel('PFC delay period timescale (%)')
    plt.tight_layout()

    # show all regions
    plt.subplot(1,4,4)
    for i_r, reg in enumerate(region_labels):
        y = df_mean[reg][~np.isnan(df_mean[reg])]
        rho, pv = stats.spearmanr(x[y.index].values,y.values)
        # print(reg, rho, pv)
        plt.barh(len(region_labels)-i_r, rho, fc=C_ORD[i_r], ec='k', lw=1, alpha=0.8)
        s = sig_str(rho, pv)

    plt.xticks([-1,1]); plt.xlabel(r'$\rho$', labelpad=-15)
    plt.yticks(range(1,5), region_labels[::-1]);

    findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.

    ### normalized averages timescales over parcels per patient
    df_patient = pd.read_csv('./data/fig4E_patient.csv', index_col=0)
    tau_subj_normed = pd.read_csv('./data/fig4E.csv', index_col=0)
    n_parcel_thresh = 10
    n_main = n_parcel_thresh
    x = df_patient['age'][df_patient['coverage']>=n_parcel_thresh]
    y = tau_subj_normed[df_patient['coverage']>=n_parcel_thresh]
    rho, pv = stats.spearmanr(x,y)
    s = sig_str(rho, pv, form='text') +'\nn = %i'%len(x)
    m,b,_,_,_ = stats.linregress(x.values, y['0'].values)

    plt.figure(figsize=(8,4))
    plt.subplot(1,2,1)
    plt.plot(x, y, 'ok', alpha=0.5, ms=5);
    XL= np.array(plt.xlim())
    plt.plot(XL,XL*m+b, 'k--', lw=2, alpha=0.8)
    plt.xlabel('age (years)'); plt.ylabel('parcel-normalized timescale');
    plt.xticks([10, 35, 60],None); plt.yticks([0,1])
    plt.annotate(s, xy=(0.05, 0.025), xycoords='axes fraction');

    ### distribution of correlations
    parcel_corr = pd.read_csv('./data/fig4F.csv', index_col=0)
    x = parcel_corr.values
    # t-test on whether parcels have positive or negative correlation in aggregate
    tval, tt_pv = stats.ttest_1samp(x, 0, nan_policy='omit')
    plt.subplot(1,2,2)
    plt.hist(x[~np.isnan(x)], bins=np.arange(-1,1.1,0.1), color='k', alpha=0.5)
    plt.axvline(np.nanmean(x), color=C_ORD[3], ls='--')
    plt.xlabel(r'$\rho$');plt.ylabel('parcels (%)', labelpad=-15);
    s = sig_str(np.nanmean(x), tt_pv, form='text', corr_letter=r'$\bar{\rho}$')+'\nn = %i'%len(x[~np.isnan(x)])
    plt.annotate(s, xy=(0.45, 0.8), xycoords='axes fraction')
    plt.xticks([-1,0,1],None);plt.yticks(np.array([0,15])/100*sum(~np.isnan(x)), ['0%', '15%']);
    plt.tight_layout()

Figure 4—figure supplement 1.

#### Spectral correlates of working memory performance.

(A) Difference between delay and baseline periods for 1/f-exponent, timescale (same as main Figure 4C but absolute units on y-axis, instead of percentage), theta power, and high-frequency power. (B) Spearman correlation between spectral feature difference and working memory accuracy across participants, same features as in (A). *p<0.05, **p<0.01, ***p<0.005 in (A, B). (C) Scatter plot of other significantly correlated spectral features from (B).

Figure 4—figure supplement 2.

#### Parameter sensitivity for timescale-aging analysis.

(A) Cortex-averaged timescale is independent of parcel coverage across participants. (B) Sensitivity analysis for the number of valid parcels a participant must have in order to be included in analysis for main Figure 4E (red). As threshold increases (more stringent), fewer participants satisfy the criteria (right) but correlation between participant age and timescale remains robust (left). (C) Sensitivity analysis for the number of valid participants a parcel must have in order to be included in analysis for main Figure 4F. As threshold increases (more stringent), fewer parcels satisfy the criteria (right) but average correlation across all parcels remains robust (left, error bars denote s.e.m of distribution as in Figure 4F).

While timescales are consistent with the anatomical and gene expression hierarchy at a snapshot, brain structure itself is not static over time, undergoing many slower, neuroplastic changes during early development and throughout aging in older populations. In particular, aging is associated with a broad range of functional and structural changes, such as working memory impairments 89Voytek et al.2015 96Wang et al.2011, as well as changes in neuronal dynamics 89Voytek et al.2015 90Voytek and Knight2015 96Wang et al.2011 and cortical structure 18de Villers-Sidani et al.2010, such as the loss of slow-deactivating NMDA receptor subunits 73Pegasiou et al.2020. Since neuronal timescales may support working memory maintenance, we predict that timescales would shorten across the lifespan, in agreement with the observed cognitive and structural deteriorations. To this end, we leverage the wide age range in the MNI-iEEG dataset (13–62 years old) and probe average cortical timescales for each participant as a function of age. Since ECoG coverage is sparse and nonuniform across participants, simply averaging across parcels within individual participants confounds the effect of aging with the spatial effect of cortical hierarchy. Instead, we first normalize each parcel by its max value across all participants before averaging within participants, excluding those with fewer than 10 valid parcels (N = 71 of 106 subjects remaining; results hold for a large range of threshold values, Figure 4—figure supplement 2B). We observe that older adults have faster neuronal timescales (ρ = −0.31, p=0.010; Figure 4E), and that timescales shorten with age in most areas across the cortex (Figure 4F, t = −7.04, p<0.001; 114 out of 189 parcels where at least six participants have data, see Figure 4—figure supplement 2C). This timescale compression is especially prominent in sensorimotor, temporal, and medial frontal regions. These results support our hypothesis that neuronal timescales, estimated from transmembrane current fluctuations, can rapidly shift in a functionally relevant manner, as well as slowly—over decades—in healthy aging.

## Discussion

Theoretical accounts and converging empirical evidence predict a graded variation of neuronal timescales across the human cortex 13Chaudhuri et al.2015 49Huntenburg et al.2018 98Wang2020, which reflects functional specialization and implements hierarchical temporal processing crucial for complex cognition 57Kiebel et al.2008. This timescale gradient is thought to emerge as a consequence of cortical variations in cytoarchitecture, as well as both macroscale and microcircuit connectivity, thus serving as a bridge from brain structure to cognitive function 55Kanai and Rees2011. In this work, we infer the timescale of non-rhythmic transmembrane current fluctuations from invasive human intracranial recordings and test those predictions explicitly. We discuss the implications and limitations of our findings below.

### Multiple quantities for neuronal timescale and anatomical hierarchy

We first find that neuronal timescales vary continuously across the human cortex and coincide with the anatomical hierarchy, with timescales increasing from primary sensory and motor to association regions. While we use the continuous T1w/T2w gradient as a surrogate measure for anatomical hierarchy, there are multiple related but distinct perspectives on what ‘cortical hierarchy’ means, including, for example, laminar connectivity patterns from tract tracing data 26Felleman and Van Essen1991 85Vezoli et al.2020, continuous (and latent-space) gradients of gene expression and microarchitectural features 49Huntenburg et al.2018, and network connectivity scales (see review of 46Hilgetag and Goulas2020)—with most of these following a graded sensorimotor-to-association area progression. Similarly, it is important to note that there exist many different quantities that can be considered as characteristic neuronal timescales across several spatial scales, including membrane potential and synaptic current timescales 22Duarte et al.2017, single-unit spike-train timescales 68Murray et al.2014, population code timescales 76Runyan et al.2017, and even large-scale circuit timescales measured from the fMRI BOLD signal 100Watanabe et al.2019. We show here that timescales inferred from ECoG are consistently approximately 10 times faster than single-unit spiking timescales in macaques, corroborating the fact that field potential signals mainly reflect fast transmembrane and synaptic currents 12Buzsáki et al.2012, whose timescales are related to, but distinct from, single-unit timescales measured in previous studies 21Dotson et al.2018 70Ogawa and Komatsu2010 99Wasmuht et al.2018.

Because field potential fluctuations are driven by currents from both locally generated and distal inputs, our results raise questions on how and when these timescales interact to shape downstream spiking dynamics. Furthermore, while we specifically investigate here the aperiodic timescale, which corresponds to the exponential decay timescale measured in previous studies, recent work has shown a similar gradient of oscillatory timescales (i.e., frequency) along the anterior–posterior axis of the human cortex 63Mahjoory et al.2020. Based on the similarity of these gradients and known mechanisms of asynchronous and oscillatory population dynamics (e.g., balance of excitation and inhibition in generating gamma oscillations and the asynchronous irregular state in cortical circuits 8Brunel2000 9Brunel and Wang2003), we speculate that timescales of oscillatory and aperiodic neural dynamics may share (at least partially) circuit mechanisms at different spatial scales, analogous to the relationship between characteristic frequency and decay constant in a damped harmonic oscillator model.

### Collinearity and surrogate nature of postmortem gene expression gradients

Using postmortem gene expression data as a surrogate for protein density, transcriptomic analysis uncovers the potential roles that transmembrane ion transporters and synaptic receptors play in establishing the cortical gradient of neuronal timescales. The expression of voltage-gated potassium channel, chloride channel, and GABAergic receptor genes, in particular, are strongly associated with the spatial variation of neuronal timescale. Remarkably, we find that electrophysiology-transcriptomic relationships discovered at the single-cell level, through patch-clamp recordings and single-cell RNA sequencing, are recapitulated at the macroscale between bulk gene expression and timescales inferred from ECoG. That being said, it is impossible to make definitive causal claims with the data presented in this study, especially considering the fact that several microanatomical features, such as gray matter myelination and cortical thickness, follow similar gradients across the cortex 10Burt et al.2018. To discover genes specifically associated with timescale while accounting for the contribution of the overlapping anatomical hierarchy, we linearly regress out the T1w/T2w gradient from both timescale and gene expression gradients. Although this procedure does not account for any nonlinear contributions from anatomy, gene enrichment control analysis using T1w/T2w instead of timescales further demonstrates that the discovered genes—transmembrane ion transporters and inhibitory synaptic receptors—are more specifically associated with the timescale gradient, over and above the level predicted by anatomical hierarchy alone. From these results, we infer that potassium and chloride ion channels, as well as GABAergic receptors, may play a mechanistic role in altering the timescale of transmembrane currents at the macroscopic level.

However, this interpretation rests on the key assumption that mRNA expression level is a faithful representation of the amount of functional proteins in a given brain region. In general, gene expression levels are highly correlated with the percentage of cells expressing that gene within brain regions 60Lein et al.2007. Therefore, on a population level, the regional density of a particular ion channel or receptor protein is high if bulk mRNA expression is high. Furthermore, recent works have shown that neurotransmitter receptor density measured via autoradiography in postmortem brains follows similar cortical gradients 41Goulas et al.2020, and that gene expression levels of neurotransmitter receptors (e.g., 5HT) are strongly correlated with ligand binding potential measured via PET 42Gryglewski et al.2018. Thus, as a first order approximation, receptor gene expression is an adequate surrogate for receptor protein density in the brain at the macroscale, though the relationship between mRNA expression and their transport and translation into channel proteins, the process of incorporating those proteins into membranes and synapses, and how these gene expression maps can be related to other overlapping macroscopic gradients are complex issues (see e.g., 28Fornito et al.2019; 62Liu et al.2016). Thus, our analyses represent an initial data-mining process at the macroscopic level, which should motivate further studies in investigating the precise roles voltage-gated ion channels and synaptic inhibition play in shaping functional neuronal timescales through causal manipulations, complementary to existing lines of research focusing on NMDA activation and recurrent circuit motifs.

### Structural constraints vs. behaviorally required flexibility in timescale

Finally, we show that neuronal timescales are not static, but can change both in the short and long terms. Transmembrane current timescales across multiple association regions, including parietal, frontal, and medial temporal cortices, increase during the delay period of a working memory task, consistent with the emergence of persistent spiking during working memory delay. Working memory performance across individuals, however, is predicted by the extent of timescale increase in the PFC only. This further suggests that behaviorally relevant neural activity may be localized despite widespread task-related modulation 74Pinto et al.2019, even at the level of neuronal membrane fluctuations. In the long term, we find that neuronal timescale shortens with age in most cortical regions, linking age-related synaptic, cellular, and connectivity changes—particularly those that influence neuronal integration timescale—to the compensatory posterior-to-anterior shift of functional specialization in healthy aging 17Davis et al.2008.

These results raise further questions regarding contrasting, and potentially complementary, aspects of neuronal timescale: on the one hand, task-free timescales across the cortex are shaped by relatively static macro- and microarchitectural properties (Figures 2 and 3); on the other hand, timescales are dynamic and shift with behavioral demand (Figure 4). While long-term structural changes in the brain can explain shifts in neuronal timescales throughout the aging process, properties such as ion channel protein density probably do not change within seconds during a working memory task. We speculate that structural properties may constrain dynamical properties (such as timescale) to a possible range within a particular brain region and at different spatial scales, while task requirements, input statistics, short-term synaptic plasticity, and neuromodulation can then shift timescale within this range. We posit, then, that only shifts in dynamics within the area of relevance (i.e., PFC for working memory) are indicative of task performance, consistent with recent ideas of computation-through-dynamics 91Vyas et al.2020. Nevertheless, which neuromodulatory and circuit mechanisms are involved in shifting local timescales, and how timescales at different spatial scales (e.g., synaptic, neuronal, population) interact to influence each other remain open questions for future investigation 6Breakspear2017 22Duarte et al.2017 31Freeman2000 32Freeman and Erwin2008 37Gjorgjieva et al.2016 79Shine et al.2019.

### Conclusion

In summary, we identify consistent and converging patterns between transcriptomics, anatomy, dynamics, and function across multiple datasets of different modalities from different individuals and multiple species. As a result, evidence for these relationships can be supplemented by more targeted approaches such as imaging of receptor metabolism. Furthermore, the introduction and validation of an open-source toolbox 20Donoghue et al.2020 for inferring timescales from macroscale electrophysiological recordings potentially allows for the noninvasive estimation of neuronal timescales, using widely accessible tools such as EEG and MEG. These results open up many avenues of research for discovering potential relationships between microscale gene expression and anatomy with the dynamics of neuronal population activity at the macroscale in humans.

## Materials and methods

### Inferring timescale from autocorrelation and PSD

Consistent with previous studies, we define ‘neuronal timescale’ as the exponential decay time constant (τ) of the empirical ACF, or lagged correlation 47Honey et al.2012 68Murray et al.2014. τ can be naively estimated to be the time it takes for the ACF to decrease by a factor of e when there are no additional long-term, scale-free, or oscillatory processes, or by fitting a function of the form f(t)=e−t τ and extracting the parameter τ. Equivalently, the PSD is the Fourier Transform of the ACF via Wiener–Khinchin theorem 56Khintchine1934 and follows a Lorentzian function of the form L(f)=A k+f χ for approximately exponential-decay processes, with χ=2 exactly when the ACF is solely composed of an exponential decay term, though it is often variable and in the range between 2 and 6 for neural time series 20Donoghue et al.2020 65Miller et al.2009 75Podvalny et al.2015 89Voytek et al.2015. Timescale can be computed from the parameter k as τ=1 2 π f k, where f k≈k 1/χ is approximated to be the ‘knee frequency’, at which a bend or knee in the power spectrum occurs, and equality holds when χ=2.

### Computing PSD

PSDs are estimated using a modified Welch’s method, where short-time windowed Fourier transforms (STFT) are computed from the time series, but the median is taken across time instead of the mean (in conventional Welch’s method) to minimize the effect of high-amplitude transients and artifacts 50Izhikevich et al.2018. Custom functions for this can be found in NeuroDSP 16Cole et al.2019, a published and open-source digital signal processing toolbox for neural time series (neurodsp.spectral.compute_spectrum). For simulated data, Neurotycho macaque ECoG, and MNI-iEEG datasets, we use 1 s long Hamming windows with 0.5 s overlap. To estimate single-trial PSDs for the working memory ECoG dataset (CRCNS Johnson-ECoG 51Johnson et al.2018; 52Johnson et al.2018), we simply apply Hamming window to 900 ms long epoched time series and compute the squared magnitude of the windowed Fourier transform.

### Spectral parametrization

We apply spectral parameterization 20Donoghue et al.2020 to extract timescales from PSDs. Briefly, we decompose log-power spectra into a summation of narrowband periodic components—modeled as Gaussians—and an aperiodic component—modeled as a generalized Lorentzian function centered at 0 Hz (L(f) above). For inferring decay timescale, this formalism can be practically advantageous when a strong oscillatory or variable power-law (χ) component is present, as is often the case for neural signals. While oscillatory and power-law components can corrupt naive measurements of τ as time for the ACF to reach 1/e, they can be easily accounted for and ignored in the frequency domain as narrowband peaks and 1/f-exponent fit. We discard the periodic components and infer timescale from the aperiodic component of the PSD. For a complete mathematical description of the model, see 20Donoghue et al.2020.

### Simulation and validation

We simulate the aperiodic background component of neural field potential recordings as autocorrelated stochastic processes by convolving Poisson population spikes with exponentially decaying synaptic kernels with predefined decay time constants (neurodsp.sim.sim_synaptic_current). PSDs of the simulated data are computed and parameterized as described above, and we compare the fitted timescales with their ground-truth values.

### Macaque ECoG and single-unit timescales data

Macaque single-unit timescales are taken directly from values reported in Figure 1c of 68Murray et al.2014. Whole-brain surface ECoG data (1000 Hz sampling rate) is taken from the Neurotycho repository 69Nagasaka et al.2011 103Yanagawa et al.2013, with eight sessions of 128-channel recordings from two animals (George and Chibi, four sessions each). Results reported in Figure 2E–G are from ~10 min eyes-open resting periods to match the pre-stimulus baseline condition of single-unit experiments. Timescales for individual ECoG channels are extracted and averaged over regions corresponding to single-unit recording areas from 68Murray et al.2014; Figure 2F inset and Figure 2—figure supplement 3, which are selected visually based on the overlapping cortical map and landmark sulci/gyri. Each region included between two and four electrodes (see Figure 2—figure supplement 3B for selected ECoG channel indices for each region).

### Statistical analysis for macaque ECoG and spiking timescale

For each individual recording session, as well as the grand average, Spearman rank correlation was computed between spiking and ECoG timescales. Linear regression models were fit using the python package scipy 86Virtanen et al.2020 (scipy.stats.linregress) and the linear slope was used to compute the scaling coefficient between spiking and ECoG timescales.

### Variations in neuronal timescale, T1/T2 ratio, and mRNA expression across human cortex

The following sections describe procedures for generating the average cortical gradient maps for neuronal timescale, MR-derived T1w/T2w ratio, and gene expression from the respective raw datasets. All maps were projected onto the 180 left hemisphere parcels of Human Connectome Project’s Multimodal Parcellation 38Glasser et al.2016 (HCP-MMP1.0) for comparison, described in the individual sections. Projection of T1w/T2w and gene expression maps from MNI volumetric coordinates to HCP-MMP1.0 can be found: https://github.com/rudyvdbrink/Surface_projection 82van den Brink2020.

All spatial correlations are computed as Spearman rank correlations between maps. Procedure for computing statistical significance while accounting for SA is described in detail below under the sections 'Spatial statistics' and 'SA modeling'.

### Neuronal timescale map

The MNI Open iEEG dataset consists of 1 min of resting state data across 1772 channels from 106 epilepsy patients (13–62 years old, 58 males, and 48 females), recorded using either surface strip/grid or stereoEEG electrodes, and cleaned of visible artifacts 29Frauscher et al.2018 30Frauscher et al.2018. Neuronal timescales were extracted from PSDs of individual channels, and projected from MNI voxel coordinates onto HCP-MMP1.0 surface parcellation as follows.

For each patient, timescale estimated from each electrode was extrapolated to the rest of the cortex in MNI coordinates using a Gaussian weighting function (confidence mask), w(r)=e−(r 2/α 2), where r is the Euclidean distance between the electrode and a voxel, and α is the distance scaling constant, chosen here such that a voxel 4 mm away has 50% weight (or confidence). Timescale at each voxel is computed as a weighted spatial average of timescales from all electrodes (i) of that patient:

i.e., τ voxel=∑i w(r i)τ i∑i w(r i).

Similarly, each voxel is assigned a confidence rating that is the maximum of weights over all electrodes (_wvoxel(rmin), _of the closest electrode), i.e., a voxel right under an electrode has a confidence of 1, while a voxel 4 mm away from the closest electrode has a confidence of 0.5, etc.

Timescales for each HCP-MMP parcel were then computed as the confidence-weighted arithmetic mean across all voxels that fall within the boundaries of that parcel. HCP-MMP boundary map is loaded and used for projection using NiBabel 7Brett et al.2020. This results in a 180 parcels-by-106 patients timescale matrix. A per-parcel confidence matrix of the same dimensions was computed by taking the maximum confidence over all voxels for each parcel (Figure 2—figure supplement 1A). The average cortical timescale map (gradient) is computed by taking the confidence-weighted average at each parcel across all participants. Note that this procedure for locally thresholded and weighted average is different from projection procedures used for the mRNA and T1w/T2w data due to region-constrained and heterogeneous ECoG electrode sites across participants. While coverage is sparse and idiosyncratic in individual participants, it does not vary as a function of age, and when pooling across the entire population, 178 of 180 parcels have at least one patient with an electrode within 4 mm, with the best coverage in sensorimotor, temporal, and frontal regions (Figure 2—figure supplement 1).

### T1w/T2w ratio and cortical thickness maps

As a measure of structural cortical hierarchy, we used the ratio between T1- and T2-weighted structural MRI, referred to as T1w/T2w map in main text, or the myelin map 10Burt et al.2018 39Glasser and Van Essen2011. Since there is little variation in the myelin map across individuals, we used the group average myelin map of the WU-Minn HCP S1200 release (N = 1096, March 1, 2017 release) provided in HCP-MMP1.0 surface space. For correlation with other variables, we computed the median value per parcel, identical to the procedure for mRNA expression below. Cortical thickness map was similarly generated.

### mRNA expression maps

We used the Allen Human Brain Atlas (AHBA) gene expression dataset 44Hawrylycz et al.2015 43Hawrylycz et al.2012 that comprised postmortem samples of six donors (one female and five males) that underwent microarray transcriptional profiling. Spatial maps of mRNA expression were available in volumetric 2 mm isotropic MNI space, following improved nonlinear registration and whole-brain prediction using variogram modeling as implemented by 42Gryglewski et al.2018. We used whole-brain maps available from 42Gryglewski et al.2018 rather than the native sample-wise values in the AHBA database to prevent bias that could occur due to spatial inhomogeneity of the sampled locations. In total, 18,114 genes were included for analyses that related to the dominant axis of expression across the genome.

We projected the volumetric mRNA expression data onto the HCP-MMP cortical surface using the HCP workbench software (v1.3.1 running on Windows OS 10) with the ‘enclosing’ method and custom MATLAB code (github.com/rudyvdbrink/surface_projection) 82van den Brink2020. The enclosing method extracts for all vertices on the surface the value from enclosing voxels in the volumetric data. Alternative projection methods such as trilinear 3D linear interpolation of surrounding voxels, or ribbon mapping that constructs a polyhedron from each vertex's neighbors on the surface to compute a weighted mean for the respective vertices, yielded comparable values, but less complete cortical coverage. Moreover, the enclosing method ensured that no transformation of the data (nonlinear or otherwise) occurred during the projection process and thus the original values in the volumetric data were preserved.

Next, for each parcel of the left hemisphere in HCP-MMP, we extracted the median vertex-wise value. We used the median rather than the mean because it reduced the contribution of outliers in expression values within parcels. Vertices that were not enclosed by voxels that contained data in volumetric space were not included in the parcel-wise median. This was the case for 539 vertices (1.81% of total vertices). Linear interpolation across empty vertices prior to computing median parcel-wise values yielded near-identical results (r = 0.95 for reconstructed surfaces). Lastly, expression values were mean and variance normalized across parcels to facilitate visualization. Normalization had no effect on spatial correlation between gene expression and other variables since the spatial distribution of gene expression was left unaltered.

### Selection of brain-specific genes

Similar to 10Burt et al.2018; 25Fagerberg et al.2014; 36Genovese et al.2016, N = 2429 brain-specific genes were selected based on the criteria that expression in brain tissues were four times higher than the median expression across all tissue types, using Supplementary Dataset 1 of 25Fagerberg et al.2014. PC1 result shown in Figure 3A is computed from brain-specific genes, though findings are similar when using all genes (ρ = −0.56 with timescale map, Figure 3—figure supplement 1).

### Spatial statistics

All correlations between spatial maps (timescale, T1w/T2w, gene principal component [PC], and individual gene expressions) were computed using Spearman rank correlation. As noted in 11Burt et al.2020; 10Burt et al.2018; 87Vos de Wael et al.2020, neural variables vary smoothly and continuously across the cortical surface, violating the assumption of independent samples. As a result, when correlating two variables, each with nontrivial SA, the naive p-value is artificially lowered since it is compared against an inappropriate null hypothesis, i.e., randomly distributed or shuffled values across space. Instead, a more appropriate null hypothesis introduces SA-preserving null maps, which destroys any potential correlation between two maps while respecting their SAs. For all spatial correlation analyses, we generated N = 1000 null maps of one variable (timescale map unless otherwise noted), and the test statistic, Spearman correlation (ρ), is computed against the other variable of interest to build the null distribution. Two-tailed significance is then computed as the proportion of the null distribution that is less extreme than the empirical correlation value. All regression lines were computed by fitting a linear regression to log-timescale and the structural feature maps.

### SA modeling

To generate SA-preserving null maps, we used Moran Spectral Randomization (MSR) 92Wagner and Dray2015 from the python package BrainSpace 87Vos de Wael et al.2020. Details of the algorithm can be found in the above references. Briefly, MSR performs eigendecomposition on a spatial weight matrix of choice, which is taken here to be the inverse average geodesic distance matrix between all pairs of parcels in HCP-MMP1.0. The eigenvectors of the weight matrix are then used to generate randomized null feature maps that preserves the autocorrelation of the empirical map. We used the singleton procedure for null map generation. All significance values reported (Figures 2D and 3A–C) were adjusted using the above procedure.

We also compare two other methods of generating null maps: spatial variogram fitting (VF) 11Burt et al.2020 and spin permutation 1Alexander-Bloch et al.2018. Null maps were generated for timescale using spatial VF, while for spin permutation they were generated for vertex-wise T1w/T2w and gene PC1 maps before parcellation, so as to preserve surface locations of the parcellation itself. All methods perform similarly, producing comparable SA in the null maps, assessed using spatial variogram, as well as null distribution of spatial correlation coefficients between timescale and T1w/T2w (Figure 2—figure supplement 2).

### Principal component analysis (PCA) of gene expression

We used scikit-learn 72Pedregosa2011 PCA (sklearn.decomposition.PCA) to identify the dominant axes of gene expression variation across the entire AHBA dataset, as well as for brain-specific genes. PCA was computed on the variance-normalized average gene expression maps, X, an N × P matrix where N = 18,114 (or N = 2429 brain-specific) genes, and P = 180 cortical parcels. Briefly, PCA factorizes X such that X = USVT, where U and V are unitary matrices of dimensionality N × N and P × P, respectively. S is the same dimensionality as X and contains non-negative descending eigenvalues on its main diagonal (Λ). Columns of V are defined as the PCs, and the dominant axis of gene expression is then defined as the first column of V, whose proportion of variance explained in the data is the first element of Λ divided by the sum over Λ. Results for PC1 and PC2-10 are shown in Figure 3A and Figure 3—figure supplement 1, respectively.

### Comparison of timescale-transcriptomic association with single-cell timescale genes

Single-cell timescale genes were selected based on data from Table S3 of 81Tripathy et al.2017 and Online Table 1 of 5Bomkamp et al.2019. Using single-cell RNA sequencing data and patch-clamp recordings from transgenic mice cortical neurons, these studies identified genes whose expression significantly correlated with electrophysiological features derived from generalized linear integrate and fire (GLIF) model fits. We selected genes that were significantly correlated with membrane time constant (tau), input resistance (Rin or ri), or capacitance (Cm or cap) in the referenced data tables, and extracted the level of association between gene expression and those electrophysiological feature (correlation ‘DiscCorr’ in 81Tripathy et al.2017 and linear coefficient ‘beta_gene’ in 5Bomkamp et al.2019).

To compare timescale-gene expression association at the single-cell and macroscale level, we correlated the single-cell associations extracted above with the spatial correlation coefficient (macroscale ρ) between ECoG timescale and AHBA microarray expression data for those same genes, restricting to genes with p<0.05 for macroscale correlation (results identical for non-restrictive gene set). Overall association for all genes, as well as split by quintiles of their absolute macroscale correlation coefficient, are shown in Figure 3D. Example ‘single-cell timescale’ genes shown in Figure 3B and C are genes showing the highest correlations with those electrophysiology features reported in Table 2 of 5Bomkamp et al.2019.

### T1w/T2w-removed timescale and gene expression residual maps

To remove anatomical hierarchy as a potential mediating variable in timescale–gene expression relationships, we linearly regress out the T1w/T2w map from the (log) timescale map and individual gene expression maps. T1w/T2w was linearly fit to log-timescale, and the error between T1w/T2w-predicted timescale and empirical timescale was extracted (residual); this identical procedure was applied to every gene expression map to retrieve the gene residuals. SA-preserving null timescale residual maps were similarly created using MSR.

### PLS regression model

Due to multicollinearity in the high-dimensional gene expression dataset (many more genes than parcels), we fit a PLS model to the timescale map with one output dimension (sklearn.cross_decomposition.PLSRegression) to estimate regression coefficient for all genes simultaneously, resulting in N = 18,114 (or N = 2429 brain-specific) PLS weights 84Vértes et al.2016 101Whitaker et al.2016. To determine significantly associated (or ‘enriched’) genes, we repeated the above PLS-fitting procedure 1000 times but replaced the empirical timescale map (or residual map) with null timescale maps (or residual maps) that preserved its SA. Genes whose absolute empirical PLS weight was greater than 95% of its null weight distribution was deemed to be enriched, and submitted for GOEA.

### Gene ontology enrichment analysis

The Gene Ontology (GO) captures hierarchically structured relationships between GO items representing aspects of biological processes (BP), cellular components (CC), or molecular functions (MF). For example, ‘synaptic signaling’, ‘chemical synaptic transmission’, and ‘glutamatergic synaptic transmission’ are GO items with increasing specificity, with smaller subsets of genes associated with each function. Each GO item is annotated with a list of genes that have been linked to that particular process or function. GOEA examines the list of enriched genes from above to identify GO items that are more associated with those genes than expected by chance. We used GOATOOLS 58Klopfenstein et al.2018 to perform GOEA programmatically in python.

The list of unranked genes with significant empirical PLS weights was submitted for GOEA as the ‘study set’, while either the full ABHA list or brain-specific gene list was used as the ‘reference set’. The output of GOEA is a list of GO terms with annotated genes that are enriched or purified (i.e., preferentially appearing or missing in the study list, respectively) more often than by chance, determined by Fisher’s exact test.

Enrichment ratio is defined as follows: given a reference set with N total genes, and n were found to be significantly associated with timescale (in the study set), for a single GO item with B total genes annotated to it, where b of them overlap with the study set, then. Statistical significance is adjusted for multiple comparisons following Benjamini–Hochberg procedure (false discovery rate q-value reported in Figure 3F), and all significant GO items (q < 0.05) are reported in Figure 3F, in addition to some example items that did not pass significance threshold. For a detailed exposition, see 4Bauer et al.2017. Figure 3F shows results using brain-specific genes. The GO items that are significantly associated are similar when using the full gene set, but typically with larger q-values (Supplementary file 1 and 2) due to a much larger set of (non-brain-specific) genes. Control analysis was conducted using T1w/T2w, with 1000 similarly generated null maps, instead of timescale.

### Working memory ECoG data and analysis

The CRCNS fcx-2 and fcx-3 datasets include 17 intracranial ECoG recordings in total from epilepsy patients (10 and 7, respectively) performing the same visuospatial working memory task 54Johnson2019 53Johnson2018 51Johnson et al.2018 52Johnson et al.2018. Subject 3 (s3) from fcx-2 was discarded due to poor data quality upon examination of trial-averaged PSDs (high noise floor near 20 Hz), while s5 and s7 from fcx-3 correspond to s5 and s8 in fcx-2 and were thus combined. Together, data from 14 unique participants (22–50 years old, five females) were analyzed, with variable and overlapping coverage in PC (n = 14), PFC (n = 13), OFC (n = 8), and MTL (n = 9). Each channel was annotated as belonging to one of the above macro regions.

Experimental setup is described in 54Johnson2019; 53Johnson2018; 51Johnson et al.2018, 52Johnson et al.2018 in detail. Briefly, following a 1 s pre-trial fixation period (baseline), subjects were instructed to focus on one of two stimulus contexts (‘identity’ or ‘relation’ information). Then two shapes were presented in sequence for 200 ms each. After a 900 or 1150 ms jittered precue delay (delay1), the test cue appeared for 800 ms, followed by another post-cue delay period of the same length (delay2). Finally, the response period required participants to perform a 2-alternative forced choice test based on the test cue, which varied based on trial condition. For our analysis, we collapsed across the stimulus context conditions and compared neuronal timescales during the last 900 ms of baseline and delay periods from the epoched data, which were free of visual stimuli, in order to avoid stimulus-related event-related potential effects. Behavioral accuracy for each experimental condition was reported for each participant, and we average across both stimulus context conditions to produce a single working memory accuracy per participant.

Single-trial power spectra were computed for each channel as the squared magnitude of the Hamming-windowed Fourier Transform. We used 900 ms of data in all three periods (pre-trial, delay1, and delay2). Timescales were estimated by applying spectral parameterization as above, and the two delay-period estimates were averaged to produce a single delay period value. For comparison, we computed single-trial theta (3–8 Hz) and high-frequency activity (high gamma 67Mukamel et al.2005, 70–100 Hz) powers as the mean log-power within those frequency bins, as well as spectral exponent (χ). Single-trial timescale difference between delay and baseline was calculated as the difference of the log timescales due to the non-normal distribution of single-trial timescale estimates. All other neural features were computed by subtracting baseline from the delay period.

All neural features were then averaged across channels within the same regions, then trials, for each participant, to produce per-participant region-wise estimates, and finally averaged across all participants for the regional average in Figure 4B and C. Two-sided Mann–Whitney U-tests were used to test for significant differences in baseline timescale between pairs of regions (Figure 4B). Two-sided Wilcoxon rank-sum tests were used to determine the statistical significance of timescale change in each region (Figure 4C), where the null hypothesis was no change between baseline and delay periods (i.e., delay is 100% of baseline). Spearman rank correlation was used to determine the relationship between neural activity (timescale; theta; high-frequency; χ) change and working memory accuracy across participants (Figure 4D and Figure 4—figure supplement 1).

### Per-subject average cortical timescale across age

Since electrode coverage in the MNI-iEEG dataset is sparse and nonuniform across participants (Figure 2—figure supplement 1), simply averaging across parcels within individuals to estimate an average cortical timescale per participant confounds the effect of age with the spatial effect of cortical hierarchy. Therefore, we instead first normalize each parcel by its max value across all participants before averaging within participants, excluding those with fewer than 10 valid parcels (71 of 106 subjects remaining; results hold for a range of threshold values; Figure 4—figure supplement 2B). Spearman rank correlation was used to compute the association between age and average cortical timescale.

### Age–timescale association for individual parcels

Each cortical parcel had a variable number of participants with valid timescale estimates above the consistency threshold, so we compute Spearman correlation between age and timescale for each parcel, but including only those with at least five participants (114 of 180 parcels, result holds for a range of threshold values; Figure 4—figure supplement 2C). Spatial effect of age-timescale variation is plotted in Figure 4F, where parcels that did not meet the threshold criteria are grayed out. Mean age–timescale correlation from individual parcels was significantly negative under one-sample t-test.

### Data and materials’ availability

All data analyzed in this manuscript are from open data sources. All code used for all analyses and plots are publicly available on GitHub at https://github.com/rdgao/field-echos 35Gao2020 and https://github.com/rudyvdbrink/surface_projection 82van den Brink2020. See Tables 1 and 2 for details.

## References

1. 1. AFAlexander-Bloch
   2. HShou
   3. SLiu
   4. TDSatterthwaite
   5. DCGlahn
   6. RTShinohara
   7. SNVandekar
   8. ARaznahan
   2018On testing for spatial correspondence between maps of human brain structure and function178NeuroImage540551
   - doi10.1016/j.neuroimage.2018.05.070
   - pmid29860082
2. 1. DBadre
   2. MD'Esposito
   2009Is the rostro-caudal Axis of the frontal lobe hierarchical?10Nature Reviews Neuroscience659669
   - doi10.1038/nrn2667
   - pmid19672274
3. 1. CBaldassano
   2. JChen
   3. AZadbood
   4. JWPillow
   5. UHasson
   6. KANorman
   2017Discovering event structure in continuous narrative perception and memory95Neuron709721
   - doi10.1016/j.neuron.2017.06.041
   - pmid28772125
4. 1. SBauer
   2. CDessimoz
   3. NŠkunca
   2017Gene-Category Analysis175188Springer
   - doi10.1007/978-1-4939-3743-1
5. 1. CBomkamp
   2. SJTripathy
   3. CBengtssonGonzales
   4. JHjerling-Leffler
   5. AMCraig
   6. PPavlidis
   2019Transcriptomic correlates of electrophysiological and morphological diversity within and across excitatory and inhibitory neuron classes15PLOS Computational Biology
   - doi10.1371/journal.pcbi.1007113
   - pmid31211786
6. 1. MBreakspear
   2017Dynamic models of large-scale brain activity20Nature Neuroscience340352
   - doi10.1038/nn.4497
   - pmid28230845
7. 1. MBrett
   2. CJMarkiewicz
   3. MHanke
   4. M-ACôté
   5. BCipollini
   6. PMcCarthy
   7. DJarecka
   8. CPCheng
   9. YOHalchenko
   10. MCottaar
   11. SGhosh
   12. ELarson
   13. DWassermann
   14. SGerhard
   15. GRLee
   16. H-TWang
   17. EKastman
   18. JKaczmarzyk
   19. RGuidotti
   20. ODuek
   21. ARokem
   22. CMadison
   23. FCMorency
   24. BMoloney
   2020 Freecnipy/nibabel
8. 1. NBrunel
   2000Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons8Journal of Computational Neuroscience183208
   - doi10.1023/a:1008925309027
   - pmid10809012
9. 1. NBrunel
   2. XJWang
   2003What determines the frequency of fast network oscillations with irregular neural discharges? I. synaptic dynamics and excitation-inhibition balance90Journal of Neurophysiology415430
   - doi10.1152/jn.01095.2002
   - pmid12611969
10. 1. JBBurt
    2. MDemirtaş
    3. WJEckner
    4. NMNavejar
    5. JLJi
    6. WJMartin
    7. ABernacchia
    8. AAnticevic
    9. JDMurray
    2018Hierarchy of transcriptomic specialization across human cortex captured by structural neuroimaging topography21Nature Neuroscience12511259
    - doi10.1038/s41593-018-0195-0
    - pmid30082915
11. 1. JBBurt
    2. MHelmer
    3. MShinn
    4. AAnticevic
    5. JDMurray
    2020Generative modeling of brain maps with spatial autocorrelationbioRxiv
    - doi10.1101/2020.02.18.955054
12. 1. GBuzsáki
    2. CAAnastassiou
    3. CKoch
    2012The origin of extracellular fields and currents--EEG, ECoG, LFP and spikes13Nature Reviews Neuroscience407420
    - doi10.1038/nrn3241
    - pmid22595786
13. 1. RChaudhuri
    2. KKnoblauch
    3. MAGariel
    4. HKennedy
    5. XJWang
    2015A Large-Scale circuit mechanism for hierarchical dynamical processing in the primate cortex88Neuron419431
    - doi10.1016/j.neuron.2015.09.008
    - pmid26439530
14. 1. RChaudhuri
    2. HBj
    3. X-JWang
    2017Random recurrent networks near criticality capture the broadband power distribution of human ECoG dynamics8Cerebral Cortex113
    - doi10.1093/cercor/bhx233
15. 1. H-YSChien
    2. CJHoney
    2020Constructing and Forgetting Temporal Context in the Human Cerebral Cortex106Neuron675686
    - doi10.1016/j.neuron.2020.02.013
16. 1. SCole
    2. TDonoghue
    3. RGao
    4. BVoytek
    2019NeuroDSP: a package for neural digital signal processing4Journal of Open Source Software
    - doi10.21105/joss.01272
17. 1. SWDavis
    2. NADennis
    3. SMDaselaar
    4. MSFleck
    5. RCabeza
    2008Que PASA? The Posterior-Anterior Shift in Aging18Cerebral Cortex12011209
    - doi10.1093/cercor/bhm155
18. 1. EdeVillers-Sidani
    2. LAlzghoul
    3. XZhou
    4. KLSimpson
    5. RCSLin
    6. MMMerzenich
    2010Recovery of functional and structural age-related changes in the rat primary auditory cortex with operant training107PNAS1390013905
    - doi10.1073/pnas.1007885107
19. 1. MDemirtaş
    2. JBBurt
    3. MHelmer
    4. JLJi
    5. BDAdkinson
    6. MFGlasser
    7. DCVanEssen
    8. SNSotiropoulos
    9. AAnticevic
    10. JDMurray
    2019Hierarchical heterogeneity across human cortex shapes Large-Scale neural dynamics101Neuron11811194
    - doi10.1016/j.neuron.2019.01.017
    - pmid30744986
20. 1. TDonoghue
    2. MHaller
    3. EJPeterson
    4. PVarma
    5. PSebastian
    6. RGao
    7. TNoto
    8. AHLara
    9. JDWallis
    10. RTKnight
    11. AShestyuk
    12. BVoytek
    2020Parameterizing neural power spectra into periodic and aperiodic components23Nature Neuroscience16551665
    - doi10.1038/s41593-020-00744-x
    - pmid33230329
21. 1. NMDotson
    2. SJHoffman
    3. BGoodell
    4. CMGray
    2018Feature-Based Visual Short-Term Memory Is Widely Distributed and Hierarchically Organized99Neuron215226
    - doi10.1016/j.neuron.2018.05.026
22. 1. RDuarte
    2. ASeeholzer
    3. KZilles
    4. AMorrison
    2017Synaptic patterning and the timescales of cortical dynamics43Current Opinion in Neurobiology156165
    - doi10.1016/j.conb.2017.02.007
23. 1. RDuarte
    2. AMorrison
    2019Leveraging heterogeneity for neural computation with fading memory in layer 2/3 cortical microcircuits15PLOS Computational Biology
    - doi10.1371/journal.pcbi.1006781
24. 1. MDEyre
    2. MRenzi
    3. MFarrant
    4. ZNusser
    2012Setting the time course of inhibitory synaptic currents by mixing multiple GABA(A) receptor α subunit isoforms32Journal of Neuroscience58535867
    - doi10.1523/JNEUROSCI.6495-11.2012
    - pmid22539847
25. 1. LFagerberg
    2. BMHallström
    3. POksvold
    4. CKampf
    5. DDjureinovic
    6. JOdeberg
    7. MHabuka
    8. STahmasebpoor
    9. ADanielsson
    10. KEdlund
    11. AAsplund
    12. ESjöstedt
    13. ELundberg
    14. CA-KSzigyarto
    15. MSkogs
    16. JOTakanen
    17. HBerling
    18. HTegel
    19. JMulder
    20. PNilsson
    21. JMSchwenk
    22. CLindskog
    23. FDanielsson
    24. AMardinoglu
    25. ÅsaSivertsson
    26. KvonFeilitzen
    27. MForsberg
    28. MZwahlen
    29. IOlsson
    30. SNavani
    31. MHuss
    32. JNielsen
    33. FPonten
    34. MUhlén
    2014Analysis of the human Tissue-specific expression by Genome-wide integration of transcriptomics and Antibody-based proteomics13Molecular & Cellular Proteomics397406
    - doi10.1074/mcp.M113.035600
26. 1. DJFelleman
    2. DCVanEssen
    1991Distributed hierarchical processing in the primate cerebral cortex1Cerebral Cortex147
    - doi10.1093/cercor/1.1.1
    - pmid1822724
27. 1. ACFlint
    2. USMaisch
    3. JHWeishaupt
    4. ARKriegstein
    5. HMonyer
    1997NR2A Subunit Expression Shortens NMDA Receptor Synaptic Currents in Developing Neocortex17The Journal of Neuroscience24692476
    - doi10.1523/JNEUROSCI.17-07-02469.1997
28. 1. AFornito
    2. AArnatkevičiūtė
    3. BDFulcher
    2019Bridging the Gap between Connectome and Transcriptome23Trends in Cognitive Sciences3450
    - doi10.1016/j.tics.2018.10.005
29. 1. BFrauscher
    2. NvonEllenrieder
    3. RZelmann
    4. IDoležalová
    5. LMinotti
    6. AOlivier
    7. JHall
    8. DHoffmann
    9. DKNguyen
    10. PKahane
    11. FDubeau
    12. JGotman
    2018Atlas of the normal intracranial electroencephalogram: neurophysiological awake activity in different cortical areas141Brain11301144
    - doi10.1093/brain/awy035
30. 1. BFrauscher
    2. NvonEllenrieder
    3. RZelmann
    4. CRogers
    5. DKNguyen
    6. PKahane
    7. FDubeau
    8. JGotman
    2018High-Frequency Oscillations in the Normal Human Brain84Annals of Neurology374385
    - doi10.1002/ana.25304
31. 1. WJFreeman
    2000Mesoscopic neurodynamics: From neuron to brain94Journal of Physiology-Paris303322
    - doi10.1016/S0928-4257(00)01090-1
32. 1. WFreeman
    2. HErwin
    2008Freeman K-set3Scholarpedia
    - doi10.4249/scholarpedia.3238
33. 1. PGanupuru
    2. ABGoldring
    3. RHarun
    4. TDHanks
    2019Flexibility of Timescales of Evidence Evaluation for Decision Making29Current Biology20912097
    - doi10.1016/j.cub.2019.05.037
34. 1. RGao
    2. EJPeterson
    3. BVoytek
    2017Inferring synaptic excitation/inhibition balance from field potentials158NeuroImage7078
    - doi10.1016/j.neuroimage.2017.06.078
35. 1. RGao
    2020 Timescales Project Code
36. 1. GGenovese
    2. MFromer
    3. EAStahl
    4. DMRuderfer
    5. KChambert
    6. MLandén
    7. JLMoran
    8. SMPurcell
    9. PSklar
    10. PFSullivan
    11. CMHultman
    12. SAMcCarroll
    2016Increased burden of ultra-rare protein-altering variants among 4,877 individuals with schizophrenia19Nature Neuroscience14331441
    - doi10.1038/nn.4402
37. 1. JGjorgjieva
    2. GDrion
    3. EMarder
    2016Computational implications of biophysical diversity and multiple timescales in neurons and synapses for circuit performance37Current Opinion in Neurobiology4452
    - doi10.1016/j.conb.2015.12.008
38. 1. MFGlasser
    2. TSCoalson
    3. ECRobinson
    4. CDHacker
    5. JHarwell
    6. EYacoub
    7. KUgurbil
    8. JAndersson
    9. CFBeckmann
    10. MJenkinson
    11. SMSmith
    12. DCVanEssen
    2016A multi-modal parcellation of human cerebral cortex536Nature171178
    - doi10.1038/nature18933
39. 1. MFGlasser
    2. DCVanEssen
    2011Mapping Human Cortical Areas In Vivo Based on Myelin Content as Revealed by T1- and T2-Weighted MRI31Journal of Neuroscience1159711616
    - doi10.1523/JNEUROSCI.2180-11.2011
40. 1. JIGold
    2. MNShadlen
    2007The Neural Basis of Decision Making30Annual Review of Neuroscience535574
    - doi10.1146/annurev.neuro.29.051605.113038
41. 1. AGoulas
    2. J-PChangeux
    3. KWagstyl
    4. KAmunts
    5. NPalomero-Gallagher
    6. CCHilgetag
    2020The natural Axis of transmitter receptor distribution in the human cerebral cortexbioRxiv
    - doi10.1101/2020.09.28.316646
42. 1. GGryglewski
    2. RSeiger
    3. GMJames
    4. GMGodbersen
    5. AKomorowski
    6. JUnterholzner
    7. PMichenthaler
    8. AHahn
    9. WWadsak
    10. MMitterhauser
    11. SKasper
    12. RLanzenberger
    2018Spatial analysis and high resolution mapping of the human whole-brain transcriptome for integrative analysis in neuroimaging176NeuroImage259267
    - doi10.1016/j.neuroimage.2018.04.068
43. 1. MJHawrylycz
    2. ESLein
    3. ALGuillozet-Bongaarts
    4. EHShen
    5. LNg
    6. JAMiller
    7. LNvandeLagemaat
    8. KASmith
    9. AEbbert
    10. ZLRiley
    11. CAbajian
    12. CFBeckmann
    13. ABernard
    14. DBertagnolli
    15. AFBoe
    16. PMCartagena
    17. MMChakravarty
    18. MChapin
    19. JChong
    20. RADalley
    21. BDDaly
    22. CDang
    23. SDatta
    24. NDee
    25. TADolbeare
    26. VFaber
    27. DFeng
    28. DRFowler
    29. JGoldy
    30. BWGregor
    31. ZHaradon
    32. DRHaynor
    33. JGHohmann
    34. SHorvath
    35. REHoward
    36. AJeromin
    37. JMJochim
    38. MKinnunen
    39. CLau
    40. ETLazarz
    41. CLee
    42. TALemon
    43. LLi
    44. YLi
    45. JAMorris
    46. CCOverly
    47. PDParker
    48. SEParry
    49. MReding
    50. JJRoyall
    51. JSchulkin
    52. PASequeira
    53. CRSlaughterbeck
    54. SCSmith
    55. AJSodt
    56. SMSunkin
    57. BESwanson
    58. MPVawter
    59. DWilliams
    60. PWohnoutka
    61. HRZielke
    62. DHGeschwind
    63. PRHof
    64. SMSmith
    65. CKoch
    66. SGNGrant
    67. ARJones
    2012An anatomically comprehensive atlas of the adult human brain transcriptome489Nature391399
    - doi10.1038/nature11405
44. 1. MHawrylycz
    2. JAMiller
    3. VMenon
    4. DFeng
    5. TDolbeare
    6. ALGuillozet-Bongaarts
    7. AGJegga
    8. BJAronow
    9. C-KLee
    10. ABernard
    11. MFGlasser
    12. DLDierker
    13. JMenche
    14. ASzafer
    15. FCollman
    16. PGrange
    17. KABerman
    18. SMihalas
    19. ZYao
    20. LStewart
    21. A-LBarabási
    22. JSchulkin
    23. JPhillips
    24. LNg
    25. CDang
    26. DRHaynor
    27. AJones
    28. DCVanEssen
    29. CKoch
    30. ELein
    2015Canonical genetic signatures of the adult human brain18Nature Neuroscience18321844
    - doi10.1038/nn.4171
45. 1. BJHe
    2. JMZempel
    3. AZSnyder
    4. MERaichle
    2010The temporal structures and functional significance of scale-free brain activity66Neuron353369
    - doi10.1016/j.neuron.2010.04.020
    - pmid20471349
46. 1. CCHilgetag
    2. AGoulas
    2020‘Hierarchy’ in the organization of brain networks375Philosophical Transactions of the Royal Society B: Biological Sciences
    - doi10.1098/rstb.2019.0319
47. 1. CJHoney
    2. TThesen
    3. THDonner
    4. LJSilbert
    5. CECarlson
    6. ODevinsky
    7. WKDoyle
    8. NRubin
    9. DJHeeger
    10. UHasson
    2012Slow Cortical Dynamics and the Accumulation of Information over Long Timescales76Neuron423434
    - doi10.1016/j.neuron.2012.08.011
48. 1. CHuang
    2. BDoiron
    2017Once upon a (slow) time in the land of recurrent neuronal networks46Current Opinion in Neurobiology3138
    - doi10.1016/j.conb.2017.07.003
49. 1. JMHuntenburg
    2. P-LBazin
    3. DSMargulies
    2018Large-Scale Gradients in Human Cortical Organization22Trends in Cognitive Sciences2131
    - doi10.1016/j.tics.2017.11.002
50. 1. LIzhikevich
    2. RGao
    3. EPeterson
    4. BVoytek
    2018Measuring the average power of neural oscillationsbioRxiv
    - doi10.1101/441626
51. 1. ELJohnson
    2. JNAdams
    3. A-KSolbakk
    4. TEndestad
    5. PGLarsson
    6. JIvanovic
    7. TRMeling
    8. JJLin
    9. RTKnight
    2018Dynamic frontotemporal systems process space and time in working memory16PLOS Biology
    - doi10.1371/journal.pbio.2004274
52. 1. ELJohnson
    2. DKing-Stephens
    3. PBWeber
    4. KDLaxer
    5. JJLin
    6. RTKnight
    2018Spectral Imprints of Working Memory for Everyday Associations in the Frontoparietal Network12Frontiers in Systems Neuroscience
    - doi10.3389/fnsys.2018.00065
53. 1. EJohnson
    2018Intracranial EEG recordings of medial temporal, lateral frontal, and orbitofrontal regions in 10 human adults performing a visuospatial working memory task
    - doi10.6080/K06971S9
54. 1. EJohnson
    2019Intracranial EEG recordings of lateral frontal and parietal regions in 7 human adults performing a visuospatial working memory task
    - doi10.6080/K0VX0DQD
55. 1. RKanai
    2. GRees
    2011The structural basis of inter-individual differences in human behaviour and cognition12Nature Reviews Neuroscience231242
    - doi10.1038/nrn3000
    - pmid21407245
56. 1. AKhintchine
    1934Korrelationstheorie der stationären stochastischen Prozesse109Mathematische Annalen604615
    - doi10.1007/BF01449156
57. 1. SJKiebel
    2. JDaunizeau
    3. KJFriston
    2008A Hierarchy of Time-Scales and the Brain4PLOS Computational Biology
    - doi10.1371/journal.pcbi.1000209
58. 1. DVKlopfenstein
    2. LZhang
    3. BSPedersen
    4. FRamírez
    5. AWarwickVesztrocy
    6. ANaldi
    7. CJMungall
    8. JMYunes
    9. OBotvinnik
    10. MWeigel
    11. WDampier
    12. CDessimoz
    13. PFlick
    14. HTang
    2018GOATOOLS: A Python library for Gene Ontology analyses8Scientific Reports
    - doi10.1038/s41598-018-28948-z
59. 1. CKoch
    2. MRapp
    3. ISegev
    1996A Brief History of Time (Constants)6Cerebral Cortex93101
    - doi10.1093/cercor/6.2.93
60. 1. ESLein
    2. MJHawrylycz
    3. NAo
    4. MAyres
    5. ABensinger
    6. ABernard
    7. AFBoe
    8. MSBoguski
    9. KSBrockway
    10. EJByrnes
    11. LChen
    12. LChen
    13. T-MChen
    14. MChiChin
    15. JChong
    16. BECrook
    17. ACzaplinska
    18. CNDang
    19. SDatta
    20. NRDee
    21. ALDesaki
    22. TDesta
    23. EDiep
    24. TADolbeare
    25. MJDonelan
    26. H-WDong
    27. JGDougherty
    28. BJDuncan
    29. AJEbbert
    30. GEichele
    31. LKEstin
    32. CFaber
    33. BAFacer
    34. RFields
    35. SRFischer
    36. TPFliss
    37. CFrensley
    38. SNGates
    39. KJGlattfelder
    40. KRHalverson
    41. MRHart
    42. JGHohmann
    43. MPHowell
    44. DPJeung
    45. RAJohnson
    46. PTKarr
    47. RKawal
    48. JMKidney
    49. RHKnapik
    50. CLKuan
    51. JHLake
    52. ARLaramee
    53. KDLarsen
    54. CLau
    55. TALemon
    56. AJLiang
    57. YLiu
    58. LTLuong
    59. JMichaels
    60. JJMorgan
    61. RJMorgan
    62. MTMortrud
    63. NFMosqueda
    64. LLNg
    65. RNg
    66. GJOrta
    67. CCOverly
    68. THPak
    69. SEParry
    70. SDPathak
    71. OCPearson
    72. RBPuchalski
    73. ZLRiley
    74. HRRockett
    75. SARowland
    76. JJRoyall
    77. MJRuiz
    78. NRSarno
    79. KSchaffnit
    80. NVShapovalova
    81. TSivisay
    82. CRSlaughterbeck
    83. SCSmith
    84. KASmith
    85. BISmith
    86. AJSodt
    87. NNStewart
    88. K-RStumpf
    89. SMSunkin
    90. MSutram
    91. ATam
    92. CDTeemer
    93. CThaller
    94. CLThompson
    95. LRVarnam
    96. AVisel
    97. RMWhitlock
    98. PEWohnoutka
    99. CKWolkey
    100. VYWong
    101. MWood
    102. MBYaylaoglu
    103. RCYoung
    104. BLYoungstrom
    105. XFengYuan
    106. BZhang
    107. TAZwingman
    108. ARJones
    2007Genome-wide atlas of gene expression in the adult mouse brain445Nature168176
    - doi10.1038/nature05453
61. 1. YLerner
    2. CJHoney
    3. LJSilbert
    4. UHasson
    2011Topographic Mapping of a Hierarchy of Temporal Receptive Windows Using a Narrated Story31Journal of Neuroscience29062915
    - doi10.1523/JNEUROSCI.3684-10.2011
62. 1. YLiu
    2. ABeyer
    3. RAebersold
    2016On the Dependency of Cellular Protein Levels on mRNA Abundance165Cell535550
    - doi10.1016/j.cell.2016.03.014
63. 1. KMahjoory
    2. JMSchoffelen
    3. AKeitel
    4. JGross
    2020The frequency gradient of human resting-state brain oscillations follows cortical hierarchies9eLife
    - doi10.7554/eLife.53715
    - pmid32820722
64. 1. DSMargulies
    2. SSGhosh
    3. AGoulas
    4. MFalkiewicz
    5. JMHuntenburg
    6. GLangs
    7. GBezgin
    8. SBEickhoff
    9. FXCastellanos
    10. MPetrides
    11. EJefferies
    12. JSmallwood
    2016Situating the default-mode network along a principal gradient of macroscale cortical organization113PNAS1257412579
    - doi10.1073/pnas.1608282113
65. 1. KJMiller
    2. LBSorensen
    3. JGOjemann
    4. MdenNijs
    2009Power-law scaling in the brain surface electric potential5PLOS Computational Biology
    - doi10.1371/journal.pcbi.1000609
    - pmid20019800
66. 1. HMonyer
    2. NBurnashev
    3. DJLaurie
    4. BSakmann
    5. PHSeeburg
    1994Developmental and regional expression in the rat brain and functional properties of four NMDA receptors12Neuron529540
    - doi10.1016/0896-6273(94)90210-0
67. 1. RMukamel
    2. HGelbard
    3. AArieli
    4. UHasson
    5. IFried
    6. RMalach
    2005Coupling between neuronal firing, field potentials, and FMRI in human auditory cortex309Science951954
    - doi10.1126/science.1110913
    - pmid16081741
68. 1. JDMurray
    2. ABernacchia
    3. DJFreedman
    4. RRomo
    5. JDWallis
    6. XCai
    7. CPadoa-Schioppa
    8. TPasternak
    9. HSeo
    10. DLee
    11. X-JWang
    2014A hierarchy of intrinsic timescales across primate cortex17Nature Neuroscience16611663
    - doi10.1038/nn.3862
69. 1. YNagasaka
    2. KShimoda
    3. NFujii
    2011Multidimensional Recording (MDR) and Data Sharing: An Ecological Open Research and Educational Platform for Neuroscience6PLOS ONE
    - doi10.1371/journal.pone.0022561
70. 1. TOgawa
    2. HKomatsu
    2010Differential Temporal Storage Capacity in the Baseline Activity of Neurons in Macaque Frontal Eye Field and Area V4103Journal of Neurophysiology24332445
    - doi10.1152/jn.01066.2009
71. 1. OOssmy
    2. RMoran
    3. TPfeffer
    4. KTsetsos
    5. MUsher
    6. THDonner
    2013The Timescale of Perceptual Evidence Integration Can Be Adapted to the Environment23Current Biology981986
    - doi10.1016/j.cub.2013.04.039
72. 1. FPedregosa
    2011Scikit-learn: machine learning in Python12Journal of Machine Learning Research : JMLR28252830
73. 1. CMPegasiou
    2. AZolnourian
    3. DGomez-Nicola
    4. KDeinhardt
    5. JARNicoll
    6. AIAhmed
    7. GVajramani
    8. PGrundy
    9. MBVerhoog
    10. HDMansvelder
    11. VHPerry
    12. DBulters
    13. MVargas-Caballero
    2020Age-Dependent Changes in Synaptic NMDA Receptor Composition in Adult Human Cortical Neurons30Cerebral Cortex42464256
    - doi10.1093/cercor/bhaa052
74. 1. LPinto
    2. KRajan
    3. BDePasquale
    4. SYThiberge
    5. DWTank
    6. CDBrody
    2019Task-Dependent Changes in the Large-Scale Dynamics and Necessity of Cortical Regions104Neuron810824
    - doi10.1016/j.neuron.2019.08.025
75. 1. EPodvalny
    2. NNoy
    3. MHarel
    4. SBickel
    5. GChechik
    6. CESchroeder
    7. ADMehta
    8. MTsodyks
    9. RMalach
    2015A unifying principle underlying the extracellular field potential spectral responses in the human cortex114Journal of Neurophysiology505519
    - doi10.1152/jn.00943.2014
76. 1. CARunyan
    2. EPiasini
    3. SPanzeri
    4. CDHarvey
    2017Distinct timescales of population coding across cortex548Nature9296
    - doi10.1038/nature23020
77. 1. MSarafyazd
    2. MJazayeri
    2019Hierarchical reasoning by neural circuits in the frontal cortex364Science
    - doi10.1126/science.aav8911
    - pmid31097640
78. 1. NESavaskan
    2. AUBrauer
    3. RNitsch
    2004Molecular cloning and expression regulation of PRG-3, a new member of the plasticity-related gene family19European Journal of Neuroscience212220
    - doi10.1046/j.1460-9568.2003.03078.x
79. 1. JMShine
    2. MBreakspear
    3. PTBell
    4. KAEhgoetzMartens
    5. RShine
    6. OKoyejo
    7. OSporns
    8. RAPoldrack
    2019Human cognition involves the dynamic integration of neural activity and neuromodulatory systems22Nature Neuroscience289296
    - doi10.1038/s41593-018-0312-0
80. 1. BTeleńczuk
    2. NDehghani
    3. MLeVanQuyen
    4. SSCash
    5. EHalgren
    6. NGHatsopoulos
    7. ADestexhe
    2017Local field potentials primarily reflect inhibitory neuron activity in human and monkey cortex7Scientific Reports
    - doi10.1038/srep40211
81. 1. SJTripathy
    2. LToker
    3. BLi
    4. C-LCrichlow
    5. DTebaykin
    6. BOMancarci
    7. PPavlidis
    2017Transcriptomic correlates of neuron electrophysiological diversity13PLOS Computational Biology
    - doi10.1371/journal.pcbi.1005814
82. 1. RLvandenBrink
    2020 First Release of Surface Projection Code
83. 1. BvanVugt
    2. TvanKerkoerle
    3. DVartak
    4. PRRoelfsema
    2020The Contribution of AMPA and NMDA Receptors to Persistent Firing in the Dorsolateral Prefrontal Cortex in Working Memory40The Journal of Neuroscience24582470
    - doi10.1523/JNEUROSCI.2121-19.2020
84. 1. PEVértes
    2. TRittman
    3. KJWhitaker
    4. RRomero-Garcia
    5. FVáša
    6. MGKitzbichler
    7. KWagstyl
    8. PFonagy
    9. RJDolan
    10. PBJones
    11. IMGoodyer
    12. ETBullmore
    2016Gene transcription profiles associated with inter-modular hubs and connection distance in human functional magnetic resonance imaging networks371Philosophical Transactions of the Royal Society B: Biological Sciences
    - doi10.1098/rstb.2015.0362
85. 1. JVezoli
    2. LMagrou
    3. X-JWang
    4. KKnoblauch
    5. MVinck
    6. HKennedy
    2020Cortical hierarchy and the dual counterstream architecturebioRxiv
    - doi10.1101/2020.04.08.032706
86. 1. PVirtanen
    2. RGommers
    3. TEOliphant
    4. MHaberland
    5. TReddy
    6. DCournapeau
    7. EBurovski
    8. PPeterson
    9. WWeckesser
    10. JBright
    11. SJvanderWalt
    12. MBrett
    13. JWilson
    14. KJMillman
    15. NMayorov
    16. ARJNelson
    17. EJones
    18. RKern
    19. ELarson
    20. CJCarey
    21. İPolat
    22. YFeng
    23. EWMoore
    24. JVanderPlas
    25. DLaxalde
    26. JPerktold
    27. RCimrman
    28. IHenriksen
    29. EAQuintero
    30. CRHarris
    31. AMArchibald
    32. AHRibeiro
    33. FPedregosa
    34. PvanMulbregt
    35. SciPy 1.0 Contributors
    2020SciPy 1.0: fundamental algorithms for scientific computing in Python17Nature Methods261272
    - doi10.1038/s41592-019-0686-2
    - pmid32015543
87. 1. RVosdeWael
    2. OBenkarim
    3. CPaquola
    4. SLariviere
    5. JRoyer
    6. STavakol
    7. TXu
    8. SJHong
    9. GLangs
    10. SValk
    11. BMisic
    12. MMilham
    13. DMargulies
    14. JSmallwood
    15. BCBernhardt
    2020BrainSpace: a toolbox for the analysis of macroscale gradients in neuroimaging and connectomics datasets3Communications Biology
    - doi10.1038/s42003-020-0794-7
    - pmid32139786
88. 1. BVoytek
    2. ASKayser
    3. DBadre
    4. DFegen
    5. EFChang
    6. NECrone
    7. JParvizi
    8. RTKnight
    9. MD'Esposito
    2015Oscillatory dynamics coordinating human frontal networks in support of goal maintenance18Nature Neuroscience13181324
    - doi10.1038/nn.4071
89. 1. BVoytek
    2. MAKramer
    3. JCase
    4. KQLepage
    5. ZRTempesta
    6. RTKnight
    7. AGazzaley
    2015Age-Related Changes in 1/f Neural Electrophysiological Noise35Journal of Neuroscience1325713265
    - doi10.1523/JNEUROSCI.2332-14.2015
90. 1. BVoytek
    2. RTKnight
    2015Dynamic network communication as a unifying neural basis for cognition, development, aging, and disease77Biological Psychiatry10891097
    - doi10.1016/j.biopsych.2015.04.016
    - pmid26005114
91. 1. SVyas
    2. MDGolub
    3. DSussillo
    4. KVShenoy
    2020Computation Through Neural Population Dynamics43Annual Review of Neuroscience249275
    - doi10.1146/annurev-neuro-092619-094115
92. 1. HHWagner
    2. SDray
    2015Generating spatially constrained null models for irregularly spaced data using Moran spectral randomization methods6Methods in Ecology and Evolution11691178
    - doi10.1111/2041-210X.12407
93. 1. XJWang
    1999Synaptic basis of cortical persistent activity: the importance of NMDA receptors to working memory19The Journal of Neuroscience95879603
    - doi10.1523/JNEUROSCI.19-21-09587.1999
    - pmid10531461
94. 1. XJWang
    2002Probabilistic decision making by slow reverberation in cortical circuits36Neuron955968
    - doi10.1016/S0896-6273(02)01092-9
    - pmid12467598
95. 1. XJWang
    2008Decision making in recurrent neuronal circuits60Neuron215234
    - doi10.1016/j.neuron.2008.09.034
    - pmid18957215
96. 1. MWang
    2. NJGamo
    3. YYang
    4. LEJin
    5. X-JWang
    6. MLaubach
    7. JAMazer
    8. DLee
    9. AFTArnsten
    2011Neuronal basis of age-related working memory decline476Nature210213
    - doi10.1038/nature10243
97. 1. MWang
    2. YYang
    3. CJWang
    4. NJGamo
    5. LEJin
    6. JAMazer
    7. JHMorrison
    8. XJWang
    9. AFArnsten
    2013NMDA receptors subserve persistent neuronal firing during working memory in dorsolateral prefrontal cortex77Neuron736749
    - doi10.1016/j.neuron.2012.12.032
    - pmid23439125
98. 1. XJWang
    2020Macroscopic gradients of synaptic excitation and inhibition in the neocortex21Nature Reviews Neuroscience169178
    - doi10.1038/s41583-020-0262-x
    - pmid32029928
99. 1. DFWasmuht
    2. ESpaak
    3. TJBuschman
    4. EKMiller
    5. MGStokes
    2018Intrinsic neuronal dynamics predict distinct functional roles during working memory9Nature Communications
    - doi10.1038/s41467-018-05961-4
100. 1. TWatanabe
     2. GRees
     3. NMasuda
     2019Atypical intrinsic neural timescale in autism8eLife
     - doi10.7554/eLife.42256
     - pmid30717827
101. 1. KJWhitaker
     2. PEVértes
     3. RRomero-Garcia
     4. FVáša
     5. MMoutoussis
     6. GPrabhu
     7. NWeiskopf
     8. MFCallaghan
     9. KWagstyl
     10. TRittman
     11. RTait
     12. COoi
     13. JSuckling
     14. BInkster
     15. PFonagy
     16. RJDolan
     17. PBJones
     18. IMGoodyer
     19. ETBullmore
     2016Adolescence is associated with genomically patterned consolidation of the hubs of the human brain connectome113PNAS91059110
     - doi10.1073/pnas.1601745113
102. 1. KWimmer
     2. DQNykamp
     3. CConstantinidis
     4. ACompte
     2014Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory17Nature Neuroscience431439
     - doi10.1038/nn.3645
103. 1. TYanagawa
     2. ZCChao
     3. NHasegawa
     4. NFujii
     2013Large-Scale Information Flow in Conscious and Unconscious States: an ECoG Study in Monkeys8PLOS ONE
     - doi10.1371/journal.pone.0080845
104. 1. RZeraati
     2. TAEngel
     3. ALevina
     2020Estimation of autocorrelation timescales with approximate bayesian computationsbioRxiv
     - doi10.1101/2020.08.11.245944
105. 1. JZylberberg
     2. BWStrowbridge
     2017Mechanisms of Persistent Activity in Cortical Circuits: Possible Neural Substrates for Working Memory40Annual Review of Neuroscience603627
     - doi10.1146/annurev-neuro-070815-014006