Relationships Between Neuronal Oscillatory Amplitude and Dynamic Functional Connectivity (original) (raw)

Journal Article

Prejaas Tewarie ,

Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK

Address correspondence to Prejaas Tewarie, Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, University Park, Nottingham, UK. Email: [email protected].

Search for other works by this author on:

Benjamin A E Hunt ,

Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK

Search for other works by this author on:

George C O’Neill ,

Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK

Search for other works by this author on:

Aine Byrne ,

School of Mathematical Sciences, University of Nottingham, Nottingham, UK

Search for other works by this author on:

Kevin Aquino ,

Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK

Search for other works by this author on:

Markus Bauer ,

School of Psychology, University of Nottingham, University Park, Nottingham, UK

Search for other works by this author on:

Karen J Mullinger ,

Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK

Search for other works by this author on:

Stephen Coombes ,

School of Mathematical Sciences, University of Nottingham, Nottingham, UK

Search for other works by this author on:

Matthew J Brookes

Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK

Search for other works by this author on:

Stephen Coombes and Matthew J. Brookes contributed equally.

Author Notes

Received:

22 October 2017

Revision received:

12 May 2018

Cite

Prejaas Tewarie, Benjamin A E Hunt, George C O’Neill, Aine Byrne, Kevin Aquino, Markus Bauer, Karen J Mullinger, Stephen Coombes, Matthew J Brookes, Relationships Between Neuronal Oscillatory Amplitude and Dynamic Functional Connectivity, Cerebral Cortex, Volume 29, Issue 6, June 2019, Pages 2668–2681, https://doi.org/10.1093/cercor/bhy136
Close

Navbar Search Filter Mobile Enter search term Search

Abstract

Event-related fluctuations of neural oscillatory amplitude are reported widely in the context of cognitive processing and are typically interpreted as a marker of brain “activity”. However, the precise nature of these effects remains unclear; in particular, whether such fluctuations reflect local dynamics, integration between regions, or both, is unknown. Here, using magnetoencephalography, we show that movement induced oscillatory modulation is associated with transient connectivity between sensorimotor regions. Further, in resting-state data, we demonstrate a significant association between oscillatory modulation and dynamic connectivity. A confound with such empirical measurements is that increased amplitude necessarily means increased signal-to-noise ratio (SNR): this means that the question of whether amplitude and connectivity are genuinely coupled, or whether increased connectivity is observed purely due to increased SNR is unanswered. Here, we counter this problem by analogy with computational models which show that, in the presence of global network coupling and local multistability, the link between oscillatory modulation and long-range connectivity is a natural consequence of neural networks. Our results provide evidence for the notion that connectivity is mediated by neural oscillations, and suggest that time–frequency spectrograms are not merely a description of local synchrony but also reflect fluctuations in long-range connectivity.

Introduction

Neural oscillations reflect rhythmic electrophysiological activity, synchronized across large neuronal assemblies in the mammalian brain. Detectable using either invasive (electrocorticography) or noninvasive (electroencephalographic (EEG)/magnetoencephalographic (MEG)) methods, these oscillations dominate electrophysiological recordings and have been assessed in a vast number of neuroscientific studies. Robust modulation of oscillatory amplitudes by sensory or cognitive tasks, neuromodulation, pharmacological intervention, or even hypercapnia has been well characterized (Fries et al. 2001; Ward 2003; Hall et al. 2011) and is typically taken as being reflective of neural “activity”. Importantly, the perturbation of those oscillatory signatures has been shown in illnesses ranging from developmental disorders to severe psychoses and neurodegeneration (Uhlhaas and Singer 2006). The importance of studying neural oscillations is therefore well-established. However, despite over a century of research, the precise role of spectral perturbations remains poorly understood. Here, we probe the role of oscillations in dynamic coordination of spatially separate brain areas.

Neuronal oscillations can be characterized by their frequency, phase, and amplitude. The amplitude of low frequency (alpha/beta) oscillations has been shown to decrease on stimulation, and increase on stimulus cessation, leading to a hypothesis that low-frequency amplitude increase is a marker of inhibition (Cassim et al. 2001). This is supported by studies suggesting that poststimulus increases are higher in individuals with a high concentration of the inhibitory neurotransmitter, gamma-aminobutyric acid (GABA) (Gaetz et al. 2011). Pharmacological intervention studies support this, since modulation of the GABAergic system generates robust changes in beta amplitude (Hall et al. 2011; Muthukumaraswamy and Singh 2013). By contrast, high frequency (gamma) oscillations increase with excitatory input to a population (e.g., stimulation (Muthukumaraswamy and Singh 2013)) or pharmacological manipulation. However, recruitment of inhibitory interneurons is also necessary for the emergence of gamma oscillations and GABAergic agonizts boost stimulus-induced gamma (Thiele et al. 2012). This suggests that the gamma rhythm is a result of a delicate balance of excitatory and inhibitory input. Understanding the neurochemical basis of oscillations represents an important step in their characterization. However, it does not necessarily add to our understanding of their role in supporting cognitive function.

In recent years, both theoretical and neurobiological studies have shown how neural oscillations are coordinated across brain regions (Nolte et al. 2004; Fries 2005). This coordination has been probed using metrics that assess similarity in amplitude envelopes and phase synchronization between areas. Amplitude correlations have revealed spatial signatures of several resting-state networks, (Brookes et al. 2011; Hipp et al. 2012) and have also allowed tracking of dynamic changes in connectivity between regions. For example, O’Neill et al. (2017) showed that increases in amplitude envelope correlations within a sensorimotor network were present around the time of movement. Though instructive, these studies provide little insight on the origins of the fluctuations driving amplitude envelopes. In contrast, more distinct hypotheses exist concerning the role of phase synchronization in interareal communication (Fries 2005). Specifically, the “communication via coherence” hypothesis suggests that phase synchrony between 2 regions provides optimum windows of high electrical potential which facilitate action potentials, and therefore the passage of information between areas. For resting-state data, phase synchronization metrics have been applied in the context of static connectivity, where a link between oscillatory amplitude and static phase-based connectivity has been previously shown (Daffertshofer and van Wijk 2011; Moon et al. 2015). Specifically, regions with large amplitude oscillations tend to exhibit stronger static connectivity. However, this fails to provide an understanding on the role of dynamic amplitude modulations often measured in EEG/MEG spectral analysis (Ward 2003). In principle, more insight can be gained by the study of dynamic connectivity (Baker et al. 2014; Brookes et al. 2014; O’neill et al. 2015; Abeysuriya et al. 2018). For example, Vidaurre et al. (2016) used dynamic connectivity to show that phase coherence between left and right motor cortex increased during a post-movement increase in beta oscillations (Vidaurre et al. 2016). This observation makes it tempting to speculate that dynamic increases in the amplitude of (e.g., beta) oscillations in a specific brain area might be driven by an increased level of functional connectivity between that area and the rest of the brain. Indeed, this would resonate with the communication by coherence hypothesis; for example, if one considers a single (hub) region, which is connected to many others in the same network, if the other regions transiently synchronize and coherent action potentials arrive synchronously at the hub, the likelihood is an increase in phase locked rhythmic dendritic currents, and therefore an increase in oscillatory amplitude within the hub.

In this paper, we test a hypothesis that temporal modulations in oscillatory amplitude relate to modulations in connectivity. We first acquire MEG data during a motor task and reconstruct brain activity within multiple regions (or parcels) defined based upon the AAL atlas. We then show that long-range connectivity (by which we mean connectivity between AAL regions) from the sensorimotor cortex to the rest of the brain mirrors oscillatory dynamics, providing evidence for a link between connectivity and amplitude fluctuations. Second, we generalize this to spontaneous brain activity measured in the resting state, showing that during periods of high oscillatory amplitude, functional connectivity in established networks is higher than average. In order to keep track of the fast dynamics of amplitude envelopes, a high temporal resolution phase connectivity measure is required and we employ a technique introduced by Breakspear et al. (2004). This allows for maximum temporal resolution in connectivity assessment (Thatcher et al. 2009; Lackner et al. 2014). Unfortunately, empirical measurements (i.e., MEG) are confounded by the fact that increased amplitude necessarily means increased signal-to-noise ratio (SNR): this means that the question of whether amplitude and connectivity are genuinely coupled, or whether increased connectivity is observed purely due to increased SNR remains unanswered. To counter this problem, we combine our analysis of empirical MEG data with biophysically informed computational models to test whether the observed relationship between local oscillatory amplitude and global phase connectivity dynamics can be explained using models of neuronal dynamics.

Methods

Three datasets were acquired at the Sir Peter Mansfield Imaging Centre, University of Nottingham (2 datasets are reported here—the third dataset is reported in SI).

For all cohorts, the studies were approved by the University of Nottingham Medical School Ethics Committee, and all subjects gave written informed consent prior to participation.

MEG Data Acquisition and Preprocessing

MEG data were acquired using the third order synthetic gradiometer configuration of a 275-channel CTF MEG system (MISL, Coquitlam, Canada), at a sampling rate of 600 Hz and using a 150-Hz low pass antialiasing filter. Three coils were attached to the participant’s head as fiducial markers at the nasion, left and right preauricular points. These coils were energized continuously throughout acquisition to allow localization of the fiducial points (and thus the head) relative to the geometry of the MEG sensor array. Before MEG acquisition, the surface of the participant’s head was digitized relative to the 3 coils, using a 3D digitizer (Polhemus Inc., Vermont). Subsequent surface matching of the digitized head shape to an equivalent head shape extracted from an anatomical magnetic resonance image (acquired using a 3 T or 7 T MRI scanner running an MPRAGE sequence with a resolution of 1 mm3) allowed coregistration of brain anatomy to MEG sensor geometry. Following collection, MEG data were inspected for artefacts and trials containing excessive interference were removed. Trials in which the participant was found to have moved more than 7 mm from their starting position were also removed.

MEG Source Localization

An atlas-based beamforming approach was adopted to project MEG sensor level data into source-space (Hillebrand et al. 2012). The cortex was parcellated into 78 individual regions according to the automated anatomical labeling (AAL) atlas (Tzourio-Mazoyer et al. 2002). This was done by registering each subject’s anatomical MR image to an MNI template and labeling all cortical voxels according to the 78 cortical ROIs (Gong et al. 2009). Subsequently, an inverse registration to anatomical subject space was performed. A beamformer (Robinson and Vrba 1999) was employed to generate a single signal representative of electrophysiological activity within each of these AAL regions. To achieve this, for each region, first the centre of mass was derived. Voxels were then defined on a regular 4 mm grid covering the entire region and the beamformer estimated timecourse of electrical activity was derived for each voxel. To generate a single timecourse representing the whole region, individual voxel signals were combined in a weighted sum, with the weights defined as a Gaussian function of distance from the centre of mass (Brookes et al. 2016; Tewarie, Hillebrand et al. 2016). To calculate the individual timecourse for every voxel, a scalar beamformer was used (Robinson and Vrba 1999). Covariance was computed within a 1–150 Hz frequency window and a time window spanning the whole experiment (excluding bad trials) (Brookes et al. 2008). Regularization was applied to the data covariance matrix using the Tikhonov method with a regularization parameter equal to 5% of the maximum eigenvalue of the unregularised covariance matrix. The forward model was based upon a dipole approximation (Sarvas 1987) and a multiple local sphere head model (Huang et al. 1999) fitted to the MRI scalp surface extracted from the co-registered MRI. Dipole orientation was determined using a nonlinear search for optimum “signal-to-noise ratio” (SNR; Robinson and Vrba 1999) where SNR here is defined in terms of the pseudo _Z_-statistic SNR=(wTCw/(wTΣw))⁠; w corresponds to the beamformer weights, C to the covariance matrix, and the noise covariance matrix which was assumed to be a scaled identity matrix. This complete process resulted in 78 electrophysiological timecourses, each representative of a separate AAL region. This approach was applied to each subject individually.

Dynamic Connectivity and Amplitude Envelope Analysis

Following beamformer estimation of regional timecourses and after frequency filtering, leakage correction was applied to reduce the effect of spurious connectivity due to the ill-posed MEG inverse problem. Since leakage manifests as a zero-time lag linear summation of underlying signals, orthogonalisation of beamformer projected timecourses results in the reduction of leakage, albeit at the expense of genuine zero-lag connectivity. We applied a symmetric multivariate leakage correction as described in Colclough et al. (2015). Symmetric orthogonalisation was applied to timecourses from all 78 regions after data had been filtered in 8 frequency bands of interest (delta 1–4 Hz, theta 4–8 Hz, alpha 8–13 Hz, beta 13–30 Hz, gamma1 30–48 Hz, gamma2 48–70 Hz, gamma3 70–100 Hz, and gamma4 100–130 Hz). Following leakage correction, we computed the amplitude envelopes of the data. The amplitude envelopes were computed individually for each AAL region by first calculating the Hilbert transform of each filtered time signal. Subsequently, the modulus of the analytic signal Aˆi(t) is computed to obtain the amplitude envelopes Ei(t) of each time signal.

Dynamic connectivity was estimated using the phase-difference derivative (PDD) (Breakspear et al. 2004); a measure that captures the stability of phase relationships between 2 filtered unaveraged timecourses. We extracted the instantaneous phase of each time signal as φi(t)=atan(Im[Aˆi(t)]/Re[Aˆi(t)])⁠. The unwrapped phases were used for subsequent analysis. For each pair of regional signals i and j, we assumed that there is phase locking if their phase difference is approximately constant. If the phase difference is approximately constant, the derivative of the phase difference is approximately zero. We, therefore, choose a small cut-off near zero, such that phase locking is assumed if

|dΔφij(t)dt|≈0→|dΔφij(t)dt|<k

(1)

Here Δφij refers to the instantaneous phase difference between region i and j, |X| refers to the absolute value of X, and k to the cut-off. To define the “phase-difference derivative (pdd)”, we can then define a binary representation of dynamic phase connectivity for each pair of regions by (Breakspear et al. 2004)

pddij(t)={1if|dΔφij(t)dt|≤k0if|dΔφij(t)dt|>k,

(2)

pdd was computed for every time point, but individual time samples were only set to one if |dΔφij(t)/dt|<k for a period longer than one oscillation (i.e., we assume phase locking for less than one cycle is meaningless). This reduces the effect of noise fluctuations in the phase difference between 2 signals. This binary connectivity timecourse indicates periods of phase locking. The cut-off value for k depends on the sample frequency, the frequency band of interest and the signal-to-noise ratio of the data (Breakspear et al. 2004). The choice for k was data-driven and determined by the area under the curve in the left tail of the distribution of dΔφij(t)/dt for every subject individually. To analyze the effect of the cut-off value on our results, we chose a range of cut-off values (0.05, 0.1, 0.2, 0.3, and 0.4) (see Supplementary Figs S1 and S2). Afterwards, visual inspection was carried out to evaluate whether these cut-off values were including any noise in the data. For each region, we obtained the average connectivity from that region to all other regions as

PDDi(t)=1N∑j=1,j≠iNpddij(t),

(3)

that is, PDD_i_(t) is representative average connectivity between region i and all other regions (N = number of regions). PDD_i_(t) was calculated for each region, frequency band, and subject. The same method was applied to computational models (see below). Note that the pdd captures information about the stability of phase differences, which is the same information underlying the traditional phase locking value (Lachaux et al. 1999), although in a dynamic framework.

Analysis of Task-Based MEG Data

For the self-paced motor task, we selected a 30-s window around every button press (15 s before and 15 s after the press) to represent a single trial. Data were filtered into the beta band (13–30 Hz) due to the well-known movement induced response in this frequency range. Amplitude envelopes were baseline corrected and averaged over trials (Fig. 2B for the sensory cortex). For every AAL region, PDD_i_(t) timecourses were constructed (eq. 2) (independently for each trial) and averaged across all trials. We expected to see a relationship between amplitude envelope and PDD, in brain regions of interest.

The beta response to finger movement is well characterized, by an event-related desynchronization (loss in amplitude during movement) and a rebound (increase above baseline on movement cessation). To see a distinction in these phases, we subdivided every 30 s trial into windows: before the button press (−10 s to −5 s), during beta desynchronization (−3 s to 0.5 s), during the beta rebound (1.5 –3.5 s), and after the button press (5 –10 s). Both neural oscillatory amplitude and PDD timecourses were averaged within every window in order to track local beta dynamics and network formation. Two statistical tests were performed:

Analysis of Resting-State MEG Data

For each subject, amplitude envelopes and PDD were computed for all regions and all frequency bands. The last sample of the amplitude envelope data was discarded such that the number of data points for the PDD and the amplitude envelopes were equal. Pearson correlation coefficients C(PDD(f),E(f)) between the amplitude envelope and PDD were then computed for all regions, independently for each frequency band (f). This resulted in distributions (across subjects) of C(PDD(f),E(f)) values for every region and band. We then performed significance testing to analyze whether these correlation coefficients could be obtained by chance from a null-distribution Cnull(PDD(f),E(f))⁠, which were based on the 2 methods of surrogate data reconstruction. First, surrogate datasets for resting-state data were reconstructed by splitting every source projected time series at a random sample in 2 and swapping the 2 parts of the time series (Stam 2005). This procedure ensures that the spectrum, amplitude distribution, and nonstationarity of the data is exactly the same as the original data, but destroys phase locking and (amplitude coupling) in the original data. Second, the uniform phase randomization, as used above, was applied to resting-state data for each subject. For every subject, ten surrogate datasets were obtained. Then for every frequency band and region, the distribution of the empirical C(PDD(f),E(f)) values, consisting of _N_sub values, were compared with the null-distribution using a nonparametric Wilcoxon rank-sum test. The null-distribution consisted of _N_subx10 values. Performing this analysis for every frequency band resulted in _N_regions x _F P_values (where F is the number of frequency bands F = 5). These P values were corrected for multiple comparisons with the false discovery rate (correction for _N_regions x F tests) (Benjamini and Hochberg 1995). This procedure was repeated for every threshold k (eq. 2; see Supplementary Figs S1 and S2).

Networks of Large-Scale Neuronal Models

The above analyses test for an empirical relationship between oscillatory amplitude within a single region, and connectivity between that region and all other regions; significance in any of the statistical tests outlined above could imply such a relationship. However, an inescapable confound in empirical data is that of SNR. As amplitude increases, so does SNR. This raises the question of whether a significant relationship between connectivity and amplitude timecourses arises due to a genuine relationship between the 2 (e.g., increased connectivity “causes” increased amplitude) or whether, during periods of high SNR, we are simple better able to “see” connectivity. This, we believe, is a fundamental limitation of empirical data which cannot (at present) be solved using statistical testing via surrogate data. With this problem in mind, we employed large-scale biophysically informed neuronal models to further test for a relationship between connectivity and amplitude dynamics.

We employed 3 different large-scale models, known to mimic the neural oscillations measured using MEG. Specifically, the Wilson–Cowan model (Wilson and Cowan 1972), the Liley model (Liley et al. 2002), and the Robinson model (Robinson et al. 2001, 2002). The first 2 models include excitatory and inhibitory populations, whereas the third also includes thalamic populations (Fig. 1). These models are all biophysically informed, meaning that they offer direct insight into the biophysical parameters underlying the simulated oscillations. Mathematically speaking, all models give rise to diverse behavior (Coombes 2010), ranging from noise-induced fluctuations to self-sustained oscillations (limit cycles, see Box 1). Moreover, all models facilitate phenomena such as multistability where different states (e.g., noise-induced fluctuations and self-sustained oscillations of different frequencies and amplitudes) can coexist. Such diversity is key to mimicking the complex behavior measured empirically. The rationale for using 3 different models was simply to increase robustness by showing that our primary findings generalize across models. For example, if an optimal working point that matches MEG findings can be found near a certain bifurcation in one model, then this should generalize to other models that support the same bifurcation. Furthermore, the Liley model has a richer bifurcation diagram, and this allows us to analyze the potential role for various other bifurcations.

Neuronal models: neural units for each model are shown. The Wilson–Cowan and Liley model consist of excitatory populations (Ej) and inhibitory populations (Ij) connected by caa and Waa’, where a denotes the type of connection (excitatory (E) or inhibitory (I)), and the use of different names w and c is to distinguish the parameters between the models. The parameters Paa’ correspond to input currents to each population and Pnetwork refers to the input from the wider network to each region. The Robinson model consists of a cortical population and a thalamic population (reticular nucleus and relay nucleus), which are connected to each other by coupling parameters aa (e, i, r, s). The variables Qa and φa correspond to mean firing rates and damped firing rates, respectively.

Figure 1.

Neuronal models: neural units for each model are shown. The Wilson–Cowan and Liley model consist of excitatory populations (_E_j) and inhibitory populations (_I_j) connected by _c_aa and _W_aa’, where a denotes the type of connection (excitatory (E) or inhibitory (I)), and the use of different names w and c is to distinguish the parameters between the models. The parameters _P_aa’ correspond to input currents to each population and _P_network refers to the input from the wider network to each region. The Robinson model consists of a cortical population and a thalamic population (reticular nucleus and relay nucleus), which are connected to each other by coupling parameters aa (e, i, r, s). The variables _Q_a and _φ_a correspond to mean firing rates and damped firing rates, respectively.

Box 1.

Concepts from nonlinear dynamics

Bifurcation: a regime in which a small change in the value of a model parameter generates a qualitative change in system behavior.
Attractor: subspace in phase space to which a trajectory converges for certain initial conditions.
Stability: a system or attractor is stable when a trajectory with initial conditions sufficiently close to it evolves toward it in forward time. If all trajectories with initial conditions near to the attractor diverges from the attractor, then the attractor is said to be unstable.
Fixed point: a trajectory of a dynamical system that does not change in time. It is also often called an equilibrium or steady-state.
Noise-induced fluctuations: if the system resides in a stable fixed point regime, a small perturbation will exponentially decay to the steady-state, and thus weak noisy input can induce fluctuations around the steady-state.
Limit cycle: an isolated periodic orbit (oscillator), which is a solution for a dynamical system which repeats itself in time. A limit cycle can either be stable or unstable.
Saddle-node bifurcation: by tuning a bifurcation parameter of interest, an unstable fixed point approaches a stable fixed point until they collide and annihilate.
Hopf bifurcation: a fixed point undergoes a linear instability when a pair of complex conjugate eigenvalues of the linearized system change from negative to positive real parts, leading to the creation of a small amplitude periodic oscillation. This is said to be supercritical if the emerging orbit is stable, and subcritical otherwise.
Periodic doubling bifurcation: an oscillatory system switches to a new behavior with twice the period of the original system.
Homoclinic bifurcation: a global bifurcation typified by the collision of a limit cycle and a fixed point.
Multistability: coexistence of multiple attractors. For the models used here, a Hopf bifurcation with 2 nearby located saddle-node bifurcations gives rise to multistability of a limit cycle and a fixed point.
Bifurcation: a regime in which a small change in the value of a model parameter generates a qualitative change in system behavior.
Attractor: subspace in phase space to which a trajectory converges for certain initial conditions.
Stability: a system or attractor is stable when a trajectory with initial conditions sufficiently close to it evolves toward it in forward time. If all trajectories with initial conditions near to the attractor diverges from the attractor, then the attractor is said to be unstable.
Fixed point: a trajectory of a dynamical system that does not change in time. It is also often called an equilibrium or steady-state.
Noise-induced fluctuations: if the system resides in a stable fixed point regime, a small perturbation will exponentially decay to the steady-state, and thus weak noisy input can induce fluctuations around the steady-state.
Limit cycle: an isolated periodic orbit (oscillator), which is a solution for a dynamical system which repeats itself in time. A limit cycle can either be stable or unstable.
Saddle-node bifurcation: by tuning a bifurcation parameter of interest, an unstable fixed point approaches a stable fixed point until they collide and annihilate.
Hopf bifurcation: a fixed point undergoes a linear instability when a pair of complex conjugate eigenvalues of the linearized system change from negative to positive real parts, leading to the creation of a small amplitude periodic oscillation. This is said to be supercritical if the emerging orbit is stable, and subcritical otherwise.
Periodic doubling bifurcation: an oscillatory system switches to a new behavior with twice the period of the original system.
Homoclinic bifurcation: a global bifurcation typified by the collision of a limit cycle and a fixed point.
Multistability: coexistence of multiple attractors. For the models used here, a Hopf bifurcation with 2 nearby located saddle-node bifurcations gives rise to multistability of a limit cycle and a fixed point.
Bifurcation: a regime in which a small change in the value of a model parameter generates a qualitative change in system behavior.
Attractor: subspace in phase space to which a trajectory converges for certain initial conditions.
Stability: a system or attractor is stable when a trajectory with initial conditions sufficiently close to it evolves toward it in forward time. If all trajectories with initial conditions near to the attractor diverges from the attractor, then the attractor is said to be unstable.
Fixed point: a trajectory of a dynamical system that does not change in time. It is also often called an equilibrium or steady-state.
Noise-induced fluctuations: if the system resides in a stable fixed point regime, a small perturbation will exponentially decay to the steady-state, and thus weak noisy input can induce fluctuations around the steady-state.
Limit cycle: an isolated periodic orbit (oscillator), which is a solution for a dynamical system which repeats itself in time. A limit cycle can either be stable or unstable.
Saddle-node bifurcation: by tuning a bifurcation parameter of interest, an unstable fixed point approaches a stable fixed point until they collide and annihilate.
Hopf bifurcation: a fixed point undergoes a linear instability when a pair of complex conjugate eigenvalues of the linearized system change from negative to positive real parts, leading to the creation of a small amplitude periodic oscillation. This is said to be supercritical if the emerging orbit is stable, and subcritical otherwise.
Periodic doubling bifurcation: an oscillatory system switches to a new behavior with twice the period of the original system.
Homoclinic bifurcation: a global bifurcation typified by the collision of a limit cycle and a fixed point.
Multistability: coexistence of multiple attractors. For the models used here, a Hopf bifurcation with 2 nearby located saddle-node bifurcations gives rise to multistability of a limit cycle and a fixed point.
Bifurcation: a regime in which a small change in the value of a model parameter generates a qualitative change in system behavior.
Attractor: subspace in phase space to which a trajectory converges for certain initial conditions.
Stability: a system or attractor is stable when a trajectory with initial conditions sufficiently close to it evolves toward it in forward time. If all trajectories with initial conditions near to the attractor diverges from the attractor, then the attractor is said to be unstable.
Fixed point: a trajectory of a dynamical system that does not change in time. It is also often called an equilibrium or steady-state.
Noise-induced fluctuations: if the system resides in a stable fixed point regime, a small perturbation will exponentially decay to the steady-state, and thus weak noisy input can induce fluctuations around the steady-state.
Limit cycle: an isolated periodic orbit (oscillator), which is a solution for a dynamical system which repeats itself in time. A limit cycle can either be stable or unstable.
Saddle-node bifurcation: by tuning a bifurcation parameter of interest, an unstable fixed point approaches a stable fixed point until they collide and annihilate.
Hopf bifurcation: a fixed point undergoes a linear instability when a pair of complex conjugate eigenvalues of the linearized system change from negative to positive real parts, leading to the creation of a small amplitude periodic oscillation. This is said to be supercritical if the emerging orbit is stable, and subcritical otherwise.
Periodic doubling bifurcation: an oscillatory system switches to a new behavior with twice the period of the original system.
Homoclinic bifurcation: a global bifurcation typified by the collision of a limit cycle and a fixed point.
Multistability: coexistence of multiple attractors. For the models used here, a Hopf bifurcation with 2 nearby located saddle-node bifurcations gives rise to multistability of a limit cycle and a fixed point.

For all 3 models, we assumed 78 cortical regions, with each region modeled as a single unit (nodal oscillator) depicted in Figure 1. Regional parcellation was equivalent to that used in the MEG data (i.e., the AAL parcellation). Connection strengths between each region were based upon an available structural connectome derived from diffusion tensor weighted imaging (Gong et al. 2009). Mathematical details of the neuronal models and explanation of the goodness-of-fit measure can be found in the supplementary information.

Two type of simulations are performed: (1) tuning of global coupling strength and (2) tuning bifurcation parameters of interest for the different models (see Results section). Simulations are performed for only the alpha band and overall functional connectivity between our 78 cortical regions is tuned by scaling the global structural coupling parameter, η. In doing this, the relative strength of coupling between region pairs is maintained (i.e., is set by the structural connectome) but the overall coupling strength can be increased or decreased. For the purposes of these simulations, structural coupling strength is maintained within a range known to promote a “weakly coupled” regime (i.e., a perturbation that causes the trajectory to jump not far from the intrinsic limit cycle). This regime has been shown in previous work to adequately mimic MEG data (Tewarie, Hillebrand et al. 2016). In order to estimate how well the network model obtained, Cmodel(PDD(f),E(f)) values described the empirical MEG obtained C(PDD(f),E(f)) values (see Methods), we computed the Kullback–Leibler divergence _D_KL between these 2 distributions. _D_KL measures the dissimilarity between 2 distributions. To obtain a measure of similarity, we computed 1/_D_KL as our goodness-of-fit (GOF) estimate.

Results

Amplitude Envelopes and Dynamic Connectivity in a Sensorimotor Task

To assess whether task-induced changes in regional oscillatory amplitude are related to long-range (i.e., between AAL region) phase connectivity, we used MEG data recorded during our self-paced motor paradigm and a visuomotor task. All results pertaining to the latter are shown in Supplementary material (Supplementary Figs S3 and S4). Figure 2B (blue trace) shows the beta band amplitude envelope averaged over trials and subjects, for the right primary sensory region; note the characteristic (Jurkiewicz et al. 2006) modulation with a loss of amplitude (beta desynchronization) during movement and the increase above baseline (beta rebound) following movement cessation. In Figure 2A, mean regional amplitude is shown for all AAL regions, averaged within our 4 time windows prior to the button press (−10 s < t < −5 s); during beta desynchronization (−3 s < t < 0.5 s); during the beta rebound (1.5 < t < 3.5 s), and following beta rebound (5 < t < 10 s). Consistent with previous literature, the beta modulation was dominant in contralateral (right) sensorimotor cortex, however, note also an ipsilateral response suggesting that the poststimulus beta increase exists throughout the sensorimotor network.

Task-based change in local oscillatory dynamics and between region connectivity. (A) Regional beta band oscillatory amplitude for 78 cortical areas in 4 time windows: prior to the button press (−10 to 5 s), during beta suppression (−3 to 0.5 s), during beta rebound (1.5–3.5 s) and on return to baseline (5–10 s). Task-induced effects are strongest in contralateral sensorimotor areas, both during the beta suppression and rebound. (B) Timecourse of beta amplitude envelope (blue) and beta amplitude for surrogate data (red) in the right sensory cortex. (C) Timecourse of beta band connection strength between sensory cortex and the rest of the brain (blue). The red traces show the same calculation for surrogate data. For both B and C, color-coded blocks correspond to the windows used in Panels A, D, and E. (D) Bar chart showing peak-to-peak modulation of beta connectivity (between the right sensory cortex and the rest of the brain) for real and surrogate data. (E) Average connectivity as a function of AAL region for the 4 time windows. Note that during the beta rebound, connectivity is dominated by right sensorimotor cortex. (F) The dominant (strongest 5%) connections during the 4 windows. Note that the dominant network during the beta rebound is characterized by sensorimotor connections.

Figure 2.

Task-based change in local oscillatory dynamics and between region connectivity. (A) Regional beta band oscillatory amplitude for 78 cortical areas in 4 time windows: prior to the button press (−10 to 5 s), during beta suppression (−3 to 0.5 s), during beta rebound (1.5–3.5 s) and on return to baseline (5–10 s). Task-induced effects are strongest in contralateral sensorimotor areas, both during the beta suppression and rebound. (B) Timecourse of beta amplitude envelope (blue) and beta amplitude for surrogate data (red) in the right sensory cortex. (C) Timecourse of beta band connection strength between sensory cortex and the rest of the brain (blue). The red traces show the same calculation for surrogate data. For both B and C, color-coded blocks correspond to the windows used in Panels A, D, and E. (D) Bar chart showing peak-to-peak modulation of beta connectivity (between the right sensory cortex and the rest of the brain) for real and surrogate data. (E) Average connectivity as a function of AAL region for the 4 time windows. Note that during the beta rebound, connectivity is dominated by right sensorimotor cortex. (F) The dominant (strongest 5%) connections during the 4 windows. Note that the dominant network during the beta rebound is characterized by sensorimotor connections.

In Figure 2C, the blue trace shows the timecourse of average connectivity (assessed via PDD) for the right primary sensory AAL region. This timecourse represents connectivity between the region of interest to all other AAL regions, averaged across trials and subjects. Notice that the timecourse of average connectivity mirrors the amplitude envelope (Fig. 2B); at the time of the button press a drop in connectivity is observed, consistent with the notion that this region becomes less connected to a wider network. Following movement, connectivity is elevated suggesting a reintegration with a wider network. Most importantly, the observable relationship between amplitude envelope (Fig. 2_B_—blue trace) and average connectivity (Fig. 2_C_—blue trace) suggests a link between phase-based measures of functional connectivity and oscillatory amplitude.

The averaged amplitude envelope from the surrogate data (generated via randomization of trial start times) is shown in Figure 2B (red trace). Connectivity based upon these surrogate data is shown in Figure 2C (red trace). The bar plot in Figure 2D shows the peak-to-peak change in connectivity (i.e., averaged connectivity change between the beta desynchronization (−3 s < t < 0.5 s) and rebound (1.5 < t < 3.5 s) windows) for the real and surrogate data. The fluctuation in connectivity was significantly larger for real, compared with surrogate data (P = 0.008). Note that an equivalent Figure, with surrogate data based upon a uniform phase randomization process, is shown in Supplementary Fig. S5. Using this alternative methodology, the effect remained significant with P = 0.01.

Figure 2D shows average connectivity plotted as a function of cortical area, for the same time windows as those used in Figure 2A. Prior to the button press (−10 s to −5 s) connectivity is dominated by occipital and parietal connections with some frontal projections. Connectivity decreases during the beta desynchronization (−3 s to 0.5 s; P = 0.035; see also Supplementary Fig. S6) and is subsequently elevated during the beta rebound (P < 0.001) and focused in the contralateral sensorimotor areas. This focal pattern undergoes a change back to a more occipital and parietal distributed pattern of connectivity following rebound cessation. Finally, Figure 2E shows the actual connections that dominate during each window; note the changing spatial signature with occipito-parietal and fronto-parietal connections dominating at rest, and sensorimotor connections dominating during the beta rebound. We repeated the same analysis for the alpha band, however, we observed only a movement induced alpha suppression in the contralateral sensorimotor cortex, without rebound and dynamic connectivity timecourses showed no significant modulation (Supplementary Fig. S7).

Amplitude Envelopes and Dynamic Connectivity in the Resting State

The above analyses show that transient changes in long-range connectivity “between regions” relate to the “within region” amplitude envelope of beta oscillations, suggesting that the well-known beta rebound is generated by a network in which the primary sensorimotor region, which has operated independently during movement (and whose connectivity transiently dropped) is reintegrated into a wider network. We now seek to generalize this relationship by testing whether similar relationships can be found in resting-state data; in this way, we seek to establish a similar effect in a task-free regime, across multiple frequency bands, and across all cortical regions. Figure 3A shows a single example in which instantaneous envelope amplitude in the left inferior parietal cortex is shown on the _x_-axis, plotted against instantaneous connectivity between the same region and the rest of the brain, on the _y_-axis. Note the positive correlation. Results in Figure 3 are based on a threshold k of 0.3. This same analysis was undertaken for all brain regions and across frequency bands ranging from delta to high gamma. To assess statistical significance, surrogate data were used to construct an empirical null-distribution, and we tested whether the genuine correlations, for every region, and frequency band, were significantly higher than the null-distribution.

Resting-state data. (A) Example scatter plot extracted from the left inferior parietal cortex and showing instantaneous oscillatory envelope (x-axis) plotted against average connectivity between the same regions and all other areas. (B) Bar chart shows strength of correlation between envelope amplitude and average connectivity, collapsed over regions and subjects and plotted as a function of frequency. The inset images show the slopes of the connectivity/amplitude relationship for all brain regions in which there was a significant correlation compared with the null-distribution (first method). Brain regions with no significant relationship are shown in gray. Slopes are normalized by the maximum. Note that a spatially widespread significant relationship is observed between envelope amplitude and connectivity in frequency bands ranging from delta to low gamma (* represents P < 0.001). (C) Shows functional connectivity matrices based on windows with high amplitude oscillations and low amplitude oscillations (D). Note increased connectivity and structure when envelopes have high amplitude.

Figure 3.

Resting-state data. (A) Example scatter plot extracted from the left inferior parietal cortex and showing instantaneous oscillatory envelope (_x_-axis) plotted against average connectivity between the same regions and all other areas. (B) Bar chart shows strength of correlation between envelope amplitude and average connectivity, collapsed over regions and subjects and plotted as a function of frequency. The inset images show the slopes of the connectivity/amplitude relationship for all brain regions in which there was a significant correlation compared with the null-distribution (first method). Brain regions with no significant relationship are shown in gray. Slopes are normalized by the maximum. Note that a spatially widespread significant relationship is observed between envelope amplitude and connectivity in frequency bands ranging from delta to low gamma (* represents P < 0.001). (C) Shows functional connectivity matrices based on windows with high amplitude oscillations and low amplitude oscillations (D). Note increased connectivity and structure when envelopes have high amplitude.

In Figure 3B, the bar chart shows the strength of the correlation between envelope amplitude and connectivity, averaged over regions within individual frequency bands. Note that the relationship peaks in the theta to low gamma range although significance compared with surrogate data (denoted by ∗) extends to the delta band. The inset brain images show the cortical AAL regions that are driving this relationship for each frequency band. The slopes are plotted for regions that show a significant correlation. Note that correlations are found to be widespread across the entire cortex, but where slopes are strongest in the alpha band in occipital regions and in beta in the sensorimotor regions. As a post hoc analysis, we selected time windows of low and high amplitude oscillations (threshold based on the mean amplitude envelopes (over regions) plus/minus 2 standard deviations, respectively), and computed phase connectivity between all possible AAL region pairs, represented as a 78 × 78 weighted adjacency matrix, within these windows. Results are shown in Figure 3C and D. Note that across all frequency bands, high oscillatory amplitude coexists with strong network connectivity structure whereas low amplitude oscillations coincide with weak (unstructured) network connectivity. Note that results obtained with surrogate data based upon a uniform phase randomization process show a similar effect (Supplementary Fig. S8).

Resting-State MEG Data and Large-Scale Neuronal Models

In the previous sections, we demonstrated a relationship between local oscillatory amplitude and dynamic connectivity. Specifically, periods of high oscillatory amplitude coincide (beyond chance) with periods of strong connectivity. This finding potentially provides a novel interpretation for measured neural oscillations—with high local amplitude linked to long-range (between region) connectivity. However, it also raises further questions: Is this finding an inherent property of oscillations themselves (i.e., does high connectivity “drive” high amplitude (or vice versa))? Or are such empirical relationships merely driven by SNR? Further, assuming the relationship is genuine, does connectivity to a region drive high amplitude oscillations directly, or do high amplitude oscillations offer natural windows for strong network connectivity. Such questions are challenging to answer using MEG data alone. For this reason, to further support our MEG findings and to take them further, we turned to computational modeling.

Figure 4A shows how the amplitude envelope averaged across time and regions (denoted <_E_>), changes for increasing η. For all models, there is an increase in amplitude when coupling increases, but this increase is less pronounced for Liley’s model, which may be due to model-specific parameters, or the different intrinsic properties of the model. Figure 4B shows functional connectivity, measured using the phase-difference-derivative approach (denoted ), plotted as a function global structural coupling. For calculation of the PDD, we used exactly the same method as was used for MEG data (see eq. (1) in Methods section). Again, we observe a general increase for all 3 models (although the Wilson–Cowan and Liley models plateau for increasing coupling). The correlation between local amplitude envelope and functional connectivity is shown in Figure 4C (denoted <_C_(_E_,PDD)>); this represents our model analog of the empirical relationship predicted by our MEG data (though here we test this relationship without an SNR confound). This relationship exists with all 3 computational models, but scales with coupling strength, η. For the Robinson and Liley model, a general increase in correlation with η is observed although an inverted U-curve is seen for the Wilson–Cowan model (Fig. 4C). Finally, Figure 4D shows an estimate of goodness-of-fit (in terms of the inverse of the Kullback–Leibler divergence) between MEG and simulated data as a function of η, and suggests that the best fit with MEG data for all 3 models occurs for intermediate values of coupling. Overall, the evidence in Figure 4 suggests that a significant relationship between oscillatory amplitude and functional connectivity is an inherent property of modeled neural oscillators. However, the existence of such a relationship depends critically on overall coupling strength in the network.

The effect of global coupling on the relationship between local amplitude envelopes and dynamic connectivity in 3 computational models. (A) Mean envelope amplitude (<E>, collapsed over all space and time) versus global coupling strength, η. (B) Mean functional connectivity (<PDD>, again collapsed over all space and time) versus η. (C) Correlation between timecourses of envelope amplitude and dynamic connectivity versus η—note the empirical relationship shown in Figures 3 and 4 is an inherent property of all 3 models. (D) Shows the goodness-of-fit (GOF) between model obtained correlations and empirical resting-state data for all 3 models.

Figure 4.

The effect of global coupling on the relationship between local amplitude envelopes and dynamic connectivity in 3 computational models. (A) Mean envelope amplitude (<_E_>, collapsed over all space and time) versus global coupling strength, η. (B) Mean functional connectivity (, again collapsed over all space and time) versus η. (C) Correlation between timecourses of envelope amplitude and dynamic connectivity versus _η_—note the empirical relationship shown in Figures 3 and 4 is an inherent property of all 3 models. (D) Shows the goodness-of-fit (GOF) between model obtained correlations and empirical resting-state data for all 3 models.

We further used the models to analyze the influence of local dynamics on the correlation between amplitude envelope and dynamic connectivity. Here, we reasoned that if driving up the local amplitudes leads to a stronger correlation, this would suggest that high amplitude oscillations form natural windows for increased long-range coupling. A critical concept here is a “bifurcation”, which corresponds to a regime in which a small change in the value of a specific modeled parameter generates a substantial change in the qualitative behavior of a system. Example, bifurcation diagrams are shown in Figure 5A_–_C; note diagrams are presented for uncoupled regions only and the concepts used in this section are explained in Box 1. Taking the Robinson model as an example, the bifurcation diagram shows firing rate output of the excitatory population, _ϕ_e, plotted against a thalamocortical coupling parameter _ν_es. The diagram clearly shows 3 regimes: regime 1 is characterized by a stable fixed point. Here, amplitude fluctuations around 10 Hz are induced by input noise (see Box 1); regime 2 is characterized by a bistable state with a stable fixed point, stable and unstable oscillations, while regime 3 is characterized by stable oscillations and an unstable fixed point. Given previous studies on the importance of multistability that can occur around a Hopf bifurcation for local dynamics (Freyer et al. 2011), we have chosen parameters for the bifurcation diagrams such that we could observe the appearance and disappearance of a Hopf bifurcation and at the same time a jump in amplitude (see Methods). In the case for the Wilson–Cowan model, the parameter of interest to elicit a Hopf bifurcation is an external input current, and for the Liley model the r and k parameters correspond to the level of inhibition to both excitatory and inhibitory or to the inhibitory population only. Therefore, note that qualitatively different parameters of the models can give rise to a Hopf bifurcation.

The effect of local dynamics. The upper panels show the bifurcation diagrams for single units for all 3 neuronal models. For the Robinson model (A), we vary the thalamocortical coupling parameter (ves) and we observe fixed points (stable: continuous black, unstable: dashed black), stable (continuous red) and unstable oscillations (dashed red), with saddle and Hopf bifurcations. Regime 1 thus corresponds to a stable fixed point, regime 2 a bistable state with a stable fixed point, stable, and unstable oscillations, while regime 3 to an unstable fixed point and stable oscillations. For the Wilson–Cowan model (B), we can observe stable fixed points (red lines), unstable fixed points (black lines) with 2 saddle-node bifurcations and a Hopf bifurcation (with the amplitude of unstable emergent oscillations shown with blue circles) when the external input (Pext) to the single unit is varied. The two-parameter Liley bifurcation diagram (C) shows saddle (red), Hopf (blue), periodic doubling (green), homoclinic (yellow) for bifurcation parameters k (drive to inhibitory population), and r (drive to excitatory population). Panels D–F show (i) mean amplitude envelopes, (ii) mean dynamic connectivity, (iii) the correlation between the amplitude envelopes and dynamic connectivity, and (iv) the goodness-of-fit of this correlation to empirical data. (D) shows the case for the Robinson model, (E) shows the case for the Wilson–Cowan model and F shows the case for the Liley model.

Figure 5.

The effect of local dynamics**.** The upper panels show the bifurcation diagrams for single units for all 3 neuronal models. For the Robinson model (A), we vary the thalamocortical coupling parameter (_v_es) and we observe fixed points (stable: continuous black, unstable: dashed black), stable (continuous red) and unstable oscillations (dashed red), with saddle and Hopf bifurcations. Regime 1 thus corresponds to a stable fixed point, regime 2 a bistable state with a stable fixed point, stable, and unstable oscillations, while regime 3 to an unstable fixed point and stable oscillations. For the Wilson–Cowan model (B), we can observe stable fixed points (red lines), unstable fixed points (black lines) with 2 saddle-node bifurcations and a Hopf bifurcation (with the amplitude of unstable emergent oscillations shown with blue circles) when the external input (_P_ext) to the single unit is varied. The two-parameter Liley bifurcation diagram (C) shows saddle (red), Hopf (blue), periodic doubling (green), homoclinic (yellow) for bifurcation parameters k (drive to inhibitory population), and r (drive to excitatory population). Panels D_–_F show (i) mean amplitude envelopes, (ii) mean dynamic connectivity, (iii) the correlation between the amplitude envelopes and dynamic connectivity, and (iv) the goodness-of-fit of this correlation to empirical data. (D) shows the case for the Robinson model, (E) shows the case for the Wilson–Cowan model and F shows the case for the Liley model.

Figures 5D_–_F show (i) amplitude envelope (), (ii) dynamic connectivity (), (iii) correlation between amplitude envelope and dynamic connectivity (C(E,PDD), and (iv) goodness-of-fit (GOF)) to MEG data as a function of bifurcation parameters. For the Robinson model (Fig. 5A,D), local amplitude increases when the bifurcation parameter is increased in an almost sigmoid form and note that the Hopf bifurcation (_ν_es= 1.6) occurs before the peak amplitude. Dynamic connectivity on the other hand peaks around the Hopf bifurcation; this value also shows the strongest correlation between dynamic connectivity and amplitude envelopes. It also shows the highest goodness-of-fit with MEG data. This model illustrates nicely that the observed relationship between dynamic connectivity and amplitude envelope modulations is unlikely to be an SNR effect. The noise in the model is fixed, whereas the genuine amplitude of the signal φ e increases with an increase of the bifurcation parameter _ν_es. The maximum amplitude in the model occurs well after the Hopf bifurcation (_ν_es= 1.7; see Fig. 5D(i)), which is not located at the best fit with the empirical data. For the Wilson–Cowan (Fig. 5B,E), local amplitude increases when the bifurcation parameter is increased, and a similar pattern is seen for dynamic connectivity. However, again the maximum correlation between amplitude envelope and dynamic connectivity is seen around the Hopf bifurcation (_P_ext = 7), for which we also find the best fit with the MEG data. This means that again the maximum correlation between amplitude envelopes and dynamic connectivity does not coincide with the regime with the largest SNR, that is, again suggesting that this is not an SNR effect. The Liley model shows a rich bifurcation diagram with several other bifurcations apart from the Hopf (Fig. 5C,F), such as periodic doubling bifurcation, or a homoclinic bifurcation. This allows us to test whether the match the good match with MEG data is specific to the Hopf bifurcation or whether this can be extended to other bifurcations as well. However, a working point near the Hopf bifurcation outperforms the other types of bifurcations to explain our MEG findings. Results show that the relationship between amplitude envelopes and dynamic connectivity is not merely driven by tuning up the local amplitudes, but is optimized in a regime where the oscillators themselves exist in a multistable state.

Discussion

In this paper, we hypothesized that modulations in neural oscillatory amplitude, at some brain location of interest, reflect (at least in part) modulations in connectivity between that region and the rest of the brain. During task-induced modulation of oscillatory amplitude, and during the natural fluctuations in resting-state amplitude, empirical data suggest that this is the case. In a self-paced motor task, we showed that beta band connectivity between primary sensorimotor cortex and the rest of the brain mirrored the well-established beta band movement-related amplitude decrease and post-movement beta rebound. In the resting state, we showed significant correlations between amplitude envelopes and dynamic functional connectivity, across multiple frequency bands ranging from theta to low gamma. Our MEG findings mirror results from 3 separate computational neural models which show that correlation between oscillatory envelopes and functional connectivity is an inherent property of a coupled nodal oscillator model. However, our models also show that this coupling is optimized when nodal oscillators exist in a multistable regime. This latter point potentially offers new insights into the way in which measured oscillations might be analyzed or characterized.

Our task data show how functional connectivity, between contralateral sensory cortex and the rest of the brain, varies during a simple self-paced button press task. Temporally resolved connectivity (Fig. 2C) mirrored the previously well-characterized timecourse of beta band oscillatory envelope and showed a significant reduction during movement, followed by an increase on movement cessation. This is consistent with the idea that regions of a functional network operate in isolation during information encoding, and the poststimulus response represents a mechanism by which that isolated region reintegrates into the wider sensorimotor network. Indeed, it’s tempting to speculate that during the poststimulus rebound, a top-down inhibitory influence is placed on the primary sensorimotor regions by a broader brain network. Moreover, that influence is mediated by rhythmic fluctuations which are coherent between the wider network and the primary cortex. As coherent clusters of action potentials from other (e.g., motor planning) regions arrive, this would not only increase phase locking between the (postsynaptic driven) MEG signals originating from contralateral motor cortex and the wider network, but it would also likely drive the local population into synchrony, thus increasing the amplitude of local beta oscillations. This concept, which is nothing more than the communication by coherence hypothesis, is supported qualitatively by our connectivity graphs (Fig. 2E) which show that the dominant connections during the rebound phase are within the distributed motor network including bilateral sensory, motor and premotor cortices as well as the supplementary motor area.

It should be pointed out that this association between beta oscillations and interregional connectivity is not new. In general, beta band oscillations have been shown to be coordinated across spatially separate brain regions and spatial analyses of those connections suggest that they may also be responsible for the mediation of functional connectivity in several well-characterized resting-state networks (Brookes et al. 2011; Hipp et al. 2012). Further, beta oscillations have been shown to be linked to structural connections in the brain (Hunt et al. 2016). In sensorimotor tasks, increases in beta oscillations following task cessation have been observed robustly for many years (Cassim et al. 2001; Jurkiewicz et al. 2006). Work on the mu rhythm also suggests that amplitude modulations in the motor cortex are initiated by observing hand movements executed by a third person (Rizzolatti et al. 2001). This gives rise to the notion that distant input from other areas, for example, the visual areas, modulates oscillations in the motor cortex. However, although the link to connectivity has been suggested, direct evidence has been lacking. Here, the good time resolution of the PDD connectivity metric has allowed a direct comparison and this dynamic relationship can be elucidated.

The observable relationship between connectivity and amplitude in the resting-state extends our task-based finding, generalizing it to multiple brain regions and frequency bands. Our primary resting-state finding (Fig. 3) was that functional connectivity and amplitude envelope timecourses correlate beyond chance in frequency bands stretching from the delta to the low gamma. Spatially, this correlation was widespread, suggesting that this is a general phenomenon. However, there was some inhomogeneity across brain regions, with, for example, the slope of the correlation maximum in occipital regions for the alpha band and in sensorimotor areas for the beta band. The fact that no significant relationship was observed in the higher frequency bands likely reflects the fact that the SNR in those frequency bands is poor. For many years, fluctuations in the amplitude envelope of oscillations in a given region have been considered to be a marker of “activity” or dynamics local to that region. However, the evidence presented here shows that interpreting such fluctuations as a purely local phenomenon is inaccurate. Rather, we might speculate that high amplitude oscillations are generated by synchronization of dendritic currents across the local neural assembly, and such synchronization is brought about by an increase in interregional connectivity. Again, the relationship between oscillations and connectivity is not new, although it perhaps contrasts with views that, for example, alpha oscillations are “idling rhythms”. Indeed, our findings are in strong agreement with others (Hummel and Gerloff 2006; Dugué et al. 2011; Zumer et al. 2014) in showing that alpha-band oscillations are crucial in communications across brain areas. It is important to note however that the relationships shown here are neither causal, nor one-to-one. Whilst in general, we have shown that periods of high amplitude tend to coincide with periods of high connectivity, this may not be the case all of the time. We, therefore, advocate that future studies of neural oscillations also take into account long-range (between region) connectivity as well as local oscillatory amplitude (i.e., a time-frequency–spectral approach).

Although our results point towards a direct relationship between oscillatory amplitude and connectivity, extreme caution is necessary when interpreting connectivity information from data with diverging amplitudes, since phase locking is amplitude (or more specifically signal-to-noise ratio (SNR)) dependent (Muthukumaraswamy and Singh 2013). Specifically, it is likely that an SNR increase might drive an increase in apparent phase locking simply because we are better able to estimate instantaneous phase. Using only empirical data, this is a difficult problem to address since, whilst SNR increases likely result in better ability to measure connectivity, an alternative hypothesis is that if neural oscillations mediate connectivity, then the amplitude (hence SNR) itself is driven directly by connectivity. In other words, it would be impossible to observe a connectivity change in real brain data without a simultaneous change in amplitude, hence SNR. This therefore remains a conceptual challenge.

Here, we attempted to deal with the SNR confound by drawing an analogy with simulations, in which neural interactions can be simulated without confounds of instrumental SNR. These simulations supported our MEG findings in showing that correlation between functional connectivity and envelope amplitude is a natural property that emerges from realistic biophysically informed models of large-scale neural networks. Critically, our simulations suggest that the connectivity/amplitude relationship is dependent on the parameters chosen: Firstly, sufficient global coupling is required before empirically observed correlations between amplitude envelopes and dynamic connectivity are obtained, suggesting that modulations in connectivity to a region drive, to some extent, modulations in amplitude envelopes. Secondly, our empirically observed connectivity/amplitude correlations seem to be most consistent with an operating point for the local network regions near a Hopf bifurcation. In our case, this Hopf bifurcation is located near 2 saddle-node bifurcations, together giving rise to a regime with multistability. Such multistability allows for switching between low and high amplitude modes in the presence of noise and these high amplitude modes co-occur with periods of strong connectivity. A previous study that only analyzed the multistability of alpha amplitude modulations also emphasized the important role of residing in the vicinity of a Hopf bifurcation (Freyer et al. 2009, 2011). The current findings extend this by showing the potential role for Hopf bifurcations not merely for amplitude modulations, but for amplitude modulations in relation to dynamic connectivity. Previous modeling studies also suggest that the influence of amplitude on connectivity is a genuine effect (Daffertshofer and van Wijk 2011; Moon et al. 2015). Previous work has further shown that phase dynamics of one region in a Kuramoto-like system, derived from Wilson–Cowan oscillators, is intrinsically dependent on the amplitude dynamics of all regions that the region is connected to in a nontrivial way (Daffertshofer and van Wijk 2011). In summary, these results suggest that the connectivity/amplitude relationship is a property that we should expect to see in empirical data, and as such is unlikely to be (completely) driven by SNR confounds. However, there is no trivial relationship between the amplitude of local dynamics and the likelihood of connectivity, that is, turning up local amplitudes is not the sole requirement, but each nodal oscillator must itself exist in a multistable regime, close to a Hopf bifurcation.

Our findings from simulations suggest that the empirical amplitude/connectivity relationship that we see in real MEG data is not a sole property of SNR. However, even given this good agreement, no causal relationship has been demonstrated and statements about the role of oscillations can only be made by analogy. This is a limitation of this study and future work will be required to further this empirical demonstration. As mentioned above, in MEG data alone it is extremely challenging to separate SNR from connectivity. However, in future studies, neuromodulation might offer a means to examine how changing the amplitude of oscillations in some specific brain region changes connectivity within a broader brain network to which the modulated region belongs. Specifically, there are a number of recent studies showing that the use of either transcranial magnetic stimulation (TMS) or transcranial alternating current stimulation (tACS) can enhance the amplitude of cortical rhythms—particularly the alpha rhythm (Thut et al. 2011; Vossen et al. 2015). Further, this enhancement has been shown to persist even after cessation of stimulation. It is therefore possible to conceive a study where one measures resting-state oscillatory amplitude and connectivity, within some region of interest, before and after TMS or tACS and then relate the 2 metrics. This, in the future, may offer a means to generate a causal link between amplitude and connectivity.

Conclusion

We have demonstrated that amplitude envelope modulations during both task and resting-state mirror changes in functional connectivity. This relationship between local amplitude modulation and dynamic connectivity to or from this region may be the result of a bidirectional causality. This requires sufficient global network coupling but also multistability at the regional level that allows for switching between low and high amplitude modes that may be natural windows for synchronization. Our results facilitate a change in the way in which invasive and noninvasive electrophysiological data are interpreted; specifically, the notion that shifts in oscillatory amplitude represent only local dynamics need to be revised to an interpretation where increased oscillatory amplitude can be driven by increased long-range connectivity to that region. This new interpretation will allow for new insights not only in basic science, but also in disease, in particular for pathologies where impaired connectivity is implicated.

Notes

Conflict of Interest: None declared.

References

Abeysuriya

RG

,

Hadida

J

,

Sotiropoulos

SN

,

Becker

R

,

Hunt

BAE

,

Brookes

MJ

,

Woolrich

MW

.

2018

.

A biophysical model of dynamic balancing of excitation and inhibition in fast oscillatory large-scale networks

.

PLOS Comput. Biol

.

14

:

e1006007

.

Baker

AP

,

Brookes

MJ

,

Rezek

IA

,

Smith

SM

,

Behrens

T

,

Smith

PJP

,

Woolrich

M

.

2014

.

Fast transient networks in spontaneous human brain activity

.

Elife

.

3

:

e01867

.

Benjamini

Y

,

Hochberg

Y

.

1995

.

Controlling the false discovery rate: a practical and powerful approach to multiple testing

.

J R Stat Soc Ser B

.

57

:

289

300

.

Breakspear

M

,

Williams

LM

,

Stam

CJ

.

2004

.

A novel method for the topographic analysis of neural activity reveals formation and dissolution of “dynamic cell assemblies”

.

J Comput Neurosci

.

16

:

49

68

.

Brookes

MJ

,

Hale

JR

,

Zumer

JM

,

Stevenson

CM

,

Francis

ST

,

Barnes

GR

,

Owen

JP

,

Morris

PG

,

Nagarajan

SS

.

2011

.

Measuring functional connectivity using MEG: methodology and comparison with fcMRI

.

Neuroimage

.

56

:

1082

1104

.

Brookes

MJ

,

O’neill

GC

,

Hall

EL

,

Woolrich

MW

,

Baker

A

,

Corner

SP

,

Robson

SE

,

Morris

PG

,

Barnes

GR

.

2014

.

Measuring temporal, spectral and spatial changes in electrophysiological brain network connectivity

.

Neuroimage

.

91

:

282

299

.

Brookes

MJ

,

Tewarie

PK

,

Hunt

BAE

,

Robson

SE

,

Gascoyne

LE

,

Liddle

EB

,

Liddle

PF

,

Morris

PG

.

2016

.

A multi-layer network approach to MEG connectivity analysis

.

Neuroimage

.

132

:

425

438

.

Brookes

MJ

,

Vrba

J

,

Robinson

SE

,

Stevenson

CM

,

Peters

AM

,

Barnes

GR

,

Hillebrand

A

,

Morris

PG

.

2008

.

Optimising experimental design for MEG beamformer imaging

.

Neuroimage

.

39

:

1788

1802

.

Cassim

F

,

Monaca

C

,

Szurhaj

W

,

Bourriez

J-L

,

Defebvre

L

,

Derambure

P

,

Guieu

J-D

.

2001

.

Does post-movement beta synchronization reflect an idling motor cortex?

Neuroreport

.

12

:

3859

3863

.

Colclough

GL

,

Brookes

MJ

,

Smith

SM

,

Woolrich

MW

.

2015

.

A symmetric multivariate leakage correction for MEG connectomes

.

Neuroimage

.

117

:

439

448

.

Coombes

S

.

2010

.

Large-scale neural dynamics: simple and complex

.

Neuroimage

.

52

:

731

739

.

Daffertshofer

A

,

van Wijk

BCM

.

2011

.

On the influence of amplitude on the connectivity between phases

.

Front Neuroinform

.

5

:

6

.

Dugué

L

,

Marque

P

,

VanRullen

R

.

2011

.

The phase of ongoing oscillations mediates the causal relation between brain excitation and visual perception

.

J Neurosci

.

31

:

11889

11893

.

Freyer

F

,

Aquino

K

,

Robinson

PA

,

Ritter

P

,

Breakspear

M

.

2009

.

Bistability and non-Gaussian fluctuations in spontaneous cortical activity

.

J Neurosci

.

29

:

8512

8524

.

Freyer

F

,

Roberts

JA

,

Becker

R

,

Robinson

PA

,

Ritter

P

,

Breakspear

M

.

2011

.

Biophysical mechanisms of multistability in resting-state cortical rhythms

.

J Neurosci

.

31

:

6353

6361

.

Fries

P

.

2005

.

A mechanism for cognitive dynamics: neuronal communication through neuronal coherence

.

Trends Cogn Sci

.

9

:

474

480

.

Fries

P

,

Reynolds

JH

,

Rorie

AE

,

Desimone

R

.

2001

.

Modulation of oscillatory neuronal synchronization by selective visual attention

.

Science

.

291

:

1560

1563

.

Gaetz

W

,

Edgar

JC

,

Wang

DJ

,

Roberts

TPL

.

2011

.

Relating MEG measured motor cortical oscillations to resting γ-aminobutyric acid (GABA) concentration

.

Neuroimage

.

55

:

616

621

.

Gong

G

,

He

Y

,

Concha

L

,

Lebel

C

,

Gross

DW

,

Evans

AC

,

Beaulieu

C

.

2009

.

Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography

.

Cereb Cortex

.

19

:

524

536

.

Hall

SD

,

Stanford

IM

,

Yamawaki

N

,

McAllister

CJ

,

Rönnqvist

KC

,

Woodhall

GL

,

Furlong

PL

.

2011

.

The role of GABAergic modulation in motor function related neuronal network activity

.

Neuroimage

.

56

:

1506

1510

.

Hillebrand

A

,

Barnes

GR

,

Bosboom

JL

,

Berendse

HW

,

Stam

CJ

.

2012

.

Frequency-dependent functional connectivity within resting-state networks: an atlas-based MEG beamformer solution

.

Neuroimage

.

59

:

3909

3921

.

Hipp

JF

,

Hawellek

DJ

,

Corbetta

M

,

Siegel

M

,

Engel

AK

.

2012

.

Large-scale cortical correlation structure of spontaneous oscillatory activity

.

Nat Neurosci

.

15

:

884

890

.

Huang

MX

,

Mosher

JC

,

Leahy

RM

.

1999

.

A sensor-weighted overlapping-sphere head model and exhaustive head model comparison for MEG

.

Phys Med Biol

.

44

:

423

.

Hummel

FC

,

Gerloff

C

.

2006

.

Interregional long-range and short-range synchrony: a basis for complex sensorimotor processing

.

Prog Brain Res

.

159

:

223

236

.

Hunt

BAE

,

Tewarie

PK

,

Mougin

OE

,

Geades

N

,

Jones

DK

,

Singh

KD

,

Morris

PG

,

Gowland

PA

,

Brookes

MJ

.

2016

.

Relationships between cortical myeloarchitecture and electrophysiological networks

.

Proc Natl Acad Sci

.

113

:

13510

13515

.

Jurkiewicz

MT

,

Gaetz

WC

,

Bostan

AC

,

Cheyne

D

.

2006

.

Post-movement beta rebound is generated in motor cortex: evidence from neuromagnetic recordings

.

Neuroimage

.

32

:

1281

1289

.

Lachaux

J-P

,

Rodriguez

E

,

Martinerie

J

,

Varela

FJ

.

1999

.

Measuring phase synchrony in brain signals

.

Hum Brain Mapp

.

8

:

194

208

.

Lackner

CL

,

Marshall

WJ

,

Santesso

DL

,

Dywan

J

,

Wade

T

,

Segalowitz

SJ

.

2014

.

Adolescent anxiety and aggression can be differentially predicted by electrocortical phase reset variables

.

Brain Cogn

.

89

:

90

98

.

Liley

DTJ

,

Cadusch

PJ

,

Dafilis

MP

.

2002

.

A spatially continuous mean field theory of electrocortical activity

.

Netw Comput Neural Syst

.

13

:

67

113

.

Moon

J-Y

,

Lee

U

,

Blain-Moraes

S

,

Mashour

GA

.

2015

.

General relationship of global topology, local dynamics, and directionality in large-scale brain networks

.

PLoS Comput Biol

.

11

:

e1004225

.

Muthukumaraswamy

SD

,

Singh

KD

.

2013

.

Visual gamma oscillations: the effects of stimulus type, visual field coverage and stimulus motion on MEG and EEG recordings

.

Neuroimage

.

69

:

223

230

.

Nichols

TE

,

Holmes

AP

.

2002

.

Nonparametric permutation tests for functional neuroimaging: a primer with examples

.

Hum Brain Mapp

.

15

:

1

25

.

Nolte

G

,

Bai

O

,

Wheaton

L

,

Mari

Z

,

Vorbach

S

,

Hallett

M

.

2004

.

Identifying true brain interaction from EEG data using the imaginary part of coherency

.

Clin Neurophysiol

.

115

:

2292

2307

.

O’Nneill

GC

,

Bauer

M

,

Woolrich

MW

,

Morris

PG

,

Barnes

GR

,

Brookes

MJ

.

2015

.

Dynamic recruitment of resting state sub-networks

.

Neuroimage

.

115

:

85

95

.

O’Neill

GC

,

Tewarie

PK

,

Colclough

GL

,

Gascoyne

LE

,

Hunt

BAE

,

Morris

PG

,

Woolrich

MW

,

Brookes

MJ

.

2017

.

Measurement of dynamic task related functional networks using MEG

.

Neuroimage

.

146

:

667

678

.

Prichard

D

,

Theiler

J

.

1994

.

Generating surrogate data for time series with several simultaneously measured variables

.

Phys Rev Lett

.

73

:

951

.

Rizzolatti

G

,

Fogassi

L

,

Gallese

V

.

2001

.

Neurophysiological mechanisms underlying the understanding and imitation of action

.

Nat Rev Neurosci

.

2

:

661

670

.

Robinson

SE

,

Vrba

J

.

1999

.

Functional neuroimaging by synthetic aperture magnetometry (SAM)

. In:

Yoshimoto

T

,

Kotani

M

,

Kuriki

S

,

Karibe

H

,

Nakasato

N

, editors.

Recent advances in biomagnetism

.

Sendai

:

Tohoku University Press

, p.

302

305

.

Robinson

PA

,

Rennie

CJ

,

Rowe

DL

.

2002

.

Dynamics of large-scale brain activity in normal arousal states and epileptic seizures

.

Phys Rev E

.

65

:

41924

.

Robinson

PA

,

Rennie

CJ

,

Wright

JJ

,

Bahramali

H

,

Gordon

E

,

Rowe

DL

.

2001

.

Prediction of electroencephalographic spectra from neurophysiology

.

Phys Rev E

.

63

:

21903

.

Sarvas

J

.

1987

.

Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem

.

Phys Med Biol

.

32

:

11

.

Stam

CJ

.

2005

.

Nonlinear dynamical analysis of EEG and MEG: review of an emerging field

.

Clin Neurophysiol

.

116

:

2266

2301

.

Tewarie

P

,

Bright

MG

,

Hillebrand

A

,

Robson

SE

,

Gascoyne

LE

,

Morris

PG

,

Meier

J

,

Van Mieghem

P

,

Brookes

MJ

.

2016

a.

Predicting haemodynamic networks using electrophysiology: the role of non-linear and cross-frequency interactions

.

Neuroimage

.

130

:

273

292

.

Tewarie

P

,

Hillebrand

A

,

van Dijk

BW

,

Stam

CJ

,

O’Neill

GC

,

Van Mieghem

P

,

Meier

JM

,

Woolrich

MW

,

Morris

PG

,

Brookes

MJ

.

2016

b.

Integrating cross-frequency and within band functional networks in resting-state MEG: a multi-layer network approach

.

Neuroimage

.

142

:

324

336

.

Thatcher

RW

,

North

DM

,

Biver

CJ

.

2009

.

Self‐organized criticality and the development of EEG phase reset

.

Hum Brain Mapp

.

30

:

553

574

.

Thiele

A

,

Herrero

JL

,

Distler

C

,

Hoffmann

K-P

.

2012

.

Contribution of cholinergic and GABAergic mechanisms to direction tuning, discriminability, response reliability, and neuronal rate correlations in macaque middle temporal area

.

J Neurosci

.

32

:

16602

16615

.

Thut

G

,

Veniero

D

,

Romei

V

,

Miniussi

C

,

Schyns

P

,

Gross

J

.

2011

.

Rhythmic TMS causes local entrainment of natural oscillatory signatures

.

Curr Biol

.

21

:

1176

1185

.

Tzourio-Mazoyer

N

,

Landeau

B

,

Papathanassiou

D

,

Crivello

F

,

Etard

O

,

Delcroix

N

,

Mazoyer

B

,

Joliot

M

.

2002

.

Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain

.

Neuroimage

.

15

:

273

289

.

Uhlhaas

PJ

,

Singer

W

.

2006

.

Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology

.

Neuron

.

52

:

155

168

.

Vidaurre

D

,

Quinn

AJ

,

Baker

AP

,

Dupret

D

,

Tejero-Cantero

A

,

Woolrich

MW

.

2016

.

Spectrally resolved fast transient brain states in electrophysiological data

.

Neuroimage

.

126

:

81

95

.

Vossen

A

,

Gross

J

,

Thut

G

.

2015

.

Alpha power increase after transcranial alternating current stimulation at alpha frequency (α-tACS) reflects plastic changes rather than entrainment

.

Brain Stimul Basic, Transl Clin Res Neuromodulation

.

8

:

499

508

.

Ward

LM

.

2003

.

Synchronous neural oscillations and cognitive processes

.

Trends Cogn Sci

.

7

:

553

559

.

Wilson

HR

,

Cowan

JD

.

1972

.

Excitatory and inhibitory interactions in localized populations of model neurons

.

Biophys J

.

12

:

1

24

.

Zumer

JM

,

Scheeringa

R

,

Schoffelen

J-M

,

Norris

DG

,

Jensen

O

.

2014

.

Occipital alpha activity during stimulus processing gates the information flow to object-selective cortex

.

PLoS Biol

.

12

:

e1001965

.

Author notes

Stephen Coombes and Matthew J. Brookes contributed equally.

© The Author(s) 2018. Published by Oxford University Press. All rights reserved. For Permissions, please e-mail: [email protected]

Supplementary data