Neural mechanisms of the EEG alpha-BOLD anticorrelation

Accepted Manuscript Neural mechanisms of the EEG alpha-BOLD anticorrelation J.C. Pang, P.A. Robinson PII:

S1053-8119(18)30646-3

DOI:

10.1016/j.neuroimage.2018.07.031

Reference:

YNIMG 15122

To appear in:

NeuroImage

Received Date: 28 February 2018 Revised Date:

2 July 2018

Accepted Date: 12 July 2018

Please cite this article as: Pang, J.C., Robinson, P.A., Neural mechanisms of the EEG alpha-BOLD anticorrelation, NeuroImage (2018), doi: 10.1016/j.neuroimage.2018.07.031. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

1

Neural mechanisms of the EEG alpha-BOLD anticorrelation

2

J.C. Panga,b,∗, P.A. Robinsona,b a

b

4

5

School of Physics, University of Sydney, New South Wales 2006, Australia Center for Integrative Brain Function, University of Sydney, New South Wales 2006, Australia

RI PT

3

Abstract

15

Keywords: EEG, fMRI, alpha rhythm, BOLD, neural field theory

16

1. Introduction

11 12 13

17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

M AN U

10

TE D

9

Scalp encephalography (EEG) is widely used to noninvasively record cortical activity aggregated at spatial resolutions of few centimeters on a millisecond time scale. Measurements using EEG have shown that spontaneous brain activity includes some resonant oscillations, one of the most dominant in normal adult humans being the 7.5–12 Hz alpha rhythm (Berger, 1929), which is especially prominent during a waking resting state with the eyes closed. Since its discovery, the alpha rhythm has been used as a marker of the level of arousal, vigilance, and attention (Strijkstra et al., 2003; Babiloni et al., 2004; Niedermeyer and Lopes da Silva, 2005; Thut et al., 2006; Rihs et al., 2007; Olbrich et al., 2009). Several works have also argued that it plays an important role in the perception of external stimuli (Babiloni et al., 2006; Thut et al., 2006; Haegens et al., 2011; H¨andel et al., 2011) and many cognitive abilities (Klimesch, 1999; Ba¸sar et al., 2001; Jensen et al., 2002; Sauseng et al., 2005; Zumer et al., 2014). Because EEG is spatially coarse-grained, many recent studies have used functional magnetic resonance imaging (fMRI), based on the blood oxygen level-dependent (BOLD) signal, to probe brain dynamics at finer spatial resolution (∼1 mm) but on a longer time scale

EP

8

AC C

7

SC

14

An experimentally tested neural field theory of the corticothalamic system is used to model brain activity and resulting experimental EEG data, and to elucidate the neural mechanisms and physiological basis of alpha-BOLD anticorrelation observed in concurrent EEG and fMRI measurements. Several studies have proposed that the anticorrelation originates from a causal link between changes in the alpha power and BOLD signal. However, the results in this study reveal that fluctuations in alpha and BOLD power do not generate one another but instead respectively result from high- and low-frequency components of the same underlying cortical activity, and that they are inversely correlated via variations in the strengths of corticothalamic and intrathalamic feedback, thereby explaining their anticorrelation.

6

∗

Corresponding author. Email address: [email protected] (J.C. Pang)

Preprint submitted to NeuroImage

July 13, 2018

ACCEPTED MANUSCRIPT

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74

RI PT

39 40

SC

38

M AN U

36 37

TE D

34 35

EP

33

(seconds) corresponding to the dynamics of spontaneous BOLD fluctuations. Notably, covariances of BOLD fluctuations are often used to define functional connectivity between discrete brain regions (Fox and Raichle, 2007). Hence, there is a strong interest in exploring the links between EEG and BOLD measurements to better understand structure–function relationships. Even though EEG and fMRI investigate brain function in different scales, their signals are tightly coupled and originate from the same underlying neural dynamics. In particular, EEG is recorded after volume conduction effects have removed high spatial frequencies of the neural activity (Nunez et al., 1994; Nunez and Srinivasan, 2006). In contrast, fMRI measures BOLD signal resulting from slow hemodynamic responses (typically in frequencies .0.2 Hz) that suppress high temporal frequencies of the neural activity (Zarahn et al., 1997; Robinson et al., 2006; He et al., 2008; Nir et al., 2008; Scheeringa et al., 2011; Aquino et al., 2012; Wang et al., 2014). Thus, they can be combined to complement one another and obtain data with high spatial and temporal resolutions [see for example the reviews of Ritter and Villringer (2006) and Laufs et al. (2008)]. The fusion of EEG and fMRI has led several studies to investigate the relationship of the alpha rhythm and the BOLD signal. These studies have shown that, in various brain structures, fluctuations in the BOLD signal are correlated or anticorrelated with fluctuations in the alpha-rhythm power (Goldman et al., 2002; Laufs et al., 2003; Moosmann et al., 2003; Gon¸calves et al., 2006; Danos et al., 2001; Feige et al., 2005). However, these correlations or anticorrelations are not consistently found across different studies, which is possibly due to high inter- and intra-subject variabilities (Gon¸calves et al., 2006), or methodological factors such as differences in recording paradigms and the manner in which the alpha-rhythm power was quantified from EEG recordings. Nonetheless, the most robust finding is for the occipital cortex where alpha is strongest and modulations of its power are negatively correlated with those of the BOLD signal. This concurs with several studies using other neuroimaging modalities that show similar relationships between fluctuations in alpha-rhythm power and those in blood flow and oxygenation signals (Jacquy et al., 1980; Sadato et al., 1998; Moosmann et al., 2003). Even with the abundance of statistical evidence relating the alpha rhythm and the BOLD signal, the underlying physiological mechanisms of the observed alpha-BOLD anticorrelations have not been established. Some studies have conjectured that a decrease in occipital BOLD signal is a result of ‘idling’ of cortex and decreased neural activity (Chawla et al., 1999; Wenzel et al., 2000; Mayhew et al., 2013), or alternatively linked to other functionally coupled processes such as vigilance (Moosmann et al., 2003) and wakefulness (Horovitz et al., 2008); however, mechanistic links remain elusive and concepts such as ‘idling’ remain poorly defined. Most studies have proposed that there is a direct causal link between fluctuations in the alpha power and the BOLD signal, as shown in Fig. 1a (Goldman et al., 2002; Laufs et al., 2003; Moosmann et al., 2003; Schirner et al., 2018). This was found by employing data processing steps, including high-pass filtering the EEG signal (usually >0.5 Hz), calculating the power spectrum from the squared amplitude of the filtered signals Fourier transform, and convolving the time course of the alpha power with an a priori hemodynamic response

AC C

32

2

ACCEPTED MANUSCRIPT

79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

RI PT

78

SC

77

M AN U

76

function (HRF) before comparing with the BOLD signal. These processing steps have the following shortcomings. First, high-pass filtering the EEG signal inevitably throws away the important frequencies that mainly drive the slow activity of BOLD. Second, it is well known that the HRF can take into account the slow processes mediating neural activity and BOLD; hence the neural signal itself, not its power, should be convolved with the HRF to predict the BOLD signal (Friston et al., 2000). This is what convolution of the alpha power with the HRF described above aims to achieve, but we argue that it does not follow the actual neural-HRF-BOLD relationship. Moreover, even though the convolution culminates to a low-frequency signal that appears to be comparable to BOLD, it represents the windowed low-frequency envelope of alpha, not the original low-frequency neural signal that drives BOLD. In order to overcome the abovementioned shortcomings and to not introduce confounds caused by ad hoc data processing, we instead directly extracted the low-frequency power of the actual EEG signal because it best represents the slow activity of BOLD. In this way, we can relate the powers of the low-frequency and alpha bands, which have identical dimensions. Thus, the type of analysis in this study is certainly different from what has been done previously, but it investigates the mechanisms underlying the alpha-BOLD anticorrelation in a new way that has a more compelling physical basis. In particular, we aim to show that fluctuations in the alpha power do not drive those in the BOLD power, nor vice versa. Instead, they respectively result from high- and low-frequency components of the same underlying cortical activity, which ensue from variations in the strengths of corticothalamic and intrathalamic feedback, as shown in Fig. 1b.

a

b

alpha

BOLD alpha

variations in corticothalamic & intrathalamic feedback

alpha BOLD

EP

BOLD

TE D

75

97 98 99 100 101 102 103 104 105

AC C

Figure 1: Possible causal links between fluctuations in EEG alpha and BOLD power. (a) Previous proposals that alpha power fluctuations directly influence BOLD fluctuations or vice versa. (b) Present proposal that variations in the strengths of corticothalamic and intrathalamic feedback simultaneously influence fluctuations in alpha and BOLD power.

Here, we use an experimentally tested and verified neural field theory of the corticothalamic system (Robinson et al., 2001, 2002, 2004, 2005) to model and analyze experimental EEG data recorded for a cohort of subjects in a relaxed waking eyes-closed state. First, properties of the low-frequency (combination of the <1 Hz slow-wave and the 1–4 Hz deltawave frequencies) and alpha bands of the EEG power spectra are analyzed to show that the anticorrelation between alpha and BOLD power (proxied by the low-frequency power) can be reproduced for a large group of subjects. Then, the model is used to fit the data and obtain time series of estimated corticothalamic gains within a single subject and show the validity of the causal relationship in Fig. 1b. 3

ACCEPTED MANUSCRIPT

113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135

136 137 138 139 140 141 142 143 144 145 146

RI PT

112

2.1. Experimental EEG data In this study, we reanalyze experimental EEG data comprising recordings of the Cz electrode from 2141 healthy subjects in a relaxed waking eyes-closed state used by Abeysuriya et al. (2015) originally obtained from the Brain Resource International Database. Subjects were 1049 females and 1092 males with a mean age of 31 years (SD = 20 years) and 27 years (SD = 19 years), respectively. The EEGs were acquired as part of a battery of electrophysiological tests. A QuickCap (Neuroscan) was used to obtain EEG data from 26 cephalic sites following the 10–20 international system. The EEG data were recorded relative to the average of A1 and A2 (mastoid) electrode sites. Horizontal eye movements were recorded with electrodes placed 1.5 cm lateral to the outer canthus of each eye, whereas vertical eye movements were recorded with electrodes placed 3 mm above the middle of the left eyebrow and 1.5 cm below the middle of the left bottom eye-lid. Skin resistance was kept at <5 kΩ. Scalp and electrooculogram (EOG) potentials were amplified and digitized continuously by a system (NuAmps, SCAN 4.3) having a frequency response from DC to 100 Hz (above which attenuating by 40 dB per decade) and a sampling rate of 500 Hz. EEG data were screened visually for artifacts, normal variants and changes in alertness. EOG correction was carried out using a technique based on Gratton et al. (1983). The EEG recordings are Fourier transformed to obtain the corresponding power spectrum for each subject. However, because brain states continuously evolve over time, taking the Fourier transform over a long period would mix features from a range of different states. Thus, to minimize this effect, the power spectra are calculated over 30 s epochs. Each power Rspectrum P (f ) is then normalized with respect to the total power, i.e., P (f ) = P (f )/ P (f ), in order to compare data across multiple subjects and remove unnecessary units of measurement.

SC

111

M AN U

110

TE D

109

In this section, we describe the experimental EEG data used in this study and the method to quantify their power spectral properties. We then briefly outline key details of the neural field model used to explain and further analyze the underlying neural dynamics. Finally, the method to fit the model to the data and the process to infer BOLD power are summarized.

EP

107 108

2. Theory and Methods

AC C

106

2.2. Calculation of EEG power spectral properties For each power spectrum, we quantify the spectral properties of the low-frequency and alpha bands, which are widely defined for a typical adult human to cover the frequency ranges of 0–4 Hz and 7.5–13 Hz, respectively; these bands are represented by the shaded regions in Fig. 2. Note that we are particularly interested in the .0.2 Hz frequencies that correspond to BOLD. However, data are limited for that range, thus we combined the commonly used <1 Hz slow-wave and 1–4 Hz delta-wave bands to construct the lowfrequency band. Nonetheless, this will not affect the premise of the study because power in the slow- and delta-wave bands are correlated (Benoit et al., 2000; Dang-Vu et al., 2008; He et al., 2008). The power spectral properties listed in Table 1 are calculated using the method of a previous study (Abeysuriya et al., 2015). 4

ACCEPTED MANUSCRIPT

100

RI PT

10-1

10-2

0.1

0.5

1

2

3

SC

10-3

5

10

20

40

M AN U

Figure 2: Illustrative experimental EEG power spectrum for a typical adult human subject in a relaxed waking eyes-closed state (Abeysuriya et al., 2015). The shaded regions represent the customary low-frequency (0–4 Hz) and alpha (7.5–13 Hz) bands. The filled circle shows the frequency and power of the alpha peak, while the hollow circle shows the frequency and power of the minimum closest to the lower bound of the alpha band (see Table 1 for details).

148 149 150 151 152 153 154 155 156 157 158

2.3. Neural field model We use an experimentally tested neural field model of the corticothalamic system to further understand the underlying dynamics of the experimental data described in Sec. 2.1. It involves averages over many neurons of neuronal dynamics such as firing rate and membrane potential, and has successfully predicted many features of brain function, including EEG time series and spectra, evoked response potentials, and prevailing arousal state (Rennie et al., 2002; Robinson et al., 2002, 2005; Abeysuriya et al., 2015). The key details of the model are summarized below and readers are advised to refer to previous articles for detailed discussions and derivations (Robinson et al., 2001, 2002, 2004, 2005; Abeysuriya et al., 2015). The model describes four neural populations whose connectivities are shown in Fig. 3. The populations are cortical excitatory (e) and inhibitory (i) neurons and thalamic reticular

AC C

147

EP

TE D

Table 1: Definitions of power spectral properties of the low-frequency and alpha bands used in this study.   Note that the points [fαmax , P (fαmax )] and fαmin , P (fαmin ) are marked in Fig. 2 for clarity. Power spectral property Definition Plow Total power in the low-frequency band Pα Total power in the alpha band fαmax Frequency of the peak in the alpha band fαmin Frequency of the minimum closest to the lower bound of the alpha band P (fαmax ) Corresponding power at fαmax P (fαmin) Corresponding power at fαmin max min Strength of the peak in the alpha band Rα = P (fα )/P (fα )

5

ACCEPTED MANUSCRIPT

159

(r) and relay (s) neurons; in addition, n denotes the external input.

φe νes excitatory (e) νei

Cortex

φi νis inhibitory (i) νie νii

φe

SC

φi

φs

RI PT

νee

νrs reticular (r) νre

Thalamus

M AN U

nucleus

φr

νsr

relay nuclei (s) νse νsn

φn

161

Each population a = e, i, r, s has a mean firing rate Qa (r, t) at position r and time t that depends on its mean soma potential Va (r, t) via a sigmoid function S such that

EP

160

TE D

Figure 3: Schematic diagram of the neural field model. The populations are cortical excitatory (e) and inhibitory (i) neurons and thalamic reticular (r) and relay (s) neurons. In addition, n denotes the external input. The quantity φa represents the output field of population a. On the other hand, νab quantifies the connection strength to population a from population b. Excitatory connections are shown with pointed arrowheads, while inhibitory connections are shown with rounded arrowheads.

Qmax , (1) 1 + exp {− [Va (r, t) − θ] /σ ′ } √ where θ is the mean threshold voltage, σ ′ π/ 3 is the standard deviation of the threshold distribution, and Qmax is the maximum firing rate. The mean membrane potential of population a is the sum of contributions of potentials Vab (r, t) from afferent populations b, with X Va (r, t) = Vab (r, t). (2)

162 163 164 165

AC C

Qa (r, t) = S[Va (r, t)] =

b

166 167

Each Vab responds to the afferent neural activity fields after accounting for synaptic, dendritic, and soma dynamics, with

6

ACCEPTED MANUSCRIPT

Dα (t)Vab (r, t) = νab φb (r, t − τab ),   1 1 d 1 d2 + + + 1, Dα (t) = αβ dt2 α β dt

171 172 173 174 175 176 177 178 179

RI PT

170

SC

169

Da (r, t)φa (r, t) = Qa (r, t), 1 ∂ 2 ∂ Da (r, t) = 2 2 + + 1 − ra2 ∇2 , γa ∂t γa ∂t 181 182 183 184 185 186

(6)

where γa = va /ra is the effective damping rate and ra is the characteristic axonal range. The model uses the local interaction approximation (Robinson et al., 1997), where only re is large enough to produce significant propagation effects, while other populations have localized activities with ra ≈ 0. In this work, we focus on the excitatory cortical field φe because excitatory pyramidal neurons have been shown to have a dominant contribution to EEG and the alpha rhythm (Lopes da Silva and Storm van Leeuwen, 1977). Setting all the time and space derivatives in Eqs (1) to (6) to zero yields spatially uniform (0) steady states of the system. The steady state value φe of φe is given by solutions of 
−1 (0) (0) S φe − (νee + νei ) φe =νes S νse φ(0) e   νrs  −1 (0)
(0) (0) + νsr S νre φe + S φe − (νee + νei ) φe νes  (0) + νsn φn , (7)

AC C

EP

187

(5)

TE D

180

(4)

where the connectivity strength νab is the product of the average number of synapses per neuron in population a from those in population b and their average synaptic strength, φa (r, t) is the output field of population a, τab is the time delay for signals to propagate to population a from population b, and α and β are the characteristic rates of soma voltage decay and rise, respectively. Note that the only nonzero discrete time delays in the model are τes = τis = τse = τre = t0 /2, where t0 is the corticothalamic loop delay time. When each population fires, its activity field propagate along the axons and axonal tree to provide incoming fields to other neural populations various distances away. Assuming a characteristic axonal propagation velocity va and an isotropic distribution of axons in the continuum limit, the outward propagation of the field obeys a damped wave equation generated by the source Qa (r, t) (Jirsa and Haken, 1996; Robinson et al., 1997, 2001, 2002) as

M AN U

168

(3)

(0)

188 189 190 191 192

where φn is the steady state component of the input stimulus. Small perturbations of the steady states allow us to linearize the equations and study dynamics in the linear regime (Robinson et al., 1997, 2001, 2002, 2003, 2004, 2005). Equations (2), (3), and (5) are already linear and will remain as is, whereas we can linearize Eq. (1) by retaining only the first term of its Taylor expansion, giving
  ′ (0) Qa (r, t) − Q(0) Va (r, t) − Va(0) , (8) a = S Va 7

ACCEPTED MANUSCRIPT

194

195 196 197

where S ′ is the derivative of the sigmoid function with respect to its argument. Thus, its linear approximation is Qa (r, t) = ρa Va (r, t), (9)   (0) with ρa = S ′ Va . By taking the Fourier transforms of Eqs (2)–(6) and (9), we are able to express the excitatory cortical field φe in terms of the external input φn in frequency space. The resulting transfer function is given by (Robinson et al., 2005)

RI PT

193

199 200 201 202 203 204 205 206

(11) (12)

where k is the spatial wave vector (with magnitude k = 2π/λ, where λ is the wavelength in m), ω is the angular frequency (ω = 2πf , where f is the usual frequency in Hz), L embodies the low-pass filter characteristics of synaptodendritic dynamics, and the gain Gab is the response of population a due to unit input from population b. The quantities Gese = Ges Gse , Gesre = Ges Gsr Gre , and Gsrs = Gsr Grs correspond to the overall gains for the corticalrelay nuclei, cortical-reticular-relay nuclei, and intrathalamic loops, respectively. The input stimulus φn is taken to be a white noise with |φn (k, ω)| = φn = constant for nonzero k and ω (Robinson et al., 1997, 2001, 2004, 2005). The parameters fixed in the model together with their nominal values are shown in Table 2.

TE D

198

(10)

M AN U

SC

Ges Gsn L2 eiωt0 /2 φe (k, ω) = , φn (k, ω) (1 − Gsrs L2 ) (1 − Gei L) (k 2 re2 + q 2 re2 )  2   iω [Gese L2 + Gesre L3 ] eiωt0 1 2 2 q re = 1 − Gee L + , − γe 1 − Gei L 1 − Gsrs L2 −1  −1  iω iω 1− , L(ω) = 1 − α β

AC C

EP

Table 2: Nominal model parameter values for a waking eyes-closed state (Rennie et al., 2002; Roberts and Robinson, 2012; Robinson et al., 2002, 2004; Abeysuriya et al., 2015). In each row, the columns list the parameter, its symbol, its value derived from physiology, and its unit. Quantity Symbol Value Unit Rate of synaptodendritic decay α 83 s−1 Rate of synaptodendritic rise β 769 s−1 Corticothalamic loop delay time t0 0.085 s Maximum neuron firing rate Qmax 340 s−1 Cortical damping rate γe 116 s−1 Excitatory axonal range re 0.086 m Threshold spread σ′ 3.8 mV Mean threshold voltage θ 12.9 mV Volume conduction filter cutoff k0 10 m−1 Cortical sheet length Lx 0.5 m Cortical sheet width Ly 0.5 m 207 208

Using Eq. (10) and the white noise approximation for φn , the EEG power spectrum P (ω) can be derived by integrating φe (k, ω) over all k as Z P (ω) = |φe (k, ω)|2F (k)d2 k, (13) 8

ACCEPTED MANUSCRIPT

212 213 214 215 216 217 218

219 220 221

222

RI PT

211

SC

210

where F (k) is an additional filter function to approximate the low-pass spatial filtering that arises due to volume conduction by the cerebrospinal fluid, skull, and scalp. Its functional 2 2 form is approximated as F (k) = e−k /k0 , where k0 is the low-pass filter cutoff. Thus, by fitting Eq. (13) to experimental EEG spectrum, relevant parameters, especially the gains Gab , can be estimated, giving us valuable insights into the underlying physiology. Further details of this fitting procedure are discussed in Sec. 2.4. Aside from the ability to estimate the EEG power spectrum, the model also has a reduced representation that has been previously used to parametrize brain states and stability at low frequencies. To derive the reduced model parameters, we first note that linear waves in the model obey the dispersion relation

1 − Gsrs L2 (1 − Gei L) k 2 re2 + q 2 re2 = 0. (14)

The main low-frequency dynamics (ω → 0) can be retained by writing L ≈ 1 in Eq. (14) except in the (1 − Gsrs L2 ) term, where the ω dependence is essential (Robinson et al., 2002). In this case, the dispersion relation becomes # " 2 
Gese + Gesre G iω ee 1 − Gsrs L2 − (1 − Gsrs ) eiωt0 = 0. − k 2 re2 + 1 − γe 1 − Gei (1 − Gsrs )(1 − Gei ) (15) By absorbing the gain terms within three quantities

M AN U

209

Gee , 1 − Gei Gese + Gesre Y = , (1 − Gsrs )(1 − Gei ) αβ , Z = −Gsrs (α + β)2

223

(16) (17) (18)

which parametrize the effective intracortical, corticothalamic, and intrathalamic loop gains, respectively (Robinson et al., 2002, 2005), Eq. (15) can be simplified to " #  2    iω (α + β)2 iωt0 (α + β)2 2 2 2 k re + 1 − e = 0. (19) L −Y 1+Z −X 1+Z γe αβ αβ

225 226 227 228 229 230 231 232 233 234 235

AC C

EP

224

TE D

X=

The solutions of the model are stable when all roots ω of the dispersion relations in Eqs (14) and (19) have negative imaginary part. On the other hand, instability occurs when Im(ω) changes sign. In terms of the fixed point’s eigenvalues, this means that an eigenvalue (for ω = 0) or complex-conjugate pair of eigenvalues (for ω > 0) crosses the imaginary axis, which correspond to saddle-node and Hopf bifurcations, respectively (Robinson et al., 2002, 2005; Roberts and Robinson, 2012). Previous works have shown that X, Y , and Z can accurately capture the main feedbacks that underlie brain states, including both wake and sleep stages (Robinson et al., 2002, 2005; Abeysuriya et al., 2015). It is important to note that X, Y , and Z can be obtained by estimating just the five gains Gee , Gei , Gese , Gesre , and Gsrs . We use these loop gains below to shed light on the physiological basis of cortical activity fluctuations. 9

ACCEPTED MANUSCRIPT

237 238 239 240 241

2.4. Fitting the experimental EEG power spectrum Equation (13) can be used to fit a power spectrum for a subject from the data described in Sec. 2.1. However, the model does not account for pericranial muscle artifacts that result in an electromyogram (EMG) signal that increases EEG power at high frequencies (>25 Hz) before fitting. This is allowed for by including an EMG component PEMG in the model power spectrum (Rowe et al., 2004): Pmodel (ω) = P (ω) + PEMG (ω), 2

248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272

SC

M AN U

247

TE D

245 246

(20) (21)

where P (ω) is given by Eq. (13), AEMG is a normalization factor, and fEMG is the characteristic EMG frequency with a value ranging from 10 to 50 Hz (Rowe et al., 2004; Abeysuriya and Robinson, 2016). Thus, instead of Eq. (13), Eq. (20) is used in fitting the data via the Metropolis-Hastings algorithm (Metropolis et al., 1953), which minimizes the squared difference between the model and data (Abeysuriya et al., 2015; Abeysuriya and Robinson, 2016). Example experimental and fitted power spectra for a subject in a relaxed waking eyes-closed state are shown in Fig. 4, which demonstrates the accuracy of the model. Most importantly, the model is able to explicitly capture the features of the low-frequency (0– 4 Hz) and alpha (7.5–13 Hz) bands, which are the bands of interest in this study. It also shows the continuity of the slow- and delta-wave bands, implying that their powers are correlated, as previous studies have found (Benoit et al., 2000; Dang-Vu et al., 2008; He et al., 2008), which is why we combined them into a single low-frequency band in our analysis. Note also that the shape of the spectrum is unaffected by the EMG correction, but there is a slight enhancement relative to the data at high frequencies. Moreover, limited resolution of experimental data in the f < 1 Hz frequency range causes an artificial discrepancy for single-subject fitting; however, this would improve for higher data resolution. Instead of applying the fitting algorithm to just a single power spectrum derived from EEG data recorded over a long duration (see Sec. 2.1), it is used to fit a sequence of spectra at successive times to track temporal changes in parameter values (Abeysuriya et al., 2015; Abeysuriya and Robinson, 2016). This is done by segmenting the data into 4-s overlapping step function windows that are moved in 1-s increments; the power spectrum is then estimated for each window. Equation (20) is fitted to each power spectrum, yielding time series of estimated model parameters, especially the five gains Gee , Gei , Gese , Gesre , and Gsrs . The resulting fits are found to be robust to noise and to different initial parameter values (Abeysuriya et al., 2015; Abeysuriya and Robinson, 2016). Using the fitted power spectra, we can robustly calculate the spectral properties of the low-frequency and alpha bands, as discussed in Sec. 2.2, vs. time. Similarly, the estimates of the five gains are used to calculate time series of the effective loop gains X, Y , and Z [via Eqs (16) to (18), respectively]. These calculations give us the time evolution of the power spectral properties and corticothalamic loop gains, allowing us to determine temporal correlations between them.

EP

243 244

[ω/ (2πfEMG )] PEMG (ω) = AEMG  2 , 1 + [ω/ (2πfEMG )]2

AC C

242

RI PT

236

10

ACCEPTED MANUSCRIPT

100

RI PT

10-1

10-2

10-4 0.1

0.5

1

2

3

5

SC

10-3

10

20

40

276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296

TE D

275

EP

274

2.5. Inferring BOLD power from EEG In previous EEG-fMRI studies, alpha-BOLD anticorrelations were found by explicitly comparing the time course of the BOLD signal with the peak alpha power time course convolved with an HRF [see for example the methodology of Goldman et al. (2002)]. This was done to take into account hemodynamic delays; however, it is only obtaining the lowfrequency envelope of alpha. Here, in order to not introduce confounds caused by ad hoc data processing, we instead want to investigate the observed anticorrelations by comparing the power of both alpha and BOLD. Even without available BOLD data, we can still infer its power from EEG data by using previous experimental and theoretical investigations that show that the BOLD signal for a typical adult human is slow and corresponds to low frequencies (Zarahn et al., 1997; Robinson et al., 2006; He et al., 2008; Nir et al., 2008; Scheeringa et al., 2011; Aquino et al., 2012; Wang et al., 2014). In fact, it is known that the BOLD signal is a result of the convolution of a neural signal and a HRF, which has low-pass filter properties and only retains significant power for frequencies .0.2 Hz, as a result of slow cerebrovascular responses including blood flow and oxygenation (Aquino et al., 2012, 2014; Pang et al., 2017). This is evident in Fig. 5 that shows significant BOLD HRF spectral content (value ≥5% of maximum) in the range f . 0.2 Hz, where the spectral content is calculated by taking the magnitude of the HRF’s Fourier transform (Aquino et al., 2014; Pang et al., 2016). Moreover, experiments have shown that fluctuations in the resting-state BOLD signal are positively correlated with those in the power of low-frequency EEG oscillations, defined in the 0–4 Hz (He et al., 2008), 1–4 Hz (Lu et al., 2007), or 2–4 Hz (Portnova et al., 2018) ranges, demonstrating that the low-frequency EEG power directly affects both the amplitude and power of the BOLD signal. Thus, these findings enable us to simply use the total power of the defined

AC C

273

M AN U

Figure 4: Sample experimental power spectra from a 30-s epoch EEG signal (solid line) and fitted power spectra with (dashed line) and without (dotted line) EMG correction for a subject in a relaxed waking eyes-closed state.

11

ACCEPTED MANUSCRIPT

297

0–4 Hz low-frequency EEG band as a proxy for BOLD power; i.e., PBOLD ∝ Plow .

RI PT

10-1

10-2

10-3 0.01

0.5

M AN U

0.1

SC

Spectral content of BOLD HRF (arbitrary units)

100

Figure 5: Spectral content of BOLD HRF for a typical adult human subject.

301 302 303 304 305 306

307 308 309 310 311 312 313 314 315 316 317 318 319

Here, we analyze the relaxed waking eyes-closed dataset discussed in Sec. 2.1 to test whether the anticorrelations between fluctuations in alpha and BOLD power can be obtained by analyses of EEG power spectral properties. Then, we fit the data with the neural field model, as discussed in Sec. 2.4, to obtain time series of the power spectral properties and corticothalamic loop gains of a single subject. Finally, temporal correlations between the mentioned quantities are calculated to substantiate our hypothesis that variations in the strengths of corticothalamic and intrathalamic feedback are inversely related to fluctuations in alpha and BOLD power, thereby explaining their anticorrelation.

TE D

300

EP

299

3. Results and Discussion

3.1. Alpha-BOLD anticorrelation Group analysis of the 2141 individual experimental spectra from Sec. 2.1 is performed by first calculating the power spectral properties of the low-frequency and alpha bands listed in Table 1. Only the 5th to 95th percentile values of each power spectral property are included in the analysis to remove confounding effects of unusual spectra, leaving us with 1927 data points. Then, correlations between the total power of the low-frequency band (Plow ) and the power spectral properties of the alpha band are calculated, as shown in Fig. 6. Figure 6 shows that Plow is negatively correlated with Pα , P (fαmax), and log(Rα ) (pvalues <10−14 ). This demonstrates that the low-frequency power decreases when the alpha power increases, and vice versa, regardless of the spectral property used to quantify alpha. What is then the significance of this anticorrelation? From our discussion in Sec. 2.5, the low-frequency EEG power can be used as a proxy for BOLD power. Thus, we can infer from this group analysis that, in general, the alpha power is negatively correlated with the

AC C

298

12

ACCEPTED MANUSCRIPT

a

b r = -0.67

0.4

c r = -0.54

0.4 0.3

0.3

0.2

0.2

0.2 0.1

0.1 0.2

0.3

0.4

0.5

0.6

RI PT

0.3

0.1

r = -0.42

0.4

0.1

0.2

0.3

0.4

0.5

0.2 0.4 0.6 0.8 1 1.2 1.4

324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345

3.2. Role of corticothalamic system in alpha-BOLD anticorrelation The results in Sec. 3.1 show that the alpha-BOLD anticorrelation is consistently found for a large cohort of subjects in a relaxed waking eyes-closed state; this relies on the premise that low-frequency EEG power corresponds to BOLD power. Here, we aim to provide evidence that the anticorrelation seen is a direct consequence of the dynamics of cortical activity fluctuations in frequency ranges that drive the respective signals, and that it is driven by changes of neural interactions in the corticothalamic system. We randomly choose one subject from the cohort and fit his/her EEG power spectra with the neural field model, as discussed in Sec. 2.4, resulting in a spectrogram in time-frequency space. The resulting experimental and fitted spectrograms are shown in Fig. 7. Figure 7 shows that features of the experimental data, such as the peaks and asymptotic low- and high-frequency behaviors, are reproduced by the model fit. The agreement is especially close in the range f . 14 Hz, which covers the low-frequency and alpha bands crucial for our study. Thus, the dynamics of the fitted model can be used to accurately interpret the underlying physiology of the experimental data. Using the fitted spectrogram, we quantify the power spectral properties listed in Table 1. Then, we use the estimated gains in producing the fitted spectra to infer effective gains of the intracortical, corticothalamic, and intrathalamic loops, as discussed in Sec. 2.4. These calculations yield time series of the power spectral properties and the feedback loop gains, as shown in Fig. 8. Figures 8a to 8d qualitatively show that Plow varies inversely relative to Pα , P (fαmax ), and log(Rα ), supporting the inverse relationship of low-frequency and alpha power demon-

TE D

323

EP

321 322

BOLD power. This shows that we can reproduce the alpha-BOLD anticorrelation found in the literature (Goldman et al., 2002; Laufs et al., 2003; Moosmann et al., 2003) by doing a simple analyses of the features of EEG power spectrum, which reflects the neural activity that drives both electrical and hemodynamic responses.

AC C

320

M AN U

SC

Figure 6: Correlations between the total power of the low-frequency band and the power spectral properties of the alpha band. For each panel, the dots represent the experimental data for each subject and the solid line represents a linear fit. The correlation coefficients (r) are also shown. (a) Plow vs. Pα . (b) Plow vs. P (fαmax ). (c) Plow vs. log(Rα ). Note that Plow , Pα , P (fαmax ), and Rα = P (fαmax )/P (fαmin ) are defined in Table 1 and that we used log(Rα ) instead of Rα because it is found to correlate better with Plow .

13

ACCEPTED MANUSCRIPT

a

b 10 0

40

10 0

30

10 -1

30

10 -1

20

10 -2

20

10

10 -3

10

10

20

30

40

10 -4

50

1

10

20

30

40

50

10 -2

10 -3

10 -4

SC

1

RI PT

40

a

0.35

M AN U

Figure 7: Spectrograms in time-frequency space. Note that the spectrograms are shown in the same logarithmic color bar for ease of comparison. (a) Experimental data. (b) Fit using the neural field model with EMG correction. 0.32

r = 0.97

0.12

0.11 0.93

r = 0.91

0.51 2.7

r = 0.67

EP

1.54 0.37

d

0.69

2.56

r = 0.64

AC C

1.85 1.15

1

10

20

30

40

f

0.31

0.57

0.26

0.45

0.21

TE D

0.72

c

0.5 0.48 0.46

0.22

0.23

b

e

0.67

g 0.8

0.51

0.74

0.36

0.67

1.87

h

0.12

1.27

0.09

0.67 50

0.07 1

10

20

30

40

50

Figure 8: Time series of power spectral properties and feedback loop gains. (a) Plow . (b) Pα . (c) P (fαmax ). (d) log(Rα ). (e) X. (f) Y . (g) X + Y . (h) Z. The black solid lines represent fit results, whereas the red dashed lines represent experimental data. The correlation coefficients (r) between experimental data and fit results are also shown.

346 347

strated by the group analysis in Fig. 6. Moreover, the results of the fit (black solid lines) generally follow the experimental data (red dashed lines), which is quantitatively supported 14

ACCEPTED MANUSCRIPT

349 350 351 352 353 354 355

by high correlation values (r), except for extreme peaks and troughs in P (fαmax) and log(Rα ). Meanwhile, Figs 8e to 8h show that X is almost constant, Y remains positive that is typical for waking states (Robinson et al., 2002; Abeysuriya et al., 2015), X + Y chiefly follows the profile of Y , and Z varies inversely relative to Y and X + Y . We also show the time course of X + Y in Fig. 8g because the dispersion relation in Eq. (19) has a slow-wave instability at X + Y = 1. The instability can be derived for a simple case when the state is spatially uniform (k = 0) and the intrathalamic loop is absent (Z = 0) (Robinson et al., 1997, 2002, 2005; Roberts and Robinson, 2012). This simplifies Eq. (19) to

362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386

M AN U

361

TE D

360

EP

358 359

(22)

which immediately shows that at ω = 0, the only solution is at X + Y = 1. The work of Roberts and Robinson (2012) generalized this derivation for k > 0 and Z > 0, and showed that the corticothalamic system always encounters a saddle-node bifurcation at X + Y = 1, resulting in a slow-wave instability and high levels of low-frequency power as this point is approached (Robinson et al., 1997, 2002, 2005). Thus, the value of X + Y , rather than X or Y individually, best informs us whether the current state of the corticothalamic system is approaching criticality for low frequencies. Overall, Fig. 8 demonstrates that as the low-frequency power increases and the alpha power decreases, X + Y increases and Z decreases. That is, high alpha power corresponds to low (X +Y ) and high Z, whereas high low-frequency power corresponds to high (X +Y ) and low Z. The reason for this is that proximity to the slow-wave stability boundary (X +Y = 1) enhances low-frequency power (Robinson et al., 2002). To make better sense of the abovementioned relationships, Fig. 9 shows snapshots of the power spectra and corresponding feedback loop gain combinations at times t = 10, 25, and 40 s. We then calculate correlations between the time courses of the spectral properties, X +Y , and Z to substantiate the qualitative findings discussed above. Figure 10 shows that Plow is negatively correlated with Pα , P (fαmax), and log(Rα ), which agrees with our results in Sec. 3.1 and experimental and empirical observations (Goldman et al., 2002; Laufs et al., 2003; Moosmann et al., 2003; Gon¸calves et al., 2006). Moreover, X + Y is positively correlated with Plow and negatively correlated with Pα , P (fαmax), and log(Rα ), whereas Z is negatively correlated with Plow and positively correlated with Pα , P (fαmax ), and log(Rα ). As discussed in Sec. 2.5, since the low-frequency power can serve as a proxy for BOLD power, our model results reproduce the observed anticorrelation of alpha and BOLD power, and show that this ensues from variations in the gains of corticothalamic and intrathalamic feedback loops. In particular, enhanced feedback between cortex and thalamus and decreased intrathalamic feedback simultaneously cause an increase in BOLD power and a decrease in alpha power, and vice versa. The correlations obtained in the within-subject analysis here are significantly higher than those in the group analysis in Sec. 3.1 because of the difference in timescales used to obtain the respective power spectra. The analysis in Sec. 3.1 captures averaged dynamics over a long timescale, which removes the strong coupling of spectral properties and inevitably induces a larger spread of spectral properties, whereas the analysis in this

AC C

356 357

SC

2  iω − X − Y eiωt0 = 0, 1− γe

RI PT

348

15

ACCEPTED MANUSCRIPT

100

3

RI PT

2 1

1

10-2

2

0.1 0.08

3

0.7

0.75

10-3 0.1

1

2

3

5

M AN U

0.5

0.8

SC

10-1

10

20

Figure 9: Snapshots of fitted power spectra at different times. Curves 1, 2, and 3 are power spectra at t = 10 s, 25 s, and 40 s, respectively. Similarly to Fig. 2, the shaded regions represent the customary lowfrequency (0–4 Hz) and alpha (7.5–13 Hz) bands. The inset shows the progression of (X +Y )-Z combinations from fits through time.

392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408

TE D

390 391

EP

389

section demonstrates continuous modulations of spectral properties over short timescales. Nonetheless, the results complement each other and together show the existence of alphaBOLD anticorrelations regardless of the timescale used. In summary, we have shown in this section the causal link between variations of the corticothalamic and intrathalamic feedback loops and fluctuations in the power of BOLD and alpha derived from the same cortical EEG activity. As a result, we find that the alphaBOLD anticorrelation is an emergent property of neural interactions in the corticothalamic system. Moreover, our findings are consistent across different randomly chosen subjects, as demonstrated in the Supplementary Material. We emphasize at this point that previous studies have proposed that there is a direct causal link between fluctuations in alpha power and BOLD signal (Goldman et al., 2002; Laufs et al., 2003; Moosmann et al., 2003; Schirner et al., 2018). These studies reached their findings by implementing several processing steps that have shortcomings, including filtering out low frequencies of the neural signal that that mainly drive the slow activity of BOLD and convolution with the HRF that culminates to a low-frequency signal that represent the windowed low-frequency envelope of alpha, not the original low-frequency neural signal that drives BOLD (see discussion in Secs 1 and 2.5). Hence, we instead directly used the lowfrequency EEG power and compared it with the alpha band power to understand the neural mechanisms of the alpha-BOLD anticorrelation. Despite using a different methodology, we draw support from the studies of Lu et al. (2007), He et al. (2008), and Portnova et al. (2018), who showed that the low-frequency EEG power directly affects both the amplitude and power of the BOLD signal, to directly link our results with previous experimental

AC C

387 388

16

ACCEPTED MANUSCRIPT

1 -1

-0.94 -23

)

0.93

-0.85 (8 10

-15

)

(1 10

-22

)

1

-0.93

0.98

0.96

(1 10 -22 )

(4 10 -37 )

(5 10 -27 )

-0.94

0.81

0.67

1

(7 10

-13

-7

1 0.79 -11 (1 10 )

)

(1 10 )

0.65

-0.46

-0.48

-0.46

(3 10 -7 )

(9 10 -4 )

(5 10 -4 )

(8 10 -4 )

1

SC

)

0.5

-0.65 -7

(3 10 )

M AN U

(9 10

-25

0

RI PT

(1 10

-0.5

1

1

Figure 10: Correlations of spectral properties and feedback loop gains Z and X + Y from fitted power spectra. The color indicates correlation values from −1 (blue) to 1 (red). The correlation values are also written in each box for clarity, with the numbers in parentheses representing p-values.

412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430

TE D

411

EP

410

studies. The novelty of this study is that it does not only demonstrate the existence of anticorrelation between alpha and BOLD but also provides mechanistic links to explain their relationship. Our results show that fluctuations in alpha and BOLD power do not drive one another, instead they are driven by high- and low-frequency components of the same underlying neural dynamics, and that they are inversely correlated via variations in the strengths of corticothalamic and intrathalamic feedback that underlie them, thereby explaining their anticorrelation. Significantly, we arrive at this conclusion by applying an existing neural field model, which already explains multiple phenomena of brain function (Robinson et al., 2001, 2002, 2004, 2005; Chiang et al., 2011; Abeysuriya et al., 2015) without the need to tailor it for this problem, thus avoiding overfitting. A recent modeling work by Schirner et al. (2018) used a neural mass model to explain the same alpha-BOLD anticorrelation. Similarly to what we have mentioned in Sec. 1, their analysis involved high-pass filtering the EEG signal at >1 Hz, which inevitably removes the low frequencies of neural activity that drive BOLD. Their main finding suggests that ongoing alpha power fluctuations are directly transformed into out of phase slow fMRI oscillations, which explains the empirical alpha-BOLD anticorrelations. We argue that this is not the case because we have shown that the emergence of alpha and fMRI oscillations are simultaneously derived from the fast and slow components of cortical activity, respectively, without the need to transform one to get the other. Then, they argued that fMRI oscillations result from the rectification of alpha waves by inhibitory populations that introduces a lowfrequency modulation of population firing rates and synaptic activity. This process involves

AC C

409

17

ACCEPTED MANUSCRIPT

435

4. Conclusions

440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471

SC

439

M AN U

437 438

We have analyzed the EEG spectra of a cohort of subjects in a relaxed waking eyes-closed state to explain the existence of the alpha-BOLD anticorrelation. First, we have used the low-frequency power of EEG as a proxy for BOLD power. This argument rests on the findings of previous studies showing that the BOLD signal corresponds to low frequencies (.0.2 Hz) (Zarahn et al., 1997; Robinson et al., 2006; He et al., 2008; Nir et al., 2008; Scheeringa et al., 2011; Aquino et al., 2012; Wang et al., 2014), as a result of slow cerebrovascular responses including blood flow and oxygenation (Aquino et al., 2012, 2014; Pang et al., 2017), and that fluctuations in the resting-state BOLD signal are positively correlated with those in the power of low-frequency EEG oscillations (.4 Hz) (Lu et al., 2007; He et al., 2008). Our analyses of the spectra of a large cohort of subjects have shown that robust alpha-BOLD anticorrelations are found across subjects, reproducing results of several EEG-fMRI studies (Goldman et al., 2002; Laufs et al., 2003; Moosmann et al., 2003). Second, we have fitted single-subject data with an experimentally tested and verified neural field model (Robinson et al., 2001, 2002, 2004, 2005) and have shown that the model can accurately reproduce the features of experimental power spectra, such as peaks, especially for frequencies .14 Hz that cover the low-frequency and alpha bands crucial for our study. This demonstrates that the model can provide a reliable interpretation of the underlying physiology of the experimental data. Finally, we have used the model fits to obtain time series of power spectral properties and feedback loop gains in the corticothalamic system of a single subject. Our measurements of correlations between the time series have shown that the combined effect of fluctuations in corticothalamic and intrathalamic loop gains can reproduce the observed anticorrelation between alpha and BOLD power, thus overcoming misconceptions that alpha and BOLD have a direct causal relationship in either direction. As a result, we have found that the alphaBOLD anticorrelation observed in experimental paradigms emerges from anticorrelations between high- and low-frequency neural activity caused by gain changes in the corticothalamic system. More importantly, the use of the physiologically based neural field model has allowed us to gain insights into the fundamental neural mechanisms of the alpha-BOLD anticorrelation, which goes beyond prior statistical and data-analysis methods. Moreover, our findings are consistent across different randomly chosen subjects, as demonstrated in the Supplementary Material. In the context of functional connectivity, recent EEG-fMRI studies have found that co-activating resting-state networks (RSNs) are correlated to the activity of specific EEG bands (Mantini et al., 2007; Jann et al., 2009; Scheeringa et al., 2012; Chang et al., 2013; Hutchison et al., 2015). These spatially dependent findings can be partially explained by the mechanisms explored in this work because we have shown that variations in the connectivity

TE D

436

EP

433

AC C

432

RI PT

434

nonlinearity much stronger than that seen in typical alpha dynamics. Finally, another aspect of their findings was that the anticorrelation ensues from modulations in cortical inhibitory feedback. We go beyond this to include the effects of changes in the strengths of corticothalamic and intrathalamic feedback to the alpha-BOLD anticorrelation.

431

18

ACCEPTED MANUSCRIPT

477

Appendix: Supplementary Material

473 474 475

RI PT

476

strengths of the corticothalamic system can account for the observed modulations in the lowfrequency and alpha power of cortical activity. In fact, Tu et al. (2015) already found this link, showing that changes in RSNs are influenced by changes in corticothalamic connections in schizophrenic patients. Further discussion of this topic is beyond the scope of the work but is an important direction to investigate further.

472

Supplementary Material that shows additional analyses for four other randomly chosen subjects, similar to those performed in Figs 7 to 10 of Sec. 3.2.

480

Conflicts of interest

482

The authors declare no conflicts of interest.

M AN U

481

SC

479

478

Acknowledgements

489

References

486 487

490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510

Abeysuriya, R. G., Rennie, C. J., Robinson, P. A., 2015. Physiologically based arousal state estimation and dynamics. J. Neurosci. Methods 253, 55–69. Abeysuriya, R. G., Robinson, P. A., 2016. Real-time automated EEG tracking of brain states using neural field theory. J. Neurosci. Methods 258, 28–45. Aquino, K. M., Robinson, P. A., Drysdale, P. M., 2014. Spatiotemporal hemodynamic response functions derived from physiology. J. Theor. Biol. 347, 118–136. Aquino, K. M., Schira, M. M., Robinson, P. A., Drysdale, P. M., Breakspear, M. J., 2012. Hemodynamic traveling waves in human visual cortex. PLoS Comp. Biol. 8 (3), e1002435. Babiloni, C., Miniussi, C., Babiloni, F., Carducci, F., Cincotti, F., Percio, C. D., Sirello, G., Fracassi, C., Nobre, A. C., Rossini, P. M., 2004. Sub-second “temporal attention” modulates alpha rhythms. A high-resolution EEG study. Cogn. Brain Res. 19 (3), 259–268. Babiloni, C., Vecchio, F., Bultrini, A., Romani, G. L., Rossini, P. M., 2006. Pre-and poststimulus alpha rhythms are related to conscious visual perception: A high-resolution EEG study. Cereb. Cortex 16 (12), 1690–1700. Ba¸sar, E., Ba¸sar-Eroglu, C., Karaka¸s, S., Sch¨ urmann, M., 2001. Gamma, alpha, delta, and theta oscillations govern cognitive processes. Int. J. Psychophysiol. 39 (2–3), 241–248. Benoit, O., Daurat, A., Prado, J., 2000. Slow (0.7–2 Hz) and fast (2–4 Hz) delta components are differently correlated to theta, alpha and beta frequency bands during NREM sleep. Clin. Neurophysiol. 111 (12), 2103–2106. ¨ Berger, H., 1929. Uber das elektrenkephalogramm des menschen. Eur. Arch. Psy. Clin. Neurosci. 87 (1), 527–570.

EP

485

AC C

484

TE D

488

We thank R.G. Abeysuriya for his assistance regarding data processing. We acknowledge the support of the Brain Resource International Database for EEG data acquisition and processing. This work was supported by the Australian Research Council Center of Excellence for Integrative Brain Function (ARC Center of Excellence grant CE140100007), the Australian Research Council Laureate Fellowship Grant FL140100025, and the Australian Research Council Discovery Project Grant DP170101778.

483

19

ACCEPTED MANUSCRIPT

520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561

RI PT

519

SC

518

M AN U

517

TE D

515 516

EP

513 514

Chang, C., Liu, Z., Chen, M. C., Liu, X., Duyn, J. H., 2013. EEG correlates of time-varying BOLD functional connectivity. NeuroImage 72, 227–236. Chawla, D., Lumer, E. D., Friston, K. J., 1999. The relationship between synchronization among neuronal populations and their mean activity levels. Neural Comput. 11 (6), 1389–1411. Chiang, A. K. I., Rennie, C. J., Robinson, P. A., van Albada, S. J., Kerr, C. C., 2011. Age trends and sex differences of alpha rhythms including split alpha peaks. Clin. Neurophysiol. 122 (8), 1505–1517. Dang-Vu, T. T., Schabus, M., Desseilles, M., Albouy, G., Boly, M., Darsaud, A., Gais, S., Rauchs, G., Sterpenich, V., Vandewalle, G., Carrier, J., Moonen, G., Balteau, E., Degueldre, C., Luxen, A., Phillips, C., Maquet, P., 2008. Spontaneous neural activity during human slow wave sleep. Proc. Natl. Acad. Sci. USA 105 (39), 15160–15165. Danos, P., Guich, S., Abel, L., Buchsbaum, M. S., 2001. EEG alpha rhythm and glucose metabolic rate in the thalamus in schizophrenia. Neuropsychobiology 43 (4), 265–272. Feige, B., Scheffler, K., Esposito, F., Salle, F. D., Hennig, J., Seifritz, E., 2005. Cortical and subcortical correlates of electroencephalographic alpha rhythm modulation. J. Neurophysiol. 93 (5), 2864–2872. Fox, M. D., Raichle, M. E., 2007. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nature Rev. Neurosci. 8, 700–711. Friston, K. J., Mechelli, A., Turner, R., Price, C. J., 2000. Nonlinear responses in fMRI: The Balloon model, Volterra kernels, and other hemodynamics. NeuroImage 12 (4), 466–477. Goldman, R. I., Stern, J. M., Engel Jr., J., Cohen, M. S., 2002. Simultaneous EEG and fMRI of the alpha rhythm. Neuroreport 13 (18), 2487–2492. Gon¸calves, S. I., de Munck, J. C., Pouwels, P. J. W., Schoonhoven, R., Kuijer, J. P. A., Maurits, N. M., Hoogduin, J. M., Van Someren, E. J. W., Heethaar, R. M., Lopes da Silva, F. H., 2006. Correlating the alpha rhythm to BOLD using simultaneous EEG/fMRI: Inter-subject variability. NeuroImage 30 (1), 203–213. Gratton, G., Coles, M. G. H., Donchin, E., 1983. A new method for off-line removal of ocular artifact. Electroencephalogr. Clin. Neurophysiol. 55 (4), 468–484. Haegens, S., N´ acher, V., Luna, R., Romo, R., Jensen, O., 2011. α-Oscillations in the monkey sensorimotor network influence discrimination performance by rhythmical inhibition of neuronal spiking. Proc. Natl. Acad. Sci. USA 108 (48), 19377–19382. H¨ andel, B. F., Haarmeier, T., Jensen, O., 2011. Alpha oscillations correlate with the successful inhibition of unattended stimuli. J. Cognitive Neurosci. 23 (9), 2494–2502. He, B. J., Snyder, A. Z., Zempel, J. M., Smyth, M. D., Raichle, M. E., 2008. Electrophysiological correlates of the brain’s intrinsic large-scale functional architecture. Proc. Natl. Acad. Sci. USA 105 (41), 16309–16044. Horovitz, S. G., Fukunaga, M., de Zwart, J. A., van Gelderen, P., Fulton, S. C., Balkin, T. J., Duyn, J. H., 2008. Low frequency BOLD fluctuation during resting wakefulness and light sleep: A simultaneous EEG-fMRI study. Hum. Brain Mapp. 29 (6), 671–682. Hutchison, R. M., Hashemi, N., Gati, J. S., Menon, R. S., Everling, S., 2015. Electrophysiological signatures of spontaneous BOLD fluctuations in macaque prefrontal cortex. NeuroImage 113, 257–267. Jacquy, J., Charles, P., Piraux, A., N¨ oel, G., 1980. Relationship between the electroencephalogram and the rheoencephalogram in the normal young adult. Neuropsychobiology 6, 341–348. Jann, K., Dierks, T., Boesch, C., Kottlow, M., Strik, W., Koenig, T., 2009. BOLD correlates of EEG alpha phase-locking and the fMRI default mode network. NeuroImage 45 (3), 903–916. Jensen, O., Gelfand, J., Kounios, J., Lisman, J. E., 2002. Oscillations in the alpha band (9–12 Hz) increase with memory load during retention in a short-term memory task. Cereb. Cortex 12 (8), 877–882. Jirsa, V. K., Haken, H., 1996. Field theory of electromagnetic brain activity. Phys. Rev. Lett. 77, 960. Klimesch, W., 1999. EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis. Brain Res. Rev. 29 (2–3), 169–195. Laufs, H., Daunizeau, J., Carmichael, D. W., Kleinschmidt, A., 2008. Recent advances in recording electrophysiological data simultaneously with magnetic resonance imaging. NeuroImage 40 (2), 515–528. Laufs, H., Kleinschmidt, A., Beyerle, A., Eger, E., Salek-Haddadi, A., Preibisch, C., Krakow, K., 2003. EEG-correlated fMRI of human alpha activity. NeuroImage 19 (4), 1463–1476.

AC C

511 512

20

ACCEPTED MANUSCRIPT

572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612

RI PT

570 571

SC

568 569

M AN U

566 567

TE D

564 565

Lopes da Silva, F. H., Storm van Leeuwen, W., 1977. The cortical source of the alpha rhythm. Neurosci. Lett. 6 (2–3), 237–241. Lu, H., Zuo, Y., Gu, H., Waltz, J. A., Zhan, W., Scholl, C. A., Rea, W., Yang, Y., Stein, E. A., 2007. Synchronized delta oscillations correlate with the resting-state functional MRI signal. Proc. Natl. Acad. Sci. USA 104 (46), 18265–18269. Mantini, D., Perrucci, M. G., Del Gratta, C., Romani, G. L., Corbetta, M., 2007. Electrophysiological signatures of resting state networks in the human brain. Proc. Natl. Acad. Sci. USA 104 (32), 13170– 13175. Mayhew, S. D., Ostwald, D., Porcaro, C., Bagshaw, A. P., 2013. Spontaneous EEG alpha oscillation interacts with positive and negative BOLD responses in the visual–auditory cortices and default-mode network. NeuroImage 76, 362–372. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., Teller, E., 1953. Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 1087–1092. Moosmann, M., Ritter, P., Krastel, I., Brink, A., Thees, S., Blankenburg, F., Taskin, B., Obrig, H., Villringer, A., 2003. Correlates of alpha rhythm in functional magnetic resonance imaging and near infrared spectroscopy. NeuroImage 20 (1), 145–158. Niedermeyer, E., Lopes da Silva, F. H., 2005. Electroencephalography: basic principles, clinical applications, and related fields. Lippincott Williams & Wilkins, Philadelphia, Pennsylvania, United States. Nir, Y., Mukamel, R., Dinstein, I., Privman, E., Harel, M., Fisch, L., Gelbard-Sagiv, H., Kipervasser, S., Andelman, F., Neufeld, M. Y., Kramer, U., Arieli, A., Fried, I., Malach, R., 2008. Interhemispheric correlations of slow spontaneous neuronal fluctuations revealed in human sensory cortex. Nature Neurosci. 11, 1100–1108. Nunez, P. L., Silberstein, R. B., Cadusch, P. J., Wijesinghe, R. S., Westdorp, A. F., Srinivasan, R., 1994. A theoretical and experimental study of high resolution EEG based on surface Laplacians and cortical imaging. Electroencephalogr. Clin. Neurophysiol. 90 (1), 40–57. Nunez, P. L., Srinivasan, R., 2006. Electric fields of the brain: The neurophysics of EEG, 2nd Edition. Oxford University Press, Oxford, United Kingdom. Olbrich, S., Mulert, C., Karch, S., Trenner, M., Leicht, G., Pogarell, O., Hegerl, U., 2009. EEG-vigilance and BOLD effect during simultaneous EEG/fMRI measurement. NeuroImage 45 (2), 319–332. Pang, J. C., Robinson, P. A., Aquino, K. M., 2016. Response-mode decomposition of spatio-temporal haemodynamics. J. R. Soc. Interface 13 (118), 20160253. Pang, J. C., Robinson, P. A., Aquino, K. M., Vasan, N., 2017. Effects of astrocytic dynamics on spatiotemporal hemodynamics: Modeling and enhanced data analysis. NeuroImage 147, 994–1005. Portnova, G. V., Tetereva, A., Balaev, V., Atanov, M., Skiteva, L., Ushakov, V., Ivanitsky, A., Martynova, O., 2018. Correlation of BOLD signal with linear and nonlinear patterns of EEG in resting state EEGinformed fMRI. Front. Hum. Neurosci. 11, 654. Rennie, C. J., Robinson, P. A., Wright, J. J., 2002. Unified neurophysical model of EEG spectra and evoked potentials. Biol. Cybern. 86 (6), 457–471. Rihs, T. A., Michel, C. M., Thut, G., 2007. Mechanisms of selective inhibition in visual spatial attention are indexed by α-band EEG synchronization. Eur. J. Neurosci. 25 (2), 603–610. Ritter, P., Villringer, A., 2006. Simultaneous EEG–fMRI. Neurosci. Biobehav. Rev. 30 (6), 823–838. Roberts, J. A., Robinson, P. A., 2012. Corticothalamic dynamics: Structure of parameter space, spectra, instabilities, and reduced model. Phys. Rev. E 85, 011910. Robinson, P. A., Drysdale, P. M., Van der Merwe, H., Kyriakou, E., Rigozzi, M. K., Germanoska, B., Rennie, C. J., 2006. BOLD responses to stimuli: Dependence on frequency, stimulus form, amplitude, and repetition rate. NeuroImage 31 (2), 585–599. Robinson, P. A., Rennie, C. J., Rowe, D. L., 2002. Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. Phys. Rev. E 65, 041924. Robinson, P. A., Rennie, C. J., Rowe, D. L., O’Connor, S. C., 2004. Estimation of multiscale neurophysiologic parameters by electroencephalographic means. Hum. Brain Mapp. 23, 53–72. Robinson, P. A., Rennie, C. J., Rowe, D. L., O’Connor, S. C., Gordon, E., 2005. Multiscale brain modelling.

EP

563

AC C

562

21

ACCEPTED MANUSCRIPT

623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654

RI PT

621 622

SC

620

M AN U

618 619

TE D

616 617

EP

615

Philos. Trans. R. Soc. B 360 (1457), 1043–1050. Robinson, P. A., Rennie, C. J., Wright, J. J., 1997. Propagation and stability of waves of electrical activity in the cerebral cortex. Phys. Rev. E 56, 826–840. Robinson, P. A., Rennie, C. J., Wright, J. J., Bahramali, H., Gordon, E., Rowe, D. L., 2001. Prediction of electroencephalographic spectra from neurophysiology. Phys. Rev. E 63, 021903. Robinson, P. A., Whitehouse, R. W., Rennie, C. J., 2003. Nonuniform corticothalamic continuum model of electroencephalographic spectra with application to split-alpha peaks. Phys. Rev. E 68, 21922. Rowe, D. L., Robinson, P. A., Harris, A. W., Felmingham, K. L., Lazzaro, I. L., Gordon, E., 2004. Neurophysiologically-based mean-field modelling of tonic cortical activity in post-traumatic stress disorder (PTSD), chronic schizophrenia, first episode schizophrenia and attention deficit hyperactivity disorder (ADHD). J. Integr. Neurosci. 3, 453–487. Sadato, N., Nakamura, S., Oohashi, T., Nishina, E., Fuwamoto, Y., Waki, A., Yonekura, Y., 1998. Neural networks for generation and suppression of alpha rhythm: a PET study. Neuroreport 9 (5), 893–897. Sauseng, P., Klimesch, W., Doppelmayr, M., Pecherstorfer, T., Freunberger, R., Hanslmayr, S., 2005. EEG alpha synchronization and functional coupling during top-down processing in a working memory task. Hum. Brain Mapp. 26 (2), 148–155. Scheeringa, R., Fries, P., Petersson, K.-M., Oostenveld, R., Grothe, I., Norris, D. G., Hagoort, P., Bastiaansen, M. C. M., 2011. Neuronal dynamics underlying high- and low-frequency EEG oscillations contribute independently to the human BOLD signal. Neuron 69 (3), 572–583. Scheeringa, R., Petersson, K. M., Kleinschmidt, A., Jensen, O., Bastiaansen, M. C. M., 2012. EEG alpha power modulation of fMRI resting-state connectivity. Brain Connect. 2 (5), 254–264. Schirner, M., McIntosh, A. R., Jirsa, V., Deco, G., Ritter, P., 2018. Inferring multi-scale neural mechanisms with brain network modelling. eLife 7, e28927. Strijkstra, A. M., Beersma, D. G. M., Drayer, B., Halbesma, N., Daan, S., 2003. Subjective sleepiness correlates negatively with global alpha (8–12 Hz) and positively with central frontal theta (4–8 Hz) frequencies in the human resting awake electroencephalogram. Neurosci. Lett. 340 (1), 17–20. Thut, G., Nietzel, A., Brandt, S. A., Pascual-Leone, A., 2006. α-band electroencephalographic activity over occipital cortex indexes visuospatial attention bias and predicts visual target detection. J. Neurosci. 26 (37), 9494–9502. Tu, P.-C., Lee, Y.-C., Chen, Y.-S., Hsu, J.-W., Li, C.-T., Su, T.-P., 2015. Network-specific cortico-thalamic dysconnection in schizophrenia revealed by intrinsic functional connectivity analyses. Schizophr. Res. 166 (1–3), 137–143. Wang, Y.-F., Liu, F., Long, Z.-L., Duan, X.-J., Cui, Q., Yan, J. H., Chen, H.-F., 2014. Steady-state BOLD response modulates low frequency neural oscillations. Sci. Rep. 4, 7376. Wenzel, R., Wobst, P., Heekeren, H. H., Kwong, K. K., Brandt, S. A., Kohl, M., Obrig, H., Dirnagl, U., Villringer, A., 2000. Saccadic suppression induces focal hypooxygenation in the occipital cortex. J. Cerebr. Blood Flow Metab. 20 (7), 1103–1110. Zarahn, E., Aguirre, G. K., D’Esposito, M., 1997. Empirical analyses of BOLD fMRI statistics. NeuroImage 5, 178–197. Zumer, J. M., Scheeringa, R., Schoffelen, J.-M., Norris, D. G., Jensen, O., 2014. Occipital alpha activity during stimulus processing gates the information flow to object-selective cortex. PLoS Biol. 12 (10), e1001965.

AC C

613 614

22

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

a

b alpha

BOLD

BOLD

alpha

variations in corticothalamic & intrathalamic feedback

a

B

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

-0.94

-0.5

0

0.5

1

1

AC C

(1 10 -23 )

-1

EP

1

-0.85

0.93

(8 10 -15 )

(1 10 -22 )

-0.93

0.98

0.96

(1 10 -22 )

(4 10 -37 )

-27 (5 10 )

-0.94

0.81

0.67

0.79

(9 10 -25 )

(7 10 -13 )

(1 10 -7 )

(1 10 -11 )

0.65

-0.46

-0.48

-0.46

-0.65

(3 10 -7 )

(9 10 -4 )

(5 10 -4 )

(8 10 -4 )

(3 10 -7 )

1 1 1 1

EP

10-1

AC C

100

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

10-2

10-3

0.1

0.5

1

2

3

5

10

20

4

TE D

φe

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

νee

Cortex

EP

νes excitatory (e) νei

φi

φs

Thalamus

AC C

νis inhibitory (i) νie νii

φi

νrs reticular (r) νre

nucleus

φr νsr

relay nuclei (s) νse νsn

φn

φe

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

10-2

AC C

10-1

EP

100

10-3

10-4 0.1

0.5

1

2

3

5

10

20

40

Spectral content of BOLD HRF (arbitrary units)

ACCEPTED MANUSCRIPT

AC C

EP

0.1

TE D

M AN U

SC

RI PT

100

10-1

10-2

10-3

0.01 0.5

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

a

b r = -0.67

0.4

c r = -0.54

0.4

0.4

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0.2

0.3

0.4

0.5

0.6

0.1

0.2

0.3

0.4

0.5

0.2 0.4

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

b

AC C

a 40

10 0

40

30

10 -1

30

20

10 -2

20

10

10 -3

10

1

10

20

30

40

50

10 -4

1

10

20

0.35

r = 0.97

0.23 0.11 0.93

EP

b

c

d

2.7

e 0.5 0.48 0.46

0.22 0.12 r = 0.91

AC C

0.72 0.51

0.32

TE D

a

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

r = 0.67

0.69

f

0.31

0.57

0.26

0.45

0.21

0.67

g 0.8

1.54

0.51

0.74

0.37

0.36

0.67

2.56

1.87

r = 0.64

h

0.12

1.85

1.27

0.09

1.15

0.67 50

0.07

1

10

20

30

40

1

10

20

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

EP

100

10

-1

10

-2

AC C

3 2 1

1

0.1

2

0.08

3

0.7 10-3 0.1

0.75 0.5

0.8 1

2

3

5

10

20