Relationships Between Neuronal Oscillatory Amplitude and Dynamic Functional Connectivity (original) (raw)
Journal Article
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:
Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK
Search for other works by this author on:
Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK
Search for other works by this author on:
School of Mathematical Sciences, University of Nottingham, Nottingham, UK
Search for other works by this author on:
Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK
Search for other works by this author on:
School of Psychology, University of Nottingham, University Park, Nottingham, UK
Search for other works by this author on:
Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy, University of Nottingham, Nottingham, UK
Search for other works by this author on:
School of Mathematical Sciences, University of Nottingham, Nottingham, UK
Search for other works by this author on:
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.
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).
- Fifteen healthy volunteers (9 males, age 25 ± 4 years (mean and standard deviation)) took part in a self-paced motor task (O’Neill et al. 2015, 2017). Participants were positioned supine in the MEG system and instructed to press a button with the index finger of their left hand, approximately once every 30 s but without counting the periods between button presses. Data were recorded for 25 min giving an average of 34 button presses (trials). Button presses were recorded by a response pad. Single trials were defined as 30 s in length, 15 s before each button press, and 15 s after.
- Thirty-one healthy volunteers (18 males, age 27.4 ± 6.4 years (mean and standard deviation)) took part in a resting-state recording (Tewarie, Bright et al. 2016; Tewarie, Hillebrand et al. 2016). Subjects were positioned supine in the MEG system and were asked to “think of nothing” whilst 300 s of MEG data were acquired. Subjects fixated on a red cross throughout the acquisition. For subsequent analysis, the length of a trial was the whole recording (excluding artefacts).
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:
- Functional connectivity for each window was compared with that measured in the subsequent window. This was done at the whole brain level using a Wilcoxon rank test. Statistical testing was also performed at the regional level using a nonparametric permutation test, where a null-distribution for between window differences (independent _t_-test) was derived by permuting group (window) assignment and calculating a _t_-statistic after each permutation (1000 permutations), followed by a Bonferroni correction for multiple comparisons (Nichols and Holmes 2002).
- We used 2 different methods to reconstruct surrogate data. First, surrogate data, similar in dynamics to the real data, were constructed by taking the real data for a single subject and randomizing the trial start times. We reasoned that, if no task-induced change in PDD was measurable and then the trial start times would be meaningless. In this way, static connectivity is preserved whilst fluctuations in connectivity time locked to trial onset will be averaged out. This procedure was repeated for all trials and subjects: Fifteen surrogate datasets, one for every subject, were constructed. The peak-to-peak change in PDD (beta rebound window minus beta desynchronization window) was compared between the real and surrogate distributions using a Wilcoxon rank-sum test. Second, surrogate datasets were constructed using a phase randomization procedure (Prichard and Theiler, 1994) in which the phase of the 78 regional timecourses were randomized using the same sequence of random numbers for every region (hence, maintaining the same covariance as the original data, i.e., preserving static connectivity). The same statistical tests (Wilcoxon rank-sum test) were used to compare the genuine distribution with the null-data. Note that both of these methods preserve the static functional connectivity whilst removin g any possibility of task-induced change from the surrogate data.
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.

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.

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.

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.

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.

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]