Determination of Dynamic Brain Connectivity via Spectral Analysis

Spectral analysis based on neural field theory is used to analyze dynamic connectivity via methods based on the physical eigenmodes that are the building blocks of brain dynamics. These approaches integrate over space instead of averaging over time and thereby greatly reduce or remove the temporal averaging effects, windowing artifacts, and noise at fine spatial scales that have bedeviled the analysis of dynamical functional connectivity (FC). The dependences of FC on dynamics at various timescales, and on windowing, are clarified and the results are demonstrated on simple test cases, demonstrating how modes provide directly interpretable insights that can be related to brain structure and function. It is shown that FC is dynamic even when the brain structure and effective connectivity are fixed, and that the observed patterns of FC are dominated by relatively few eigenmodes. Common artifacts introduced by statistical analyses that do not incorporate the physical nature of the brain are discussed and it is shown that these are avoided by spectral analysis using eigenmodes. Unlike most published artificially discretized “resting state networks” and other statistically-derived patterns, eigenmodes overlap, with every mode extending across the whole brain and every region participating in every mode—just like the vibrations that give rise to notes of a musical instrument. Despite this, modes are independent and do not interact in the linear limit. It is argued that for many purposes the intrinsic limitations of covariance-based FC instead favor the alternative of tracking eigenmode coefficients vs. time, which provide a compact representation that is directly related to biophysical brain dynamics.

1. Introduction

Brain activity spans many decades of spatial and temporal scale and constantly changes due to stimuli and internally generated signals (Raichle, 2011). Some of these changes are due to changed activity (e.g., evoked by stimuli) in the pre-existing brain structure, while others are due to changes in this structure that alter the activity—i.e., changes in neural connections and their strengths (Bassett et al., 2011, 2017; Deco et al., 2011; Raichle, 2011; Hutchison et al., 2013a,b; Calhoun et al., 2014; Kopell et al., 2014; Preti et al., 2017; Babaie-Janvier and Robinson, 2019, 2020). It is thus of central interest to determine how the latter dynamic readjustments of connectivity are caused by, and support, signal processing demands placed on the brain during tasks.

On learning and developmental timescales of 100 s or more, synaptic connections can be formed and broken and long-term plasticity plays a key role in optimizing brain responses (Koch, 1999; Ghosh et al., 2008; Bassett et al., 2011; Kopell et al., 2014), as do neuromodulatory effects during cognitive processes and the sleep-wake cycle (Koch, 1999; Tagliazucchi and Laufs, 2014). Effective connectivity (EC) is defined to be the strength of connections, direct or indirect, between two points. On timescales below a few seconds, changes in EC due to modification of synaptic strengths, firing thresholds, and other aspects of neural responsiveness can result from a range of biophysical processes such as adaptation and facilitation (Koch, 1999; Rennie et al., 1999, 2002; Robinson and Roy, 2015), and these affect EEG activity and evoked responses (Rennie et al., 1999; Babaie-Janvier and Robinson, 2019, 2020). There the change in connectivity can be interpreted as part of the response that implements attention (Babaie-Janvier and Robinson, 2019, 2020). In the intermediate range of roughly 5–100 s, processes such as plasticity certainly exist that can change functional brain connectivity, but their role is less understood.

What is perhaps more puzzling than the occurrence of temporal changes per se is that these often appear to be quite rapid and widespread (Britz et al., 2010; Cabral et al., 2014; Hansen et al., 2015; Preti et al., 2017; Michel and Koenig, 2018). Much of the evidence for such changes is from studies of functional connectivity (FC), which is most commonly defined to be the covariance of activity at pairs of points and is often described in terms of patterns termed resting state “networks” (Friston, 1994, 2011; Raichle et al., 2001; Bullmore and Sporns, 2009; Sporns, 2011; Zalesky et al., 2012; Hutchison et al., 2013a; Zalesky and Breakspear, 2015; Fornito et al., 2016). The basis for this definition is the assertion that positively correlated points are likely to be involved in supporting the same function, an issue to which we return later where we show that negative correlations are equally important to the dynamics (Robinson, 2019). Because FC is defined in terms of activity, it is quite possible for it to change quickly and large-scale activity patterns are observed to change on timescales as short as 50–100 ms using EEG and MEG (Britz et al., 2010; Musso et al., 2010; Van De Ville et al., 2010; Michel and Koenig, 2018), or over tens of seconds using functional MRI (fMRI) (v et al., 1995; Damoiseaux et al., 2006; Fox and Raichle, 2007; Britz et al., 2010; Chang and Glover, 2010; Deco and Jirsa, 2012; Hipp et al., 2012; Hutchison et al., 2013a,b; Calhoun et al., 2014; Mitra et al., 2014; Hansen et al., 2015; Chang et al., 2016; Cabral et al., 2017; Babaie-Janvier and Robinson, 2019; Hunyadi et al., 2019). Notably, the patterns of activity and inferred connectivity are similar, but not identical, across different states of arousal, under both spontaneous and task-based conditions, and when observed using differing measurement modalities (Peltier et al., 2005; Smith et al., 2009; Tagliazucchi and Laufs, 2014; Chang et al., 2016; Robinson, 2019).

Rapid large-scale reorganizations of FC have been reported by numerous authors working with fMRI, particularly when matching to libraries of patterns is used in its quantification (Deco et al., 2011, 2017; Calhoun et al., 2014; Hansen et al., 2015; Cabral et al., 2017). This is at least in part because matching is done moment by moment to the most similar of a finite number of patterns, which necessarily induces apparent rapid interstate jumps or “switching” even if the actual underlying dynamics are smooth. Similar sudden jumps between EEG microstates (which are characteristic large-scale patterns of scalp potential) have been asserted when a small number of microstates are used for matching (Britz et al., 2010; Gabay et al., 2018). However, the latter effect has been shown to be consistent with smooth evolution (no jumps) being projected onto the nearest discrete microstate at each point in time (Gabay et al., 2018). It would thus be advantageous to develop methods that allow for continuous evolution, with sufficient time resolution to follow rapid changes when and if they occur. The ability to discriminate, if possible, between effects of changing structure and changing activity on FC would also be valuable.

One limitation on tracking dynamic FC is that it is most often measured using fMRI, which is slow and usually involves calculating the covariance within a sliding window of tens of seconds in duration (Fox and Raichle, 2007; Bullmore and Sporns, 2009; Chang and Glover, 2010; Deco et al., 2011; Hutchison et al., 2013a,b; Allen et al., 2014; Leonardi et al., 2015; Chang et al., 2016; Hindriks et al., 2016). This limits temporal resolution and introduces artifacts due to contributions from fluctuations with periods longer than the window length (Yule, 1926; Leonardi et al., 2015; Zalesky and Breakspear, 2015; Ernst et al., 2017). This has resulted in a host of attempts to quantify FC changes and to separate true changes from windowing artifacts. These have overwhelmingly been statistically based and have involved many disjoint methods, including correlations of FC with EC or structural connectivity (SC; i.e., anatomical connectivity), abstract graph theory, clustering, PCA, ICA, matching to dictionaries of predetermined patterns and many others (Arfanakis et al., 2000; Beckmann et al., 2005; Bullmore and Sporns, 2009; Greicius et al., 2009; Honey et al., 2009, 2010; Smith et al., 2009; Pernice et al., 2011; Yeo et al., 2011; Fornito et al., 2013; Leonardi et al., 2013, 2014; Allen et al., 2014; Anderson et al., 2014; Calhoun et al., 2014; Mitra et al., 2014; Hansen et al., 2015; Yaesoubi et al., 2015; Eklund et al., 2016; Cabral et al., 2017; Preti and Van De Ville, 2019), as reviewed by a number of authors (Deco et al., 2011; Sporns, 2011; Fornito et al., 2016; Bassett et al., 2017; Preti et al., 2017). As well as many methods being ad hoc and/or assuming particular connectivity architectures such as modularity or hierarchy, the relationships between different methods and their respective results are often obscure or unknown.

We stress that phenomenological and statistical methods often yield robust differences between conditions, and can be used to narrow the range of possible links between phenomena and for classification and hypothesis testing, but their nature prevents them from providing direct links to physical mechanisms. In contrast, physically based models have increasingly exposed these links in a straightforward way, showing how EC supports activity, and how correlations and resulting FC can be calculated directly from EC (Robinson, 2012, 2019; Friston et al., 2014; Robinson et al., 2014, 2016). It is found that the total EC is equivalent to the system transfer function that relates outputs to inputs (Robinson, 2012, 2019), thereby linking such analyses directly to engineering control systems (Babaie-Janvier and Robinson, 2018, 2019, 2020). Notably, such approaches explicitly link changes in EC and/or activity with resulting FC changes, which can be especially rapid when caused by changes in activity rather than structure or gain. Moreover, work in the last decade has shown that the inverse problem of determining EC from FC can be solved in a wide variety of situations (Galán, 2008; Friston et al., 2014; Robinson et al., 2014, 2016; MacLaurin and Robinson, 2019; Robinson, 2019; Henderson et al., 2021), as can the problem of obtaining EC, and hence FC, from evoked responses (Robinson et al., 2018; Henderson et al., 2021). Figure 1 shows approximate scales spanned by EEG, fMRI, and low-order modes and highlights the fundamental difference between statistical and physical approaches, with the latter prioritizing the system that generates the activity, with observations and their correlations flowing as corollaries.

FIGURE 1

Figure 1. Schematic of brain observations and analysis. (A) Approximate spatiotemporal scales spanned by EEG, ERP, and fMRI observations, and by phenomenological resting state networks and microstates. (B) Phenomenological and statistical links between measurements and constructs. (C) Physical activity generates the links seen in (B).

Systematic physically based work on EC-FC relationships has rested on spectral analysis in which activity, EC, and FC are all expanded in series of eigenfunctions (also termed natural modes, eigenvectors, eigenmodes, or just modes) of the system, analogous to the natural modes of a violin string. This allows the dominant, large-scale activity and connectivity to be compactly expressed in terms of just a few modes, leading to dramatic simplifications in their representation. For example, the 998 × 998 FC matrix of (Honey et al., 2009) has nearly 500 000 independent elements, but can be reduced to just 998 entries that represent eigenmodes; these in turn are dominated by just a few tens of entries and can often be approximated by ~10 eigenmodes and their amplitudes in applications to EEG, evoked activity and fMRI (Nunez, 1995; Nunez et al., 2001; Robinson et al., 2014, 2016, 2018; Mukta et al., 2019, 2020; Gao and Robinson, 2020). This compact representation is the fundamental reason that fMRI and EEG experiments only robustly detect around 10 resting-state patterns or around 5 microstate patterns, respectively, against a noisy background.

Unlike statistical components extracted by some varieties of independent component analysis (ICA) or principal component analysis (PCA) (Beckmann et al., 2005; Leonardi et al., 2013; Anderson et al., 2014), or clusters of nodes grouped on the basis of similarity of correlation (Yeo et al., 2011), eigenmodes are the natural dynamical modes of the physical system, rather than statistical constructs, and thus reflect connectivity and dynamics directly (Robinson, 2012, 2019; Friston et al., 2014; Robinson et al., 2014, 2016; Gabay et al., 2018). Indeed, if the aim is to probe dynamics of the actual physical brain, rather than treat it as a black-box signal generator, one must go beyond statistics, model-free signal analysis, and phenomenology to dynamics-based methods. In symmetric cases such as those obtained by covariance-based FC, eigenmodes form a complete orthonormal basis, which means that any spatially continuous activity and connectivity whatsoever can be expressed in terms of them (Galán, 2008; Friston et al., 2014; Robinson et al., 2014; Robinson, 2019)—this includes microstates, resting state “networks” (RSNs), and ICA and PCA components (v et al., 1995; Beckmann et al., 2005; Fox and Raichle, 2007; Deco et al., 2011; Anderson et al., 2014)—a key advantage of treating the brain as a physical system.

For readers who are not familiar with eigenmodes, Figure 2A shows a simple 1D example of the first few eigenmodes of a violin string that is fixed at both ends, where displacement is the analog of changes in brain activity relative to baseline. For a uniform string, each eigenmode is a spatial sinusoid that oscillates in time between the two extremes shown; an integer number of half wavelengths must fit into the length of the string. Key points worth noting are: (i) The eigenmodes are determined by the geometry of the string and its boundary conditions—i.e., its length and fixed endpoints. The endpoints do not move (i.e., no change in activity) but are critical to the properties of the modes. (ii) Each mode occupies the whole string, not a localized region. (iii) Modes overlap in space and every point on the string is part of every mode, although some points in each mode (so-called nodes or zeros) have zero amplitude. (iv) Despite (iii), modes are independent and do not interact, unless nonlinear effects are introduced. Even in the nonlinear case, linear spatial modes remain a useful starting point for spatiotemporal analysis. (v) There is no part of the string that can be excised and said to produce one mode of oscillation. (vi) Each mode, except the lowest, involves both simultaneous positive and negative regions. Hence, anticorrelated regions are just as essential to the dynamics as correlated ones and nodal lines of zero activity change will exist. (vii) Because the lowest mode has the same sign at all points at a given instant, it tends to induce global positive correlations relative to any given point, whereas the contributions from other modes average spatially to zero, so these modes always tend to cause some other points to be anticorrelated with any given point. When this mode is deleted, negative correlations must occur, as has been reported as a result of global signal removal in fMRI (Murphy and Fox, 2017). (viii) Eigenmode expansion is illustrated in Figure 2B, which shows how a complex waveform can be expressed in terms of the amplitudes and phases of the lowest eigenmodes.

FIGURE 2

Figure 2. Eigenmodes of a violin string clamped at both ends. (A) First five modes, showing waveform at times (solid and dotted) half an oscillation apart. (B) Examples of waveforms obtained by superposition of the modes in (A) at times t₁ . . . t₅, with frequencies proportional to the mode number as in a real violin string, equal amplitudes, and random initial phases.

The features noted in the previous paragraph are based on methods that have been developed and applied over two centuries in physics and engineering, starting with Fourier (Fourier, 1822; Schiff, 1968; Zwillinger, 1989; Courant and Hilbert, 2004), and whose properties are extremely well understood. These attributes contrast with many inherent in widely used FC analyses (Robinson, 2019): (i) ICA and PCA do not take into account geometry. (ii) Resting state networks are usually constructed to be spatially nonoverlapping. (iii) Use of positive correlations to define “networks” in some analyses neglects dynamically essential negative correlations—the supposedly problematic production of negative correlations when global signals are removed from fMRI (Chen et al., 2011) is actually an automatic consequence of removing activity in the lowest mode, which is the only one with global all-positive correlations. (iv) Thresholding of correlations deletes modal zeros, which have no activity but are critical in determining dynamical properties, and selects a few local regions as being responsible for dynamics. This is akin to retaining just the highest-amplitude regions of the string in Figure 2A, which would certainly not produce a note. (v) The common practice of coarsely discretizing the brain then treating thresholded links between these regions as being an actual brain network with graph-theoretic properties such as degree, clustering, modularity, and hierarchy is a category error (Ryle, 1949) because there is no discrete cortical network with such properties at macroscopic scales. (vi) Many methods ignore the dimensionality and units of the quantities involved, leading to neglect of relative areas of regions of interest, for example.

In this paper we use spectral analysis to quantify dynamic connectivity and constrain its dependence on actual EC changes, transient stimuli, and windowing effects. Specifically, we synthesize multiple results from the literature on eigenmode analysis of brain structure and function to decompose brain activity in terms of long-term average eigenmodes of the covariance of activity. These modes form a complete orthonormal basis, which we demonstrate to be robust to perturbations—so there can be no rapid wholesale changes in these dynamical building blocks. We then use these basis functions to compactly decompose instantaneous brain activity, EC, and FC in terms of mode amplitudes. These are then used to probe short term changes in connectivity, plus the effects of windowing, averaging, natural activity bandwidth, and other phenomena. Similarly, we show how to express evoked responses and the system transfer function in terms of these modes and to approximate them via sparse measurements of evoked responses.

The structure of the paper is as follows: Section 2 briefly overviews the background theory for an interdisciplinary readership. It generalizes key results and demonstrates that eigenmodes of the corticothalamic system have stable spatial structure under moderate perturbations. Section 3 then discusses how to analyze dynamic FC in terms of eigenmodes, examines intrinsic limitations of covariance-based approaches, and points out advantages of a more direct approach via eigenmode coefficients. Examples are provided to illustrate the key results in simple situations, with a summary and discussion in section 5.

2. Materials and Methods

In this section we provide the necessary background to our analysis, generalizing and adapting it as required. We first outline neural field theory and its relationship to EC, then introduce eigenmodes and expansions of activity, EC, and FC in terms of them. This summarizes and elucidates prior work, cited below, extends prior results, expresses them in alternative notations that are useful in different contexts, and clarifies several misconceptions. To keep the scope manageable we restrict attention to corticocortical EC and FC but allow for local dynamics, which can include corticothalamic feedbacks (Robinson et al., 2016; Robinson, 2019).

2.1. Neural Field Theory

The spatial connectivity scales of interest here range from ~1 mm up to the whole brain, with large-scale patterns seen in microstates and resting state networks having sizes of several cm and up. Connectivity is continuous, not discrete, on the scales of interest (i.e., neurons are not resolved). Hence, we do not need to analyze individual neurons and neural field theory (NFT) is well suited to address activity and connectivity (Nunez, 1995; Deco et al., 2008; Coombes et al., 2014; Robinson, 2019). NFT averages over neural populations to obtain equations for local means of quantities such as afferent activity, soma potential, and firing rate, which underlie observable signals (Nunez, 1995; Jezzard et al., 2001). Some reviews of its many successful applications to predict and analyze experimental results include (Deco et al., 2008; Coombes et al., 2014).

Large scale connectivity is usually assumed to be approximately linear, which is a good approximation because axonal outputs are generally proportional to their inputs even when the points they connect undergo nonlinear dynamics. Hence, we employ the linear limit of NFT here and make as few other assumptions as possible about the specific form of the theory to keep the results as general as possible. We note that if the equations of NFT are discretized (e.g., for numerical analysis or coarse-grained approximation), this yields neural mass theory (NMT), which is often incorrectly asserted to be distinct from NFT. By proceeding from NFT to NMT one obtains correctly weighted connections between discretized spatial neural masses, which is not the case if NMs are simply added to a network such as an observed connectome.

2.1.1. Integral Formulation

We assume that experimental signals are linear functions of perturbations Q(r, t) of brain activity from an average steady state. In the linear regime, one has (Robinson, 2019).

$\begin{array}{ll}Q(r,t)={\displaystyle \int {\displaystyle \int \Lambda }}(r,t;r^{\prime },t^{\prime })Q(r^{\prime },t^{\prime })d^2{r}^{\prime }d{t}^{\prime }+N(r,t),& (1)\end{array}$

where signals travel to a position r at time t from r′ at time t′ with a strength Λ, which is the direct effective propagator or direct EC tensor, and N represents external input (Robinson, 2012, 2019; Robinson et al., 2014). Figure 3 shows the physical situation corresponding to Equation (1). In Equation (1) local dynamics, including interactions between short-range populations such as interneurons, are absorbed into the structure of Λ (Robinson, 2019), so only long range propagation via white matter fibers is explicit, as in macroscopic observations.

FIGURE 3

Figure 3. Schematic of the terms in Equation (1) for propagation of activity to r = (r,t) from r′ = (r′, t′). The total propagator T is the sum of terms representing direct input I plus, on the right, direct propagation Λ, propagation via neural interaction at one intermediate location r″, and so forth. Most commonly arrows indicate propagation between cortical locations via white matter fibers, while vertexes denote points where local interactions occur, including via corticothalamic loops.

Before proceeding, we note that the linear assumption in Equation (1) is well justified as a starting point: (i) Although the steady states of the brain may be determined nonlinearly, perturbations from them can be linearly approximated. Linear approximations have been successful in describing a wide range of normal phenomena, as noted in the Introduction. Dynamics of phenomena such as seizures require a fully nonlinear approach, but are not considered in the present work. (ii) Correlations of relatively weak activity changes suffice to establish functional connections. (iii) Spatial connectivity is often approximately linear even when local temporal dynamics are nonlinear—output spike rates from each axon closely follow the rate at which spikes are generated at its axonal hillock. (iv) Unless the linear case is understood, it is likely to be premature to try to understand nonlinear aspects of the problem.

One can re-express Equation (1) in the form

$\begin{array}{ll}Q(r,t)={\displaystyle \int {\displaystyle \int T}}(r,t;r^{\prime },t^{\prime })N(r^{\prime },t^{\prime })d^2{r}^{\prime }d{t}^{\prime },& (2)\end{array}$

where T is the system transfer function, or total EC (Robinson, 2012, 2019). Note that if N is a delta function, the response is just T, whence we note that T is the overall Green function of the system, and corresponds to the widely studied evoked response (Luck and Kappenman, 2012; Robinson, 2012, 2019; Robinson et al., 2018; Henderson et al., 2021).

In a system with static structure (we generalize this in section 3), the correlation function can be used to define a generalized FC via

$\begin{array}{ll}C(r,r^{\prime },\tau )={\displaystyle \int Q}(r,t+\tau )Q(r^{\prime },t)dt,& (3)\end{array}$

which is equivalent to selecting the zero-frequency Fourier component in the spectral domain. In practice, one would perform this integral over a window long enough that correlations between activity at opposite ends can be ignored. If the system structure does not depend on time, Λ and T depend on time only via the difference τ, so we can use the notations Λ(r, r′, τ) and T(r, r′, τ). The temporal integrals in Equations (1, 2) then become convolutions and we can Fourier transform vs. time to obtain (Robinson, 2012, 2019).

$\begin{array}{ll}Q(r,\omega )={\displaystyle \int \Lambda }(r,r^{\prime },\omega )Q(r^{\prime },\omega )d^2{r}^{\prime }+N(r,\omega ),& (4)\end{array}$

$\begin{array}{ll}={\displaystyle \int T}(r,r^{\prime },\omega )N(r^{\prime },\omega )d^2{r}^{\prime },& (5)\end{array}$

where ω is the angular frequency. Setting τ = 0 in Equation (3) yields the covariance, which is widely used to define functional connectivity in fMRI experiments.

In the commonly considered case of “resting state” activity, in which no task is imposed by an experimenter and the subject is in relaxed surroundings (although the brain is not resting but is constantly performing background functions Raichle, 2011), background stimuli span a broad range of spatial and temporal scales after passing through the peripheral nervous system, which also tends to whiten them to make best use of available bandwidth. Consequently, numerous applications to experimental data have shown that background perturbations N(r, ω) can be approximated by spatially uncorrelated white noise (Robinson et al., 1997, 2002; Deco et al., 2008; van Albada et al., 2010; Abeysuriya et al., 2015), with

$\begin{array}{ll}{\displaystyle \int \int }N(\text{r},t+\tau )N^{*}({\text{r}}^{\prime },t)d^2\text{r}dt=N_0^2{\delta }^2(\text{r}-{\text{r}}^{\prime })\delta (\tau ),& (6)\end{array}$

with normalization N₀. We then evaluate the Fourier transform of (3); i.e., the cross-spectral density:

$\begin{array}{ll}C(\text{r},{\text{r}}^{\prime },\omega )=N_0^2{\displaystyle \int }T(\text{r},{\text{r}}^{″},\omega )T^{*}({\text{r}}^{\prime },{\text{r}}^{″},\omega )d^2{\text{r}}^{″},& (7)\end{array}$

where the asterisk denotes complex conjugation and the integral is over all r″.

A key issue that is largely ignored in the literature is that most of the quantities in the above equations have physical dimensions and units. If these are not included, serious errors can occur (Robinson, 2019), such as equating quantities with different dimensions (e.g., area = volume) or failing to include relative areas when comparing different discretized regions. These dimensions are summarized in Table 1 on the assumption that Q(r, t) has the dimensions of a firing rate.

TABLE 1

Table 1. Quantities in Equations (1–14) and their SI units.

2.1.2. Matrix Approximation

It is often useful to discretize the above equations on a fine enough scale that the discretized quantities faithfully represent the continuous ones. We follow this approach as being the most convenient to present a number of key results used below; but revert to the continuous representation when required. Most importantly, we do not make the mistake of viewing this discretization as corresponding to an actual discrete macroscopic cortical network that can be analyzed via graph theory—no such network exists. Data are always recorded in discretized form, but not necessarily resolved finely enough to represent the underlying dynamics faithfully.

For now we limit attention to Equations (4–8) and define a discrete set of positions r^j which are viewed as elements of a column vector r (we use superscripts to denote vector and matrix elements to avoid confusion with subscripts, which are used below to denote different eigenvectors and their amplitudes). The corresponding Q(r, ω) and N(r, ω) are also column vectors so we write Equation (4) as (Robinson, 2019).

$\begin{array}{ll}Q(\omega )\approx \widehat{\Lambda }(\omega )Q(\omega )+N(\omega ),& (8)\end{array}$

where the elements of the square matrix $\widehat{\Lambda }$ are

$\begin{array}{ll}{\widehat{\Lambda }}^{jk}(\omega )={\Lambda }^{jk}\Delta S^k,& (9)\end{array}$

ΔS^k being the discrete piece of cortical surface area represented by r^k. The approximation (9) becomes more accurate as the discretization becomes finer; an analogous definition is used for $\widehat{T}$. Strictly, $\widehat{\Lambda }$ is a tensor because its rows and columns are indexed by 2D position, but this does not affect the analysis below. In a similar way, if we divide by the normalization $N_0^2$ of the assumed white noise to obtain a normalized version Ĉ of C, Equations (5, 7) become

$\begin{array}{ll}Q(\omega )=\widehat{T}(\omega )N(\omega ),& (10)\end{array}$

$\begin{array}{ll}Ĉ(\omega )=\widehat{T}(\omega ){\widehat{T}}^{†}(\omega ),& (11)\end{array}$

where the dagger indicates the Hermitian conjugate (Galán, 2008; Pernice et al., 2011; Robinson, 2012; Pinotsis et al., 2013; Friston et al., 2014; Robinson et al., 2014). The covariance, which is the most commonly used form of FC is given by Equation (3) with τ = 0, which is proportional to C(r, r′, ω = 0). Most commonly, the time average is subtracted from the signals before the covariance is computed and sometimes the result is normalized.

Using the above notation, Equations (8, 10) yield (Galán, 2008; Robinson et al., 2014, 2016; Robinson, 2019).

$\begin{array}{ll}\widehat{T}(\omega )={\left[I-\widehat{\Lambda }(\omega )\right]}^{-1},& (12)\end{array}$

$\begin{array}{ll}=I+\widehat{\Lambda }(\omega )+{\widehat{\Lambda }}^2(\omega )+\dots ,& (13)\end{array}$

$\begin{array}{ll}\widehat{\Lambda }=I-{\widehat{T}}^{-1},& (14)\end{array}$

where I is the unit matrix, the superscript −1 indicates the matrix inverse, and (13) applies within the radius of convergence. The coordinate equivalent of the formal expansion in Equation (13) is Equation (14) of Robinson (2019). In Equation (12), $\widehat{T}$ is the total response, whereas $\widehat{\Lambda }$ is the part of the response that travels directly from source to destination, ${\widehat{\Lambda }}^2$ has one intermediate interaction, and so forth (Galán, 2008; Adachi et al., 2012; Robinson, 2012, 2019; Mehta-Pandejee et al., 2017; Tewarie et al., 2020); $\widehat{T}$ and $\widehat{\Lambda }$ are thus total and direct ECs.

2.1.3. Perturbation Analysis

As noted above EC can include transient modifications that are due to gain changes that can be approximated as part of the linear response to stimuli (Koch, 1999; Rennie et al., 1999, 2002; Robinson and Roy, 2015; Babaie-Janvier and Robinson, 2020). We treat such effects via perturbation analysis and, in this subsection only, superscript quantities (0) for unperturbed values and (1) to denote perturbations. In this case, we can generalize Equation (8) to yield the following relationship:

$\begin{array}{ll}{Q}^{(1)}(\omega )={\widehat{\Lambda }}^{(0)}(\omega )Q^{(1)}(\omega )+{\widehat{\Lambda }}^{(1)}(\omega )Q^{(0)}(\omega )+N^{(1)}(\omega ).& (15)\end{array}$

This can be understood if we note that Equation (1) is written on the assumption of static $\widehat{\Lambda }={\Lambda }^{(0)}$, so perturbations Q⁽¹⁾ arise from the effect of ${\widehat{\Lambda }}^{(0)}$ on perturbations Q⁽¹⁾ and direct drive by the inputs N⁽¹⁾. However, if $\widehat{\Lambda }$ is perturbed by an amount ${\widehat{\Lambda }}^{(1)}$, there will be another contribution to Q⁽¹⁾ from the effect of the perturbation on the steady state activity Q⁽⁰⁾.

We now suppose that $\widehat{\Lambda }$ is affected by spatially local feedbacks due to habituation, facilitation, and other processes. In coordinate notation, we assume that changes in Λ⁽¹⁾(r, r′, τ) are driven by local activity arriving at r from r′ and obey a temporal differential equation. Then, writing indexes explicitly, one has

$\begin{array}{l}{\left[{\widehat{\Lambda }}^{(1)}(\omega )\right]}^{ij}=Z(\omega )\{{\left[{\widehat{\Lambda }}^{(0)}(\omega )\right]}^{ij}{\left[Q^{(1)}(\omega )\right]}^j\\ +{\left[{\widehat{\Lambda }}^{(1)}(\omega )\right]}^{ij}{\left[Q^{(0)}(\omega )\right]}^j\},& (16)\end{array}$

where ${Z}(\omega )$ is the relevant differential operator, which is assumed to be the same at all points, but this assumption can be relaxed. The other factor on the right is the activity arriving at r from r′, labeled i and j, respectively. Equation (16) can be solved for ${\left[{\widehat{\Lambda }}^{(1)}\right]}^{ij}$, whence (15) can be written as in Equation (8), but with

$\begin{array}{ll}\widehat{\Lambda }={\widehat{\Lambda }}^{(0)}+\widehat{{L}},& (17)\end{array}$

$\begin{array}{ll}{\left[\widehat{{L}}\right]}^{ij}=\frac{Z(\omega ){\left[{\widehat{\Lambda }}^{(0)}\right]}^{ij}{\left[Q^{(0)}\right]}^j}{1-Z(\omega ){\left[Q^{(0)}\right]}^j},& (18)\end{array}$

and the transfer function is modified to

$\begin{array}{ll}\widehat{T}={\left[I-{\widehat{\Lambda }}^{(0)}-\widehat{{L}}\right]}^{-1}.& (19)\end{array}$

Figure 4 shows an example of cortical evoked responses and transfer functions from inputs to cortical activity, with and without accompanying EC changes due to gain evolution. These curves are obtained from a full analysis whose details can be found in Babaie-Janvier and Robinson (2020), showing that evoked EC changes (and, by implication, evoked FC changes) can be substantial. This result implies that one cannot separate the contributions from ${\widehat{\Lambda }}^{(0)}$ and $\widehat{{L}}$ without a model that allows the form and dynamics of ${Z}$ to be estimated.

FIGURE 4

Figure 4. Evoked responses and transfer functions for cases with static EC (green) and dynamic EC (red) starting from the same initial state. (A) Model cortical evoked responses. The dashed black line shows the contribution of the EC dynamics and adds to the green curve to give the red one. (B) Magnitudes of the corresponding transfer functions.

Equations (11) and (12) explicitly show the connections between direct EC, total EC, and resting state FC, thereby obviating much of the huge literature devoted to probing these connections via statistical methods such as correlations between connectivity patterns or connection matrix entries (Greicius et al., 2009; Sporns, 2011; Messe et al., 2014, 2015; Fornito et al., 2016; Preti et al., 2017). Equations (15–19) show that the connections between structural connectivity (SC; i.e., anatomy) and EC are more complex because it is not only the presence of connections that is relevant, but also their time-dependent strength (Britz et al., 2010; Bassett et al., 2011; Deco and Jirsa, 2012; Kopell et al., 2014; Hansen et al., 2015; Leonardi et al., 2015; Deco et al., 2017; Babaie-Janvier and Robinson, 2019, 2020). However, because only the EC is directly relevant to the dynamics and the FC, the SC-EC relationship can be viewed as a separate issue.

2.2. Eigenmode Expansion: Spectral Analysis

The matrix Ĉ in Equation (11) is Hermitian, so its spatial eigenfunctions u_j form an orthonormal basis set (Schiff, 1968; Courant and Hilbert, 2004), and are written as column vectors that approximate the continuum eigenmodes and satisfy

$\begin{array}{ll}Ĉ(\omega )u_j={\kappa }_j(\omega )u_j,& (20)\end{array}$

for some eigenvalues κ_j(ω), which are discussed below. Note that in a static structure the u_j do not depend on time, just as for the shapes of the sine waves in Figure 2.

Analogous equations to (20) apply for $\widehat{\Lambda }$ and $\widehat{T}$. Given Equations (12, 13) these quantities commute and thus have the same eigenvectors, v_j, with

$\begin{array}{ll}\widehat{\Lambda }(\omega )v_j={\lambda }_j(\omega )v_j,& (21)\end{array}$

$\begin{array}{ll}\widehat{T}(\omega )v_j={\theta }_j(\omega )v_j.& (22)\end{array}$

If the connectivities are symmetric [i.e., Λ(r, t, r′, t′) = Λ(r′, t, r, t′)] then v_j = u_j, so the eigenmodes of the white-noise-driven covariance are the same as those of the transfer function and the activity. We assume symmetry for the remainder of the present work.

Because the eigenfunctions u_j form a complete orthonormal set, one can expand arbitrary activity in terms of them (Schiff, 1968; Courant and Hilbert, 2004); i.e.,

$\begin{array}{ll}Q(\omega )={\displaystyle \sum _j}{c}_j(\omega )u_j,& (23)\end{array}$

$\begin{array}{ll}{c}_j(\omega )=Q^T(\omega )u_j,& (24)\end{array}$

$\begin{array}{ll}{u}_j^T{u}_k={\delta }_{jk},& (25)\end{array}$

where δ_jk is the Kronecker delta. Alternatively, in coordinate notation, one has

$\begin{array}{ll}Q(\text{r},t)={\displaystyle \sum _j}{c}_j(t)u_j(\text{r}),& (26)\end{array}$

$\begin{array}{ll}{c}_j(t)={\displaystyle \int }{u}_j^{*}(\text{r})Q(\text{r},t)d^2\text{r},& (27)\end{array}$

$\begin{array}{ll}{\delta }_{jk}={\displaystyle \int }{u}_j^{*}(\text{r})u_k(\text{r})d^2\text{r}.& (28)\end{array}$

One can also expand the connectivities in terms of eigenfunctions via spectral decomposition. In the symmetric case, one has

$\begin{array}{ll}\widehat{\Lambda }(\omega )=UL(\omega )U^{†},& (29)\end{array}$

$\begin{array}{ll}\widehat{T}(\omega )=U\Theta (\omega )U^{†},& (30)\end{array}$

$\begin{array}{ll}Ĉ(\omega )=UK(\omega )U^{†},& (31)\end{array}$

where U is the unitary matrix whose columns are the orthonormal eigenvectors and L, Θ, and K are diagonal matrices whose nonzero entries are the corresponding eigenvalues Robinson et al. (2014); Robinson (2019), which satisfy

$\begin{array}{ll}{\kappa }_j(\omega )=|{\theta }_j(\omega ){|}^2,& (32)\end{array}$

$\begin{array}{ll}{\theta }_j(\omega )={\left[1-{\lambda }_j(\omega )\right]}^{-1},& (33)\end{array}$

using Equations (11, 12), respectively. From Equations (10, 30) we have

$\begin{array}{ll}N={\widehat{T}}^{-1}Q=U{\Theta }^{-1}{U}^{†}Q,& (34)\end{array}$

which describes how to infer the input from measurements of the activity using T. Terms in the power series in Equation (13) have the form [Λ(ω)]^m = U[L(ω)]^mU^†.

The coordinate space equivalents of Equations (29–31) are

$\begin{array}{ll}\Lambda (r,r^{\prime },\tau )={\displaystyle \sum _j{u}_j^{\ast }}(r)u_j(r^{\prime }){\lambda }_j(\tau ),& (35)\end{array}$

$\begin{array}{ll}T(r,r^{\prime },\tau )={\displaystyle \sum _j{u}_j^{\ast }}(r)u_j(r^{\prime }){\theta }_j(\tau ),& (36)\end{array}$

$\begin{array}{ll}C(r,r^{\prime },\tau )={\displaystyle \sum _j{u}_j^{\ast }}(r)u_j(r^{\prime }){\kappa }_j(\tau ),& (37)\end{array}$

where τ = t − t′. Each eigenvector u_j, or u_j(r) in coordinate notation, extends over the whole system, with every point being part of every eigenfunction, as in the example in Figure 2. However, the diagonal representations in Equations (29–31) mean that these discrete modes are not coupled, are completely independent in the linear regime, and can only be linearly excited by input that has a component in the same mode. Eigenmodes thus form a fundamentally discrete set of dynamical building blocks, unlike statistically derived and/or phenomenological resting state patterns, for example.

Yeo et al. (2011) and others have clustered discretized brain regions based on their overall patterns of functional connectivity defined by the covariance. Their analyses also employed thresholding and sparsification of connection matrices, which we do not do here. Rather, we show that their measure of similarity of connectivity can be expressed in simple form using the present methods. The similarity S(r₁, r₂) of FC patterns emanating from r₁ and r₂ is the dot product between the relevant vectors of functional connections:

$\begin{array}{ll}S({\text{r}}_1,{\text{r}}_2)={\displaystyle \int }C({\text{r}}_1,{\text{r}}^{\prime })C({\text{r}}^{\prime },{\text{r}}_2)d{\text{r}}^{\prime },& (38)\end{array}$

where we note that the covariance is symmetric and omit the argument ω = 0. In matrix notation, the integral corresponds to multiplication, so

$\begin{array}{ll}Ŝ={Ĉ}^2,& (39)\end{array}$

and the eigenvectors of Ŝ are those of Ĉ, while its jth eigenvalue is just ${\kappa }_j^2$.

The arrows in the schematic Figure 5 show how the various quantities above, and many others used in the field, are related. Note that phenomenological statistical approaches are located toward the circumference while more detailed physical ones are toward the center.

FIGURE 5

Figure 5. Schematic of some of the relationships of brain data, analyses, and applications. Clockwise from the top, categories of brain observables are shown, followed by measurements, direct and indirect data analyses, resulting models, and applications. In general, smaller radii refer to more dynamic and physical aspects, while larger radii relate to more static and phenomenological aspects. (Of course, some physical aspects are static and some dynamical analyses are phenomenological.) Arrows show a subset of the relationships used in the field.

2.3. Spatiotemporal Structure

So far, we have not introduced any particular theory of the dynamics of brain activity. To predict the spatial structure of modes, one needs to introduce a specific model and geometry; although modes can be inferred from FC data without doing this, via Equations (29–33) (Robinson et al., 2014; Robinson, 2019). We begin with the illustrative case of a violin string from Figure 2, then proceed to a spherical-cortex case, which forms the basis for many of the subsequent arguments where we generalize to a convoluted cortex. The methods used are standard approaches in the physics and engineering literature, but have not yet been widely used in connectomics; more details can be found in standard references such as Schiff (1968), Courant and Hilbert (2004), and Ogata and Yang (1970).

2.3.1. Violin String Analog System

The equation of motion of an undamped, undriven violin string in response to an impulse stimulus at x = 0 and t = 0 is

$\begin{array}{ll}\frac{{\partial }^{2}Q(x,t)}{\partial t^2}-v^2\frac{{\partial }^{2}Q(x,t)}{\partial x^2}=\delta (x)\delta (t),& (40)\end{array}$

where x is the position on the string, v is the velocity of waves, and Q here represents the displacement perpendicular to the string, which is the analog of brain activity.

The transfer function is the response to an impulse stimulus (Zwillinger, 1989; Courant and Hilbert, 2004). Hence T is obtained by Fourier transforming Equation (40) vs. position and time to give

$\begin{array}{ll}T(k,\omega )\propto \frac{1}{k^2{v}^2-{\omega }^2},& (41)\end{array}$

where k is the wave number, which is the spatial analog of ω.

Because the properties of the violin string do not change in time, we can solve Equation (40) by the method of separation of variables, in which we make the ansatz that Q(x, t) is a sum of solutions of the form

$\begin{array}{ll}Q(x,t)=u(x)\theta (t).& (42)\end{array}$

By substituting the form (42) into Equation (40) and dividing by Q(x, t), we find

$\begin{array}{ll}\frac{1}{v^2\theta (t)}\frac{{\partial }^2\theta (t)}{\partial t^2}=\frac{1}{u(x)}\frac{{\partial }^{2}u(x)}{\partial x^2}=-k^2,& (43)\end{array}$

where k is a separation constant that arises because the other two members of Equation (43) are independent of x and t, respectively, so they must both equal a common constant.

Taking the spatial part of Equation (43) yields the Helmholtz equation and we find

$\begin{array}{ll}u(x)=acos(kx)+bsin(kx),& (44)\end{array}$

for some constants a and b. The only orthonormal solutions that are zero at x = 0, L are

$\begin{array}{ll}{u}_j(x)={(2/L)}^{1/2}sin(k_{j}x),& (45)\end{array}$

where L is the length of the string. The permitted solutions correspond to the eigenvalues

$\begin{array}{ll}{k}_j=j\pi /L,& (46)\end{array}$

where j = 1, 2, 3, ….

Once the eigenvalues have been determined, they can be substituted into the temporal part of Equation (43), the dispersion equation, to yield

$\begin{array}{ll}{\theta }_j(t)=sin({\omega }_{j}t+{\psi }_j),& (47)\end{array}$

$\begin{array}{ll}{\omega }_j=k_{j}v,& (48)\end{array}$

where ψ_j is the phase at t = 0. Hence, the displacement corresponding to a given eigenfunction has a profile that oscillates sinusoidally in time as well as space. The coefficients of the modes can be determined from Equation (27) at t = 0. The standing waves given by Equations (42, 44, 47) can equally well be interpreted as superpositions of equal-amplitude left- and right-propagating traveling waves.

This simple example epitomizes the fact that boundary conditions—where the displacement is zero in this case—break the symmetry to select a discrete set of eigenfunctions with specific spatial structure. All modes except the lowest have both positive and negative regions, so methods that seek patterns on the basis of clustering around points of high positive correlation ignore fundamental features of the actual brain dynamics, just as a children's seesaw without both its “negatively correlated” ends cannot function. It is also worth stressing that the 1D Laplacian operator ∂²/∂x² appears in Equation (40), which is a wave equation. A common error in the literature is to presume that the appearance of this operator means that the dynamics is diffusive Atasoy et al. (2016, 2018); although the spatial eigenfunctions are the same in both cases, the dynamics are very different and must not be confused. Indeed the diffusion equation and the wave equation are in different families of partial differential equations—parabolic and hyperbolic, respectively (Zwillinger, 1989; Courant and Hilbert, 2004). Spatial spreading of activity via wave propagation can also be treated by matrix methods that include all possible multistep paths with appropriate time delays (Galán, 2008; Robinson et al., 2016; Mehta-Pandejee et al., 2017; Robinson, 2019; Tewarie et al., 2020).

2.3.2. Spherical Case

Many common neural field brain models have a transfer function of the form

$\begin{array}{ll}T(\text{k},\omega )=\frac{A(\omega )}{k^2{v}^2+q^2(\omega )v^2},& (49)\end{array}$

which generalizes Equation (41) to embody an overall frequency envelope A(ω) plus Laplacian spatial coupling and complicated local dynamics described by q(ω), which can involve multiple populations of cortical and thalamic neurons (Nunez, 1995; Robinson et al., 2002; Deco et al., 2008; Robinson, 2019); the details of A(ω) and q(ω) do not concern us here but can be found in the references cited. It is worth noting that the dynamics described by this type of transfer function have been widely verified against experiment (see Jirsa and Haken, 1996; Robinson et al., 1997, 2002, 2014, 2016; Deco et al., 2008 and the references cited therein) and the spatial coupling implied by the Laplacian operator closely matches that seen in the brain (Nunez, 1995; Robinson et al., 1997, 2016; Nunez et al., 2001; Mukta et al., 2019, 2020) and yields eigenfunctions that are very similar to those obtained from connectivity matrices and anatomical studies (Robinson et al., 1997, 2016; Braitenberg and Schüz, 1998). Moreover, the topology of a sphere is the same as that of a brain hemisphere, aside from the 0.5% lacuna where the corpus callosum passes, so it forms a useful starting point for mapping and perturbation analysis (Jirsa et al., 2001; Robinson et al., 2016; Gabay et al., 2018).

To obtain spatial eigenmodes on a sphere of radius R_s we introduce the usual spherical coordinates ϑ and φ and separate variables in these coordinates and time. This yields discrete spatial eigenmodes u_j on the sphere. In this case, the u_j are termed spherical harmonics and written in the notation Y_lm, with

$\begin{array}{ll}{Y}_{lm}(\vartheta ,\phi )=c_{lm}{P}_l^{\left|m\right|}(cos\vartheta )e^{im\phi },& (50)\end{array}$

$\begin{array}{ll}{c}_{lm}={\left[\frac{(l-|m|)!(2l+1)}{4\pi (l+|m|)!}\right]}^{1/2},& (51)\end{array}$

where l = 0, 1, 2, … and m = −l, −l + 1, …, l, and $P_l^{\left|m\right|}$ is an associated Legendre function (Dunster, 2010). Real-valued eigenmodes $Y_l^m$ can be constructed as linear combinations of Y_lm and Y_{l, −m} that give sine and cosine dependences on φ (Dunster, 2010); viz,

$\begin{array}{l}{Y}_l^m={(-1)}^m(Y_{lm}+Y_{lm}^{*})/\sqrt{2},m>0;\\ =Y_{lm},m=0;\\ =-i{(-1)}^m\left(Y_{l|m|}-Y_{l|m|}^{*}\right)/\sqrt{2},m<0.& (52)\end{array}$

The first nine of these modes are shown in Figure 6 in real form, which is more convenient for plotting and interpretation, although the complex form in Equation (50) is easier for the analytic work below. Note that symmetry implies that these modes can be rotated through any angle on the sphere and still remain eigenmodes; however, nonuniformity breaks this symmetry and removes this degeneracy.

FIGURE 6

Figure 6. First nine real eigenmodes $Y_l^m$ on a spherical cortex. In each case the mode is normalized to a maximum amplitude of unity and regions of opposite sign are shown red and blue as in the color bar.

The corresponding eigenvalues are

$\begin{array}{ll}{k}_{lm}^2=l(l+1)/R_s^2,& (53)\end{array}$

which are independent of m because of spherical symmetry. Setting the denominator of T(k, ω) to zero implies that the wave dispersion equation is

$\begin{array}{ll}{k}_{lm}^2+q^2(\omega )=0.& (54)\end{array}$

Unlike Equation (49), when Equation (54) is solved for ω it does not usually give a unique solution for a given k_lm because of the complex dynamics it embodies; rather, for each mode it typically gives rise to a series of frequency resonances, to which we return below (Robinson et al., 2002; Robinson, 2003; Babaie-Janvier and Robinson, 2018, 2019; Gabay et al., 2018). Oscillations at the resonant frequency cause the modes to change back and forth in sign as time passes.

Although the functional form of the spatial eigenfunctions is more complex on a sphere than a violin string, the same key aspects remain: they are discrete with eigenvalues determined by the boundary conditions (a spherical topology and geometry here), overlap in space because each extends over the entire cortex, and have essential anticorrelated regions.

2.3.3. Convoluted Cortex

Solution of the Helmholtz equation on the highly convoluted cortex cannot be done in closed analytical form because of the complex geometry, although the topology of a cortical hemisphere still remains that of a sphere if one neglects small gaps such as where the corpus callosum passes through, as we justify below. Numerical solution of the NFT equations yields the real modes $y_l^m$ shown in Figure 7. An appropriate choice of coordinates yields $y_l^m=Y_l^m$ in the spherical limit, although more generally, the perturbed eigenfunction $y_l^m$ is a more general linear combination of the $Y_{l^{\prime }}^{m^{\prime }}$ dominated by l′ = l (Robinson et al., 2016). Alternatively, solution of the connection matrix $\widehat{\Lambda }$ obtained from diffusion imaging yields almost identical spatial structure (Robinson et al., 2016). Unlike the spherical case, the eigenvalues k_lm are not degenerate because spherical symmetry has been broken by the convolutions, which affect the Laplacian (Jirsa et al., 2001; Robinson et al., 2016; Gabay and Robinson, 2017).

FIGURE 7

Figure 7. First nine real eigenmodes on a convoluted cortical hemisphere. In each case they are labeled $y_l^m$ by analogy with the $Y_l^m$ and normalized to a maximum value of unity; regions of opposite polarity are colored red and blue as shown in the color bar.

Much insight into the symmetry breaking and resulting effects on eigenvalues has been obtained by treating the effects of the convolutions on the Laplacian as perturbations of the spherical case and applying standard perturbation theory from quantum physics to calculate the convoluted eigenmodes as perturbed versions of the spherical ones, with the perturbations expressed as sums over the unperturbed modes (Schiff, 1968; Gabay and Robinson, 2017). This proceeds by using the standard Freesurfer mapping (Fischl et al., 1999; Robinson et al., 2016; Gabay and Robinson, 2017) to map the convoluted surface to the sphere (and vice versa) so as to assign spherical coordinates to every point for the purposes of calculation. The Laplacian on the convoluted surface is likewise mapped to the sphere and the difference from the corresponding spherical value is used as a perturbation.

Figure 8 shows the curvature on one cortical hemisphere plotted as the ratio of mean curvature κ_O of an average cortical template surface (Glasser et al., 2016), which is somewhat less folded than an actual individual cortex, vs. the curvature κ_s = 0.148 cm⁻¹ of a sphere of equal area, where the curvature is defined to be the mean of the two principal curvatures at each point (i.e., the mean of the reciprocals of the two principal radii of curvature). Despite the large local values of this ratio, and hence of the Laplacian, the perturbation analysis gave excellent results and showed that the form of the lowest modes is locked in by the relatively mild curvature perturbations due to the gross shape of the brain, as seen in Figure 9, rather than by the more severe local perturbations due to highly curved sulci and gyri (Gabay and Robinson, 2017). The reason for this can be seen from the structure of the perturbation terms that arise in the calculations, as we discuss in the next paragraphs.

FIGURE 8

Figure 8. Ratio κ_O/κ_s of one cortical hemisphere's local curvature to that of an equal-area sphere vs. position, as indicated by the color bar. Surfaces are of an average cortical template from the Human Connectome Project (Glasser et al., 2016). (A) Lateral view of the hemisphere. (B) Medial view. (C) Spherical projection corresponding to (A), mapped to spherical coordinates via the standard Freesurfer mapping (Fischl et al., 1999). (D) Spherical projection corresponding to (B). (E) Sinusoidal projection of the whole hemisphere. (F) Flat projection of the hemisphere (Fischl et al., 1999).

FIGURE 9

Figure 9. Second to ninth real eigenmodes vs. cortical folding, with each row corresponding to one mode, as labeled. Moving from left to right corresponds to the cortex being gradually folded toward its final shape, starting from an ellipsoid. In each case the modes are calculated on the surface using the formulas of Gabay and Robinson (2017) with positive and negative regions distinguished by color.

If we expand perturbations of ∇² from its spherical form as sums of Y_lm in the spherical coordinates Ω = (ϑ, φ) defined by the Freesurfer mapping, we can write the perturbation as (Gabay and Robinson, 2017).

$\begin{array}{ll}{\nabla }_{\text{pert}}^2(\Omega )={\displaystyle \sum _{l^{\prime }{m}^{\prime }}}{a}_{l^{\prime }{m}^{\prime }}{Y}_{l^{\prime }{m}^{\prime }}(\Omega ),& (55)\end{array}$

where the $a_{l^{\prime }{m}^{\prime }}$ are angular differential operators; resulting perturbations to the eigenvalues $k_{lm}^2$ and eigenfunctions Y_lm involve sums over terms with coefficients proportional to Gabay and Robinson (2017) and Schiff (1968).

$\begin{array}{ll}{\displaystyle \int }{Y}_{l_1{m}_1}^{*}(\Omega )Y_{l^{\prime }{m}^{\prime }}(\Omega )Y_{lm}(\Omega )d\Omega ,& (56)\end{array}$

which determine the contribution of mode Y_l1m1 to the perturbed state [see Equations (41, 42) of Gabay and Robinson, 2017 for details]. This integral is related to the Clebsch-Gordan coefficients of quantum theory (Schiff, 1968; Maximon, 2010); for it to be nonzero, perturbations must satisfy

$\begin{array}{ll}{m}^{\prime }=m_1-m,& (57)\end{array}$

$\begin{array}{ll}|l_1-l|\le l^{\prime }\le l_1+l& (58)\end{array}$

This means that only perturbations that satisfy (Equations 57, 58) can affect the eigenfunctions or their eigenvalues. Such selection rules are familiar from conservation of angular momentum in quantum mechanics (Schiff, 1968) and mean that the perturbed lm mode can only involve a contribution from Y_l1m1 if perturbations have a component $Y_{l^{\prime }{m}^{\prime }}$ that satisfies Equations (57, 58). Hence, because activity and connectivity patterns are dominated by low-order modes with l ≲ 3, only low-order perturbations can perturb them significantly. Indeed, it has been found that the dominant changes to low-order modes are contributed by modes with the same l (Robinson et al., 2016; Gabay and Robinson, 2017), which corresponds to coordinate rotation in the spherical case. Figure 7 shows the first nine real eigenmodes calculated on the convoluted cortex.

In Gabay and Robinson (2017) it was shown that eigenfunctions with small l and m in the spherical limit are pinned to specific orientations by large-scale curvature of the cortex on similar ~10 cm scales to the spatial structure of the eigenfunctions themselves, whereas the strong, short-scale curvature due to gyri and sulci on scales of ~2 cm had little effect, despite being larger in amplitude; this is consistent with the arguments in the previous paragraph. Departures from the mean spherical value are large, as seen in Figure 9, but the perturbation matrix element remains small because of the mismatch of spatial scales (Gabay and Robinson, 2017). The circumference of one brain hemisphere is ~60 cm and the $Y_l^m$ have l “wavelengths” (i.e., the typical separation of successive maximums) in this distance, so a 2 cm scale corresponds to l ≈ 15 if it represents a half-wavelength. The lack of effect reflects the selection rules discussed in the previous paragraph and that the strong local curvature of sulci and gyri contributes only slightly to low-order modes of the perturbation.

The number of spherical eigenmodes up to and including the level l is (l + 1)². Hence, to account for the typical numbers of <10 microstates or RSNs, we need to focus on eigenmodes with l ≲ 2; activity is also dominated by these modes (Nunez, 1995; Nunez et al., 2001; Robinson et al., 2001, 2016, 2018; Mukta et al., 2017, 2019). Note that when the bihemispheric brain is considered, the number of modes is doubled by including combinations of the present modes that are symmetric and antisymmetric between hemispheres (Robinson et al., 2016).

Before leaving this section, it is worth stressing that: (i) Even though the eigenmodes on the convoluted cortex are complicated, they can be used perfectly well to expand activity and connectivity via the general formulas in section 2.2, with sums and integrals performed numerically. (ii) The assumption of Laplacian coupling in this section is not necessary; eigenmodes found from measured FC matrices are almost identical because the underlying Green function of propagation is close to the anatomical connectivity (Robinson et al., 1997, 2016; Braitenberg and Schüz, 1998) and analogs of the above analysis in matrix notation are straightforward. We return to this point in section 3. (iii) Activity-based eigenmodes assume only linearity and automatically include the effects of unresolved short-range and subcortical dynamics and connectivity (Robinson, 2019; Gao and Robinson, 2020). (iv) These results imply that low-order modes can compactly represent complex connectivity (Robinson et al., 2016; Gao and Robinson, 2020). For example, Figure 10 shows the contributions to approximation of an experimental Ĉ from the first few modes, and the cumulative approximations obtained by summing them. This demonstrates that as few as 4 or 5 modes can give a good representation of the main features of the FC. One might wonder how the intricate structure, with block-like features, seen in Figure 10 can arise from just a few, smooth, large-scale modes. The reason is that most of the visual appearance of intricacy is an artifact of the mapping of the 2D cortex onto a 1D list of cortical locations (Henderson and Robinson, 2011). (v) Comparison of Figures 10B,C shows that these two modes are almost complements of one another. Hence, if the relative instantaneous amplitudes of activity in these modes were to change so that first one dominated then the other, the FC structure would shift substantially, leading to “switching” behavior when mapped to the instantaneous nearest of a finite set of patterns (Hansen et al., 2015).

FIGURE 10

Figure 10. Modal contributions to FC in connection-matrix form, with scale as in the color bar. (A–E) Contributions from each of the first 5 modes in order. (F–J) Cumulative sums of the frames on the left. (K) Exact FC.

It is worth noting that eigenfunctions are not direct equivalents of microstates or RSNs. Microstates are scalp projections of cortical activity, filtered through overlying tissues, so they result from activity in superpositions of multiple eigenstates (Nunez, 1995; Gabay et al., 2018). Likewise, if RSNs are defined to be spatially nonoverlapping with sharp edges and are decomposed into eigenmodes without fully taking into account the processing steps that yield them and the properties of spectral transform (Atasoy et al., 2016), they have a broad eigenspectrum due to their sharp boundaries rather than intrinsic structure in the dynamics. Hence, recent claims to have “explained” RSNs with eigenmodes (Atasoy et al., 2016, 2018) are incorrect per se, and trivial if all that is meant is that it is possible to expand RSNs in a series of eigenmodes of a structural connectivity matrix—this is mathematically guaranteed because these eigenmodes form a complete orthonormal basis that can be used to expand any well-behaved function (Schiff, 1968; Courant and Hilbert, 2004).

2.4. Control Systems Links

If eigenvalues are distinguished by a subscript j, we can approximate the part of T that arises from the jth mode as Babaie-Janvier and Robinson (2018), Babaie-Janvier and Robinson (2019), and Ogata and Yang (1970).

$\begin{array}{ll}T(k_j,\omega )={\displaystyle \sum _p}\frac{r_{jp}}{\omega -{\omega }_{jp}},& (59)\end{array}$

$\begin{array}{ll}{\omega }_{jp}={\Omega }_{jp}-i{\gamma }_{jp},& (60)\end{array}$

where Ω_jp is the real part of the frequency of pole jp, γ_jp is its damping rate, and r_jp is its weight in the expansion. The reality condition on the Fourier transform implies that any pole with nonzero Ω_jp must be paired with one at −Ω_jp with the same damping rate, and a weight that is the complex conjugate of r_jp.

The form (59) can be interpreted as a sum of PID (proportional-integral-derivative) filters (Ogata and Yang, 1970; Babaie-Janvier and Robinson, 2018), one at each pole, having the standard control-system roles of implementing prediction and attention via gain control (Babaie-Janvier and Robinson, 2018, 2019, 2020). Each pole corresponds to a resonance in the response and an enhancement in the EEG spectrum (Babaie-Janvier and Robinson, 2018, 2019, 2020), although these only yield distinguishable spectral peaks if the damping is small.

Many authors have noted similarities in the patterns of activity and connectivity observed at different timescales, ranging from the ≲ 0.1 Hz of fMRI to tens of Hz in EEG (Honey et al., 2007; Mantini et al., 2007; Ghosh et al., 2008; Van De Ville et al., 2010; Hipp et al., 2012; Marzetti et al., 2013; Gohel and Biswal, 2015; Vidaurre et al., 2018; Hunyadi et al., 2019). This is easily explained by the fact that activity at all frequencies has the same spatial eigenmodes (Robinson et al., 2016). Corticothalamic NFT has previously explained how multiple frequency resonances exist for each mode (unlike in a violin string, where there is only one frequency per mode) due to the effects of delays in loops between cortex and thalamus (Robinson et al., 2002, 2016; Deco et al., 2008; Gabay and Robinson, 2017). Likewise, the connection between weakly damped resonances and long-range correlations is a straightforward consequence of the physics—weakly damped waves are easily excited and propagate long distances (Robinson, 2003), plausibly supporting long-range in-phase communication via the communication-through-coherence hypothesis (Fries, 2005).

In applications to date it has been found that in the adult human corticothalamic system the main poles have frequencies in the traditional slow/delta, alpha, and beta bands at around 0, 10, and 20 Hz, respectively (Gabay et al., 2018). However, different modes have frequencies differing by ~1 Hz in the alpha band, for example. This means that the modes can beat with one another, giving rise to linear and circular polarization dynamics and directly explaining a number of features and selection effects in measurements of EEG microstates (Gabay et al., 2018) without needing the nonlinear mechanisms invoked by some authors (Roberts et al., 2019).

3. Results

In the previous section we reviewed and adapted the background theory and key selection rules for perturbations to affect eigenmodes. In particular, we noted that some connectivity changes can be viewed as responses to stimuli and can thus be absorbed into Λ and T (Rennie et al., 1999; Robinson and Roy, 2015; Babaie-Janvier and Robinson, 2019, 2020). We also showed that even strong perturbations that do not satisfy selection rules analogous to Equations (57, 58) do not significantly affect eigenmodes. In the present section we thus use the eigenmodes of the long-term average FC as a robust basis for expanding all dynamical quantities, using real eigenmodes without loss of generality. We then argue that these spatial modes can only be very slightly perturbed by short-term brain state changes because they are strongly pinned by overall brain curvature and any other permanent spatial parameter variations that are present, as seen in Figure 8, and because the selection rules (Equations 57, 58) prevent significant perturbations by local effects; so amplitude changes of activity in these modes are expected to dominate FC dynamics. Finally, we systematically analyze short-term dynamics by expansion in terms of the robust long-term eigenmodes.

Having laid this groundwork, we now turn our attention to treating the effects of dynamic changes in brain EC that cannot be captured as stimulus responses via equations such as Equation (16) and which are not part of the long-term (typically hours or more) structure of the brain. These putative changes have been argued to be responsible for at least part of the observed temporal evolution of FC on timescales from tens to hundreds of seconds, as observed via fMRI (Hindriks et al., 2016; Bassett et al., 2017). We do not concern ourselves here with the important and complex issues of removal of artifacts from signals—e.g., those due to head motion, breathing, and heartbeat for fMRI, and electrical interference for EEG, except to note that we assume that we start with the signals that stem from neural activity, including any global component, which will have its largest projection into the lowest mode. If the latter mode is removed, as is commonly done (examples can be found in Sporns, 2011; Fornito et al., 2016), the FC matrix will tend to develop more negative entries, as has been widely observed, because the lowest mode is the only one without negative correlations between points.

The most common approach to tracking changes in FC is to calculate covariances of temporally windowed fMRI time series. As noted in the introduction, windowed covariances have long been known to give rise to artifacts on timescales of the order of the window length (Yule, 1926; Ernst et al., 2017). The natural bandwidth (Δf ~ 0.07 Hz) of spontaneous fluctuations in fMRI signals, even with static EC, also gives rise to time variations in correlations on a timescale of ~15 s (Pang and Robinson, 2019). When trying to distinguish the effects of true EC changes we must allow for these complications and also set aside evoked changes on time scales of ≲ 15 s, which accords with observations of correlation times and simulations (Zalesky and Breakspear, 2015; Preti et al., 2017).

Numerous methods, discussed in the Introduction, have been proposed to try to deal with the above issues. Many are based on blind signal analysis and/or statistics, while others are ad hoc and often threshold and sparsify data in the interests of simplicity, but without clear justification; few have any regard to the physical nature of the system that generates the signals that are correlated to generate the FC. In contrast, the discussion of eigenmodes in section 2 shows that the physical nature of the dynamics cannot be disregarded without sundering dynamical links that one needs to retain, as illustrated schematically in Figure 1.

Many of the problems with current methods derive from the intrinsic loss of time resolution imposed by windowing—shorter windows are more affected by random variations in noise, whereas long windows lose time resolution. This is a fundamental limit imposed by Nyquist's theorem or the equivalent uncertainty principle and cannot be circumvented by any method. However, the relative effects can be estimated and one can seek to distinguish and minimize averaging effects in a systematic way.

3.1. Long-Term Eigenmodes

Stable eigenmodes that are robust with respect to transient system fluctuations are advantageous as the basis for spectral analysis. Ideally, these should be based on activity measures, such as the covariance, to ensure that the effects of spatially unresolved dynamics and connectivity are included, but estimation of these modes can be affected by noise. Modes based on solution of the Helmholtz equation on the cortical surface have also been shown to be a reasonable approximation, and are not affected by noise, but do not incorporate short-range effects nor spatial variations in parameters other than the curvature. We illustrate both types here but do not address experimental data processing issues.

fMRI signals reflect neural activity that drives changes in blood flow and oxygenation. The hemodynamics of these processes restricts the spectrum of these signals to ≲ 0.1 Hz, with a plateau below about 0.07 Hz, where most of the energy is located (Jezzard et al., 2001; Robinson et al., 2006; Drysdale et al., 2010; Pang and Robinson, 2019). Hence, if we wish to calculate mean FC in which the effects of these fluctuations on ~30-s timescales are averaged out, we need to use a window of length t_L with t_L ≳ 300 s, but not so great that state changes due to drowsiness or falling asleep occur, because these can give rise to significant alterations in FC (Tagliazucchi and Laufs, 2014), presumably due to changes in relative mode amplitudes.

Because signals are measured on an array of measurement points (e.g., voxels or discretized regions of interest), we use matrix notation for the covariance matrix and resulting connectivities. But we stress that the discretization must be sufficiently fine to resolve phenomena of interest and reiterate that one must not confuse the structure of the discrete array of measurement points with that of a (nonexistent) discrete macroscale cortical network.

Having calculated the normalized covariance matrix $\overline{C}$ by approximating Equation (3) for τ = 0 by truncating the integral to a long window of duration t_L, we obtain its orthonormal eigenfunctions ${\overline{u}}_j$ by standard matrix methods. For symmetric connectivity we can obtain $\overline{T}$ and $\overline{\Lambda }$ from Equations (29–31). More generally, even in the asymmetric case, spectral factorization methods can be used to calculate the corresponding nonorthonormal ${\overline{v}}_j$ in long-term averaged expressions analogous to Equations (17, 21) if the full correlation function is retained.

We demonstrate the above steps by using the Human Connectome Project HCP 1200 resting fMRI dataset (Glasser et al., 2013). We use postprocessed resting-state fMRI timeseries that have passed through a data processing pipeline to register the measurements onto a standard surface (Glasser et al., 2013) and denoise the data (Glasser et al., 2016) using ICA-FIX. From this release of the data, we remove subjects whose data exhibited incomplete artifact removal as identified in Elam (2019), which yielded a total of 932 subjects. Data are parcellated using the 1000 ROI Schaefer parcellation (Schaefer et al., 2018). To further reduce noise and the effects of incomplete artfact removal, we average C over all subject resting state recordings and all 932 subjects. Averaging over individuals restricts us to analyzing low order modes whose spatial variations are larger than the sizes of differences between brains of individuals and therefore relatively unaffected by this averaging compared to high order modes, but this suffices for the present illustrative purposes. Figure 11 shows that the first four modes obtained from this subject averaged fMRI $\overline{C}$ are consistent with the first four modes obtained by solving the Helmholtz equation on the curved cortical surface in section 2.3.3. These low-order $\overline{C}$ modes are obtained purely from data, thus confirming that at least the lowest-order surface-based modes are good approximations, although noise effects on the activity-based modes increase at higher mode numbers. It was previously shown that the lowest order modes were also very similar to ones based on the structural connectivity (Robinson et al., 2016), a fact that follows ultimately from the close similarity of the white matter distribution to the kernel of the Laplacian operator when written in integral form (Robinson et al., 1997). Indeed, the indexes lm of the Y_lm can be used to label the corresponding modes on the convoluted cortex.

FIGURE 11

Figure 11. First four modes for the left hemisphere. (A) Modes obtained by solving the Helmholtz equation on the cortical surface. (B) Covariance-based modes obtained from the average of $\overline{C}$ over subjects in the HCP1200 dataset; for each subject $\overline{C}$ was calculated for a 14 min recording.

Commonly discussed fMRI RSNs and EEG microstates have typical linear scales of many cm, so their decomposition must be dominated by low-order eigenmodes. As shown in Figure 8, such modes are spatially pinned by the large-scale curvature of the cortex and even strong localized perturbations do not change them significantly. Moreover, the large-scale curvature produces eigenvalue perturbations that are tens of per cent (Robinson et al., 2016; Gabay and Robinson, 2017). Hence, only large-scale perturbations of other parameters that satisfy selection rules analogous to Equations (57, 58) and have an amplitude large enough to cause similar-sized shifts could potentially change the modal structure significantly. However, early stages of sensory processing are transient and localized to specific areas of cortex, so associated changes will leave large-scale modal structure almost unchanged, but can be expected to alter the relative activity levels in these modes, which will also change the FC; such changes can also result from spontaneous background neural activity.

3.2. Instantaneous Activity and Connectivity Kernels

If we use long-term real eigenfunctions ${\overline{u}}_j$ to decompose activity as

$\begin{array}{ll}Q(\text{r},t)={\displaystyle \sum _j}{c}_j(t){\overline{u}}_j(\text{r}),& (61)\end{array}$

$\begin{array}{ll}{c}_j(t)={\displaystyle \int }{\overline{u}}_j(\text{r})Q(\text{r},t)d^2\text{r},& (62)\end{array}$

where the integral is over the cortical surface and the ${\overline{u}}_j$ are constructed to be real. In matrix notation, these equations become

$\begin{array}{ll}Q(t)={\displaystyle \sum _j}{c}_j(t){\overline{u}}_j,& (63)\end{array}$

$\begin{array}{ll}{c}_j(t)={\overline{u}}_j^{T}Q(t),& (64)\end{array}$

where the superscript T denotes the transpose. Note that Equations (62, 64) are evaluated by integrating over space at a fixed time without temporal averaging.

Figure 12A shows an illustrative example of the temporal coefficients c_j(t) of covariance modes obtained from measurements. The first nine left-hemisphere covariance matrix mode coefficients, c_j(t), are shown for a single subject, single resting state fMRI scan. Data are from the Human Connectome Project 1200 Subject Database (Glasser et al., 2013). The vertex level data are de-meaned and normalized by the mean before being parcellated into 1000 ROIs using the Schaefer 17 network parcellation (Schaefer et al., 2018). The covariance matrix is computed from the ROI timeseries using Equation (3) with τ = 0, and its eigenvectors are then calculated using standard matrix methods (e.g., MatLab's eig function) to give the spatial component of the covariance modes. The spatial components of the covariance modes are then projected back onto the ROI timeseries using Equation (64) to give the mode coefficients c_j(t). Figure 12B shows that the low order coefficients dominate the power spectrum of the measurements; relatively few modes are required to capture most of the observed dynamics, in comparison to using an ROI or vertex basis of hundreds or thousands of coefficients, respectively.

FIGURE 12

Figure 12. Temporal dependence of mode coefficients. (A) c_j(t), for the first nine covariance modes from a single resting state fMRI scan in the HCP1200 dataset (Glasser et al., 2013). For each mode j vertical axis tick marks are at −0.1, 0, 0.1. (B) The power of each covariance mode $\langle c_j{(t)}^2\rangle$.

Setting aside the integral over t, the instantaneous kernel of the correlation function in Equation (3) can now be written as

$\begin{array}{ll}\tilde{C}(\text{r},t+\tau ;{\text{r}}^{\prime },t)={\displaystyle \sum _j}{\overline{u}}_j(\text{r}){\overline{u}}_j({\text{r}}^{\prime }){\kappa }_j(\tau ,t),& (65)\end{array}$

$\begin{array}{ll}{\kappa }_j(\tau ,t)=c_j(t+\tau )c_j(t).& (66)\end{array}$

This result decomposes the kernel into a sum of terms that correspond directly to the eigenmodes and whose spatial dependence is evaluated without averaging over multiple temporal measurements. Because the evaluation of the temporal coefficients c_j via Equation (64) involves summation over many spatial points in the discrete approximation, noise is reduced without blurring the signals in time. For example, in fMRI there can be ~ 10⁴ voxels, so noise is reduced by about a factor of 100 relative to the single-voxel value at any fixed t, whereas obtaining the same reduction by using ~ 10⁴ temporal measurements would typically take hours, completely average out any short-term dynamics, and almost certainly involve connectivity changes due to drowsiness (Tagliazucchi and Laufs, 2014). If necessary, it is still possible to smooth or interpolate the c_j(t) over a small number of time points by fitting a moving average, for example, to get an improved estimate at any given time, but we do not use this step in the present work. If robust ${\overline{u}}_j$ are available, the above procedure directly projects out the dominant, large-scale modes of interest in studying resting-state fMRI and EEG activity, and suppresses fine-scale spatial features (of noise and activity) in the process, leaving them in the higher modes. This is far more straightforward and mathematically justified than trying to obtain such patterns directly from correlation matrices of voxel- or ROI-level time series that are then thresholded, clustered, and/or sparsified to reduce complexity, and time-averaged to suppress noise.

We can also define the symmetric kernel

$\begin{array}{ll}{\kappa }_j^{\prime }(\tau ,t)=c_j(t+\tau /2)c_j(t-\tau /2),& (67)\end{array}$

which enables us to define the analogs C(τ, t) and K(τ, t) of C(τ) and K(τ) in Equations (31, 37); covariance corresponds to τ = 0. Equation (67) can be Fourier transformed with respect to τ to obtain K(ω, t), the ω-dependent kernel centered at time t, and thence, via the analogs of Equations (29–33), the ECs T(ω, t, t_F) and Λ(ω, t, t_F), where these Fourier-space quantities necessarily involve integration over a temporal window of width t_F to accomplish the Fourier transform. Hence, there is an unavoidable trade-off between spectral and temporal resolution, dictated by the uncertainty principle, and κ_j(τ, t) will vary with t due to time-varying activity, even for a static brain structure. Similar results can also be defined for wavelet transforms.

It is also possible to average κ_j(τ, t) over a range Δt₀ of t₀; i.e.,

$\begin{array}{ll}{\overline{\kappa }}_j(\tau ,t,\Delta t)={\displaystyle {\int }_{-\infty }^{\infty }}{\kappa }_j(\tau ,t+t^{\prime })w(t^{\prime },\Delta t)d{t}^{\prime },& (68)\end{array}$

where w(t, Δt) is a weight function with unit integral and a characteristic width Δt in time. In the limit that Δt → ∞, ${\overline{\kappa }}_j$ approaches its long-term average and the equations in section 2.2 can be used without modification, while at shorter times, we see from the definition (67) that κ_j is second order in the amplitude c_j and thus varies on the same timescale as the correlation time (i.e., inverse bandwidth) of the activity. Correlations of the covariances measured at different times are sometimes used as a measure of dynamic FC; however, such correlations are fourth order in the amplitudes, so it cannot have a longer correlation time than the usual second-order activity correlation.

Overall, it is worth stressing that all covariance-based connectivity measures (and by analogy, ones based on other quantities such as coherence) will vary over time due to the nonzero bandwidth of the activity from which they are computed, even if the actual brain effective connectivity is static. We return to this point below.

3.3. Direct Determination of T

Having obtained the ${\overline{u}}_j$ from the covariance, for example, one can construct T in an alternate way that only uses data from time intervals of order 1 s, and then use it to obtain $\widehat{\Lambda }$. The fact that T is by definition the response to a delta function stimulus implies that the response T(r, r₀, τ) to a unit delta-function stimulus at r₀ and t = 0 is given by Equation (2) with $N({\text{r}}^{\prime },t)={\delta }^2({\text{r}}^{\prime }-{\text{r}}_0)\delta (t)$, which yields

$\begin{array}{ll}{\theta }_j(\tau )=\frac{1}{{\overline{u}}_j({\text{r}}_0)}{\displaystyle \int }{\overline{u}}_j(\text{r})T(\text{r},{\text{r}}_0,\tau )d^2\text{r}.& (69)\end{array}$

If the eigenmodes are determined via a high resolution method such as fMRI, or by solutions of the Helmholtz equation, Equation (69) then enables the coefficients θ_j(τ) in Equation (36) to be determined by EEG with good temporal resolution once T(r, r₀, τ) has been measured at a number of points. Significantly, this does not require high EEG spatial resolution because, by analogy to Fourier transforms, the lowest M/2 of the θ_j can be determined by sampling at M points, labeled k. One finds (Robinson, 2013; Robinson et al., 2018; Henderson et al., 2021).

$\begin{array}{ll}{a}_k(t)=T({\text{r}}_k,{\text{r}}_0,\tau ),& (70)\end{array}$

$\begin{array}{ll}\approx {\displaystyle \sum _{k=1}^M}{\overline{u}}_j({\text{r}}_k){\overline{u}}_j({\text{r}}_0){\theta }_j(\tau ),& (71)\end{array}$

$\begin{array}{ll}={\displaystyle \sum _{k=1}^M}{b}_{kj}{\theta }_j(\tau ),& (72)\end{array}$

$\begin{array}{ll}{b}_{kj}={\overline{u}}_j({\text{r}}_k){\overline{u}}_j({\text{r}}_0),& (73)\end{array}$

where Equation (70) defines a_k(t) and Equation (36) has been used in obtaining Equation (71). Equation (72) can be inverted to obtain the θ_j(τ), and thence T, so long as the condition number of the matrix of the b_kj is not too large, which would correspond to a measurement or stimulus point lying very near a zero of one of the modes. The simplest solution in such a case is to add or substitute more measurement points, which would also reduce the effects of noise.

This procedure enables fusion of high spatial resolution fMRI data with high temporal resolution ERP time series (see Figure 1A) via a sparse representation of modes that is analogous to the use of amplitude and phase to describe Fourier modes (Robinson, 2013). It does employ long-term eigenmodes, but does not involve any other windowing of the data, and can probe rapid changes in response via standard ERP protocols (Luck and Kappenman, 2012). In principle, by allowing sparse nonuniform sampling, it also solves the problem of the relatively coarse-grained resolution of EEG measurements and their inability to probe activity in sulci (Nunez, 1995).

3.4. Dynamic FC

Functional connectivity is most commonly defined to be the two-point covariance of signals, often normalized to constant variance and evaluated using fMRI averaged over a long window of length t_L. However, Yule showed that windowed covariance can exhibit artifactual dynamics even for a pair of sinusoidal signals of fixed relative phase, due to the interaction between the window length and the wave period (Yule, 1926). This analysis has been extended to two signals generated by Wiener processes, where the covariance has been found to be volatile even in a system of fixed structure (Yule, 1926; Ernst et al., 2017), and has also been revisited by Zalesky and Breakspear (2015), Preti et al. (2017), and Leonardi et al. (2015), and has been further clarified above. The problem is to try to separate the covariance changes that occur due to the finite bandwidth of the activity and the window length from those due to underlying EC changes.

3.4.1. Violin String Analog System

Here we analyze a simple example of the dynamics of FC under spontaneous conditions to clarify the roles of bandwidth, mode beating, and transients, and point out intrinsic shortcomings of analysis methods that seek “switching” between patterns from a precompiled dictionary or cluster-based analysis. To do this we use a simple 1D two-mode example to examine the ways in which mode frequencies and bandwidths can affect FC even when there is no change in the underlying EC. All the results here can be straightforwardly generalized to include more modes, continuous frequency spectra, 2D effects, and different boundary conditions.

Let us first consider a 1D brain akin to the violin string in Figure 2, in the absence of damping and regeneration of activity, but now with periodic boundary conditions that identify the left and right boundaries with one another and allow nonzero displacement at all points, including the two points that were fixed in Figure 2. The lowest mode here is uniform in space and k_j = 2πj/L with j = 0, 1, 2, ….

If all the activity is in the lowest two nonuniform modes, and the transfer function has two pairs of poles at slightly different frequencies we have

$\begin{array}{l}Q(x,t)=u_1(x)\left[sin({\Omega }_{11}t)+sin({\Omega }_{12}t)\right]\\ +u_2(x)\left[sin({\Omega }_{21}t)+sin({\Omega }_{22}t)\right],\\ =2u_1(x)sin({\Omega }_{1}t)cos(\Delta {\omega }_{1}t)& (74)\end{array}$

$\begin{array}{ll}+2u_2(x)sin({\Omega }_{2}t)cos(\Delta {\omega }_{2}t),& (75)\end{array}$

$\begin{array}{ll}{u}_1(x)=sin(2\pi x/L),& (76)\end{array}$

$\begin{array}{ll}{u}_2(x)=cos(2\pi x/L),& (77)\end{array}$

$\begin{array}{ll}{\Omega }_j=({\Omega }_{j1}+{\Omega }_{j2})/2,& (78)\end{array}$

$\begin{array}{ll}\Delta {\omega }_j=({\Omega }_{j1}-{\Omega }_{j2})/2& (79)\end{array}$

in the simplified case that all poles are of equal strength. Typically, nonzero mode frequencies in the corticothalamic system are of order 3–20 Hz (Babaie-Janvier and Robinson, 2018, 2019, 2020; Gabay et al., 2018) and the frequency difference between poles in the same band is of order 1 Hz or less (Gabay and Robinson, 2017; Gabay et al., 2018). Moreover, multiple modes typically have resonances close to each other in a given band (Gabay et al., 2018). In Equation (74) the mode coefficients vary sinusoidally in time, but if a form analogous to Equation (59) is used, we can separate off the sinusoidal part and fully specify the dynamics by the constant coefficients of these functions, which is even more compact. More generally, these coefficents may vary on timescales longer than the period of the oscillations (if they varied on a shorter timescale, clear oscillations would not be seen).

Equation (75) shows that each mode will display fast oscillations at each nonzero Ω_j and slow beating at Δω_j. These have been shown to account for properties of EEG microstates, including effects that lead to preferential detection of the particular patterns that occur at the extremums of the beats (Gabay et al., 2018).

If we write the instantaneous contribution of the above activity to the covariance as c(r, r′, t, τ = 0) before any temporal averaging, then

$\begin{array}{l}c(\text{r},{\text{r}}^{\prime },t,\tau =0)=u_1(x)u_1(x^{\prime }){sin}^2({\Omega }_{1}t){cos}^2(\Delta {\omega }_{1}t)\\ +u_2(x)u_2(x^{\prime }){sin}^2({\Omega }_{2}t){cos}^2(\Delta {\omega }_{2}t)\\ +\left[u_1(x)u_2(x^{\prime })+u_2(x)u_1(x^{\prime })\right]\\ \times sin({\Omega }_{1}t)sin({\Omega }_{2}t)cos(\Delta {\omega }_{1}t)cos(\Delta {\omega }_{2}t),& (80)\end{array}$

which gives a FC that is manifestly time- and frequency-dependent. Again, if averaging is done over the fast oscillations, only the first two terms on the right remain, modulated by slow oscillations that arise from the beats, and giving

$\begin{array}{l}c(\text{r},{\text{r}}^{\prime },t,\tau =0)=\frac{1}{2}{u}_1(x)u_1(x^{\prime }){cos}^2(\Delta {\omega }_{1}t)\\ +\frac{1}{2}{u}_2(x)u_2(x^{\prime }){cos}^2(\Delta {\omega }_{2}t).& (81)\end{array}$

The FC will thus vary on long timescales ~1/|Δω_j|, spending disproportionately large time intervals near one or other extremum of the cos² factors, and switching quasiperiodically between them; however, these extremums are not “metastable” (i.e., in stable local equilibrium but with a small margin of stability relative to other more stable states), contrary to a misleading usage that has crept into some corners of the literature; the term “metastable” should not be used to describe patterns that are merely transient.

Because Δω₁ ≠ Δω₂ in general, the two contributions in Equation (76) will gradually drift in and out of phase, so FC estimates on scales of ~ 1 s will vary even without any change in modes, modal structure, or mean activity in the individual modes. Further averaging over the beat timescales would replace the cos² terms by 1/2. Hence, we have three regimes of dynamics and window-length t_L for corresponding averaging: (i) t_L ≫ 2π/|Δω_j| where long-term average FC is recovered; (ii) 2π/Ω_j ≲ t_L ≲ 2π/|Δω_j| where beats of activity induce FC modulations and interact with the window length to give spurious apparent dynamics; and (iii) t_L ≲ 2π/Ω_j where fast oscillations modulate FC on timescales detectable by EEG methods.

Turning to the spatial properties of Equations (75, 81), we first note that points where u_j(x) and $u_j(x^{\prime })$ have the same sign will be positively correlated, as shown in Figures 13A,B. This means that there will be strong positive correlations near x = L/4, 3L/4 for mode 1, but these two regions will be strongly anticorrelated with one another. Hence, methods that group regions with high positive correlations into RSNs will erroneously split each mode except the lowest into two separate RSNs, as seen in Figures 13C,D, despite it being a single dynamical entity. Seed-based clustering methods that ignore negative correlations will thus tend to return four separate RSNs in place of the two modes in Equation (80). Similar comments apply to methods that cluster spatial points on the basis of similarity of patterns of incoming or outgoing connectivity (Yeo et al., 2011). Many of these methods do find robust patterns of the type they impose, but they break the underlying dynamical links and suffer from thresholding effects, artificial discretization, omission of negative correlations, assumptions of nonoverlapping patterns, and/or other restrictions built into the methods themselves, as discussed in the Introduction.

FIGURE 13

Figure 13. Schematic RSNs of a violin string with periodic boundary conditions. (A) Correlations between points at x and x′ for the mode shown in (B), with positive values shown red and negative blue. (B) One sinusoidal mode. (C) Correlations with the points L/4 (solid) and 3L/4 (dashed). (D) Thresholded correlations for the two cases in (C). (E) RSNs defined by large positive correlations between points.

The temporal variations noted above will cause the peaks and troughs of the covariance to wax and wane on timescales ~ 2π/|Δω_j|. In general, this will cause the two pairs of regions of highly correlated activity to come and go so that 0, 2, or 4 extremums will be detectable above some magnitude threshold, which would normally be set by experimental noise or imposed by fiat to keep the number of patterns small. Methods that rely on matching to a dictionary of predetermined FC patterns will thus appear to show sudden switching between patterns, despite the underlying dynamics being smooth and without jumps. Such apparent transitions are thus a spurious artifact of moment-by-moment matching of continuous evolution to the nearest pattern in a discrete list.

3.4.2. Limitations on Detecting Dynamic FC

The example in the previous subsection can be generalized to cases in which the power spectrum is continuous and many resonances are present at once. In the most relevant case, in which spontaneous activity spectra peak at zero frequency, one simply sets Ω_j = 0 and replaces the bandwidth Δω_j by the total bandwidth of the spectrum.

Figure 14 shows schematic spectra and temporal correlation functions of spontaneous fMRI BOLD signals and EEG or MEG signals, normalized to unity at low frequencies. These have bandwidths and correlation times t_c of order Δf ≈ 0.07 s and t_c ≈ 15 s for BOLD and Δf ≈ 15 Hz and t_c ≈ 0.07 s for EEG or MEG, respectively.

FIGURE 14

Figure 14. Frequencies and timescales relevant to FC dynamics. (A) Schematic EEG (solid) and BOLD (dashed) spectra, with characteristic bandwidths of ~ 15 Hz and ~ 0.1 Hz, respectively. (B) Schematic correlation functions of EEG (red) and BOLD (blue) fluctuations, with characteristic correlation times of 0.1 s and 15 s, respectively.

Two regimes are defined by the correlation functions in Figure 14B. The first regime corresponds to measurement times much less than t_c. In this case activity fluctuations (and hence fluctuations in the FC even for static EC) are essentially static during a given measurement. Here it is possible to obtain a “snapshot” of FC, but this timescale (≪0.2 s) is too short to be useful for EEG or MEG. In the case of fMRI timescales of a few seconds can be probed and such responses can be incorporated into the evoked response, as discussed in section 2.3.1.

The second regime corresponds to measurement times much greater than t_c. Here, fluctuations in FC on timescales of order t_c will tend to average out, thereby exposing modifications in FC that are due to EC changes on longer timescales. For fMRI this restricts attention to timescales of a 100 s or more, as has been argued previously (Leonardi et al., 2015; Zalesky and Breakspear, 2015). EEG and MEG methods can reach this regime on the scale of a few seconds or longer. This has the potential to enable EC changes to be monitored on timescales much shorter than those of alertness dynamics, and potentially to be probed by the evoked-response approach in section 3.3.

4. Summary and Discussion

We have used eigenfunction-based spectral methods to analyze dynamical brain connectivity in the usual linear regime, drawing together and generalizing many results from spectral analysis, propagator theory, and neural field theory in the process. The results clarify the limitations on a range of ad hoc methods for FC analysis and point the way to avoiding some of their pitfalls, although it is not possible to analyze individually the multitude that have been proposed.

Central to this work is the recognition that the brain is a physical system, understanding of which is facilitated by taking its physical nature into account. This is reflected in our use of system eigenfunctions as the physical building blocks of its dynamics and connectivity, rather than employing statistically constructed “networks.” We also avoid treating the set of links between artificially discretized brain regions as being an actual macroscale brain network—no such network exists. We do acknowledge that some such phenomenological approaches yield measures that are useful for classification and estabishment of case-control distinctions, for example, but they tend to obscure the nature of the underlying dynamics, which is the focus of the present work.

The main findings are:

(i) Links between activity, EC, and FC have been clarified, including short-term modulation of EC as part of the evoked response to a stimulus (e.g., via gain dynamics) being able to be absorbed into the definition of the transfer function (Babaie-Janvier and Robinson, 2020). The EC-FC connection is fundamental because directly observable structural connectivity (SC) does not capture short-range connectivity, especially of inhibitory neurons, and also says nothing about which neurons are active in a given situation. The SC-EC connection is thus a separate issue.

(ii) Eigenmodes of a violin-string analog system were used as a foil to elucidate the concepts, spatial brain eigenfunctions being analogous to the spatial eigenfunctions that correspond to individual notes from a musical instrument. This analogy has been used to clarify how every eigenmode spans the whole system and every point participates in each eigenmode. This invalidates classes of methods based on thresholding and/or clustering that seek to find discrete patterns (or “networks”) that have sharp boundaries and are mutually exclusive.

(iii) Eigenmode analysis emphasizes how smooth evolution of activity gives rise to time-dependent FC. Beating of modes at different frequencies will give rise to changes in activity and FC patterns without requiring nonlinear mechanisms (Gabay and Robinson, 2017). Resulting dominance of different modes at different times can yield substantial FC changes even when the EC is static. Methods that match instantaneous FC to the nearest entry in a finite dictionary of patterns will produce sudden, but artifactual, “switching” as a result.

(iv) EC and FC are dominated by a moderate number of low-order modes, which explains the robustness of a relatively small number of patterns seen in EEG and fMRI experiments—microstates and resting-state networks, respectively. The dominance of just a few modes also explains the observed similarity of FC at timescales ranging from seconds to tens of minutes.

(v) There is not a one-to-one mapping between eigenmodes and patterns such as RSNs or microstates. Generally, phenomenological patterns can be decomposed into superpositions of multiple eigenmodes, but the common steps of thresholding and imposition of mutual exclusivity usually lead to the results being noncompact and subject to serious artifact. If such steps are avoided, few-component expansions should be possible, with observed patterns being akin to musical chords, with eigenmodes corresponding to notes.

(vi) Eigenfunctions have been shown to be robust to perturbations except where these are large and of similar spatial scale to variations in the eigenfunction amplitude, as governed by selection rules analogous to those for interactions in quantum physics. Eigenmode robustness explains the general similarity of FC patterns across arousal states. Dynamic changes in FC are thus likely to be dominated by activity changes, rather than changes in EC; for example, deep sleep has enhanced low-frequency activity, which is correlated with lower-order modes and larger correlation lengths and may explain FC differences between arousal states (Nunez, 1995; Robinson, 2003; Chang et al., 2016).

(vii) It has been shown that robust low-order eigenmodes can be estimated from fMRI covariance data, for example, and that these are very similar to modes found by solving the Helmholtz equation on the cortical surface, which further supports their validity. This finding complements earlier work that showed similarities between surface-based modes and ones obtained from structural connectivity (Robinson et al., 2016). Higher-order covariance-based modes are increasingly affected by noise, at least for the data pipeline used in the HCP1200 data set. Denoising is an active area of research that promises to deliver better results in the future, but is beyond the scope of the present aper.

(viii) Projection of activity onto eigenmodes enables direct estimation of their coefficients as functions of time. For activity evoked by an impulse stimulus, this permits calculation of the transfer function, and thus the EC and implied FC, on timescales as short as a second or so using EEG or MEG. This method also has the advantage that multiple presentations under task or background conditions can be used to average out the effects of uncontrolled background activity, as is usual in evoked response experiments.

(ix) Calculation of eigenmode coefficients provides a compact few-mode representation of activity, EC, and the kernels of FC without averaging over time. This reduces noise because short-scale noise predominantly affects high-order spatial modes, removes artifacts due to windowing, and improves time resolution. Tracking mode coefficients vs. time is thus likely to be superior to covariance-based dynamic FC for many purposes.

(x) Because activity changes contribute to FC dynamics, changing FC is a mandatory consequence of the finite bandwidth of the signals even when EC is fixed, and the window length used to evaluate mean FC compounds this with artifact (Yule, 1926). These effects are not always separable, but the temporal regimes required to probe them have been elucidated, with the correlation time t_c of background activity being the key parameter. For observations on timescales much less than t_c, activity fluctuations can be viewed as static, whereas they tend to average out for times much greater than t_c, thereby potentially exposing EC changes that are not direct responses to activity.

(xi) As noted above, eigenmodes form a discrete set, with activity dominated by the lowest-order members. Moreover, the frequency response of each mode can be closely approximated by a discrete sum of resonant responses. The latter sum can be truncated to a few terms, each of which represents the dynamics of a PID filter that links the resulting activity to control-systems functions such as prediction and gain-adjustment, the latter implementing a form of attention (Babaie-Janvier and Robinson, 2018, 2019, 2020). These facts imply that eigenmodes form an actual discrete set of entities for the transmission of brain activity, as opposed to “networks” constructed artificially on the basis of arbitrary sensor locations or brain parcellations. The difference here is that the “nodes” are spatially extended eigenmodes; however, these are local in the spectral domain, where they form a discrete lattice labeled by spatial eigenvalues and resonant frequencies.

(xii) Spatial modes plausibly enable spatial communications channels with activity broadcasted into modes, especially near their peaks (Robinson et al., 2018), and available to be read out near other peaks. Significantly, as the number of active modes increases, a steep increase in the number of such channels occurs. Likewise, mode resonances enhance temporal coherence and may thus facilitate a version of the communication through coherence proposal, where in-phase activity can lead to nonlinear enhancement of local firing rates (Fries, 2005).

(xiii) Just as in fields such as nonlinear optics and nonlinear plasma theory, linear modes provide the appropriate starting point to analyze such dynamics systematically via perturbation techniques such as mode-coupling expansions (Butcher and Cotter, 1990; Melrose and McPhedran, 1991).

Overall, the above findings reinforce the advantages of using physically based approaches that work with physical objects and respect the underlying dynamics, and analyzing them via standard mathematical methods whose properties and limitations are thoroughly understood. Such approaches avoid working with unphysical objects, failure to respect dimensionality and units, category errors such as attributing discrete properties to continuous systems or quantities, and forced attribution to the brain of characteristics that are actually imposed by the analysis methods used (e.g., discreteness, graph-theoretic properties such as motifs and degree, sharp boundaries, tree-like hierarchies, or modularity, although some of these can be redefined to be valid in the limit of increasingly fine discretization Robinson, 2019). Seductive, but largely unjustified, steps such as sparsification, thresholding, and clustering are also avoided, along with the use of idiosyncratic analyses whose properties and relationships to other approaches are poorly understood. Thus, in keeping with Occam's Razor, new mathematical methods and constructs should be eschewed unless established methods are inadequate.

Systematic methods enable many problems to be attacked that could not otherwise be addressed, and the results to be unified into an overarching framework that spans activity and connectivity. However, the flip-side of this is that one must learn and acquire facility with a range of well-established mathematical and physical techniques—statistical methods alone are not enough to understand dynamics and connectivity beyond phenomenological classification.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: Human connectome project HCP 1200 Subject Resting fMRI dataset, https://www.humanconnectome.org/study/hcp-young-adult.

Author Contributions

PR conceived the project, carried out the analytic work, and wrote most of the manuscript. JH, KA, and NG contributed to the analysis. JH, KA, NG, TB-J, and XG contributed numerical work embodied in quantitative figures. All authors reviewed the MS and contributed to its final text.

Funding

Australian Research Council, Laureate Fellowship Grant, FL1401000225 Australian Research Council, Center of Excellence Grant, CE140100007.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank A. Breviario and D. Naoumenko for drafting Figures 1, 5.

References

Abeysuriya, R. G., Rennie, C. J., and Robinson, P. A. (2015). Physiologically based arousal state estimation and dynamics. J. Neurosci. Meth. 253, 55–69. doi: 10.1016/j.jneumeth.2015.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Adachi, Y., Osado, T., Sporns, O., Watanabe, T., Matsui, T., Miyamoto, K., et al. (2012). Functional connectivity between anatomically unconnected areas is shaped by collective network-level effects in the macaque cortex. Cereb. Cortex 22, 1586–1592. doi: 10.1093/cercor/bhr234

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, E. A., Damaraju, E., Plis, S. M., Erhardt, E. B., Eichele, T., and Calhoun, V. D. (2014). Tracking whole-brain connectivity dynamics in the resting state. Cereb. Cortex 24, 663–676. doi: 10.1093/cercor/bhs352

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, M., Fu, G.-S, Phylpo, R., and Adali, T. (2014). Independent vector analysis: identification conditions and performance bounds. IEEE Trans. Signal Proc. 62:4399. doi: 10.1109/TSP.2014.2333554

CrossRef Full Text | Google Scholar

Arfanakis, K., Cordes, D., Haughton, V. M., Moritz, C. H., Quigley, M. A., and Meyerand, M. E. (2000). Magn. Res. Imag. 18, 921–930. doi: 10.1016/S0730-725X(00)00190-9

CrossRef Full Text | Google Scholar

Atasoy, S., Deco, G., Kringelbach, M. L., and Pearson, J. (2018). Harmonic brain modes: a unifying framework for linking space and time in brain dynamics. Neuroscientist 24, 277–293. doi: 10.1177/1073858417728032

PubMed Abstract | CrossRef Full Text | Google Scholar

Atasoy, S., Donnelly, I., and Pearson, J. (2016). Human brain networks function in connectome-specific harmonic waves. Nat. Commun. 7, 10340. doi: 10.1038/ncomms10340

PubMed Abstract | CrossRef Full Text | Google Scholar

Babaie-Janvier, T., and Robinson, P. A. (2018). Neural field theory of corticothalamic prediction with control systems analysis. Front. Hum. Neurosci. 12:334. doi: 10.3389/fnhum.2018.00334

PubMed Abstract | CrossRef Full Text | Google Scholar

Babaie-Janvier T. and Robinson, P. A. (2019). Neural field theory of corticothalamic attention with control systems analysis. Front. Neurosci. 13:1240. doi: 10.3389/fnins.2019.01240

CrossRef Full Text | Google Scholar

Babaie-Janvier T. and Robinson, P. A. (2020). Neural field theory of evoked response potentials with attentional gain dynamics. Front. Hum. Neurosci. 14:293. doi: 10.3389/fnhum.2020.00293

CrossRef Full Text | Google Scholar

Bassett, D. S., and Sporns, O. (2017). Network neuroscience. Nat. Neurosci. 20, 353–364. doi: 10.1038/nn.4502

CrossRef Full Text | Google Scholar

Bassett, D. S., Wymbs, N. F., Porter, M. A., Mucha, P. J., Carlson, J. M., and Grafton, S. T. (2011). Dynamic reconfiguration of human brain networks during learning. Proc. Nat. Acad. Sci. U.S.A. 108, 7641–7646. doi: 10.1073/pnas.1018985108

PubMed Abstract | CrossRef Full Text | Google Scholar

Beckmann, C. F., DeLuca, M., Devlin, J. T., and Smith, S. M. (2005). Investigations into resting-state connectivity using independent component analysis. Phil. Trans. Roy. Soc. B. 360, 1001–1013. doi: 10.1098/rstb.2005.1634

PubMed Abstract | CrossRef Full Text | Google Scholar

Biswal, B., Zerrin Yetkin, F., Haughton, V. M., and Hyde, J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Res. Med. 34, 537–541. doi: 10.1002/mrm.1910340409

PubMed Abstract | CrossRef Full Text | Google Scholar

Braitenberg, V., and Schüz, A. (1998). Cortex: Statistics and Geometry of Neuronal Connectivity, 2nd Edn. Berlin: Springer.

Google Scholar

Britz, J., Van De Ville, D., and Michel, C. M. (2010). BOLD correlates of EEG topography reveal rapid resting-state network dynamics. NeuroImage 52, 1162–1170. doi: 10.1016/j.neuroimage.2010.02.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Bullmore, E. T., and Sporns, O. (2009). Complex brain networks: Graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186–198. doi: 10.1038/nrn2575

PubMed Abstract | CrossRef Full Text | Google Scholar

Butcher, P. N., and Cotter, D. (1990). The Elements of Nonlinear Optics. Cambridge: Cambridge Univ. Press.

Google Scholar

Cabral, J., Luckhoo, H., Woolrich, M., Joensson, M., Mohseni, H., Baker, A., et al. (2014). Exploring mechanisms of spontaneous functional connectivity in MEG: How delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations. NeuroImage 90, 423–435. doi: 10.1016/j.neuroimage.2013.11.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabral, J., Vidaurre, D., Marques, P., Magalhães, R., Moreira, P. S., Soares, J. M., et al. (2017). Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest. Sci. Rep. 7:5135. doi: 10.1038/s41598-017-05425-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Calhoun, V. D., Miller, R., Pearlson, G., and Adali, T. (2014). The chronnectome: time-varying connectivity networks as the next frontier in fMRI data discovery. Neuron 84, 262–274. doi: 10.1016/j.neuron.2014.10.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, C, and Glover, G. H. (2010). Time–frequency dynamics of resting–state brain connectivity measured with fMRI. NeuroImage 50, 81–98. doi: 10.1016/j.neuroimage.2009.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, C., Leopold, D. A., Schölvinck, M. L., Mandelkow, H., Picchioni, D., Liu, X., et al. (2016). Tracking brain arousal fluctuations with fMRI. Proc. Nat. Acad. Sci. USA. 113, 4518-4523. doi: 10.1073/pnas.1520613113

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, G, Chen, G., Xie, C., and Li, S.-J. (2011). Negative functional connectivity and its dependence on the shortest path length of positive network in the resting-state human brain. Brain Conn. 1, 195–206. doi: 10.1089/brain.2011.0025

PubMed Abstract | CrossRef Full Text | Google Scholar

Coombes, S., Beim Graben, P., Potthast, R., and Wright, J., (eds) (2014). Neural Fields: Theory and Applications. Berlin: Springer.

Google Scholar

Courant, R., and Hilbert, D. (2004). Methods of Mathematical Physics. Weinheim: Wiley-VCH.

Google Scholar

Damoiseaux J. S Rombouts S. A. R. B Barkhof F Scheltens P Stam C. J Smith S. M. . (2006). Consistent resting-state networks across healthy subjects. Proc. Nat. Acad. Sci. U.S.A. 103, 13848–13853. doi: 10.1073/pnas.0601417103

CrossRef Full Text | Google Scholar

Deco, G., and Jirsa, V. K. (2012). Ongoing cortical activity at rest: Criticality, multistability, and ghost attractors. J. Neurosci. 32, 3366–3375. doi: 10.1523/JNEUROSCI.2523-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Jirsa, V. K., and McIntosh, A. R. (2011). Emerging concepts for the dynamical organization of resting-state activity in the brain. Nat. Rev. Neurosci. 12, 43–56. doi: 10.1038/nrn2961

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Jirsa, V. K., Robinson, P. A., Breakspear, M., and Friston, K. (2008). The dynamic brain: From spiking neurons to neural masses and cortical fields. PLoS Comput. Biol. 4:e1000092. doi: 10.1371/journal.pcbi.1000092

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Kringelbach, M. L., Jirsa, V. K., and Ritter, P. (2017). The dynamics of resting fluctuations in the brain: Metastability and its dynamical cortical core. Sci. Rep. 7:3095. doi: 10.1038/s41598-017-03073-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Drysdale, P. M., Huber, J. P., Robinson, P. A., and Aquino, K. M. (2010). Spatiotemporal BOLD dynamics from a poroelastic hemodynamic model. J. Theor. Biol. 265, 524–534. doi: 10.1016/j.jtbi.2010.05.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunster, T. M. (2010). “Legendre and Related Functions,” in NIST Handbook of Mathematical Functions, ed. F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark (Cambridge: Cambridge University Press).

Google Scholar

Eklund, A., Nichols, T. E., and Knutsson, H. (2016). Cluster failure: Why fMRI inferences for spatial extent have inflated false-positive rates. Proc. Nat. Acad. Sci. U.S.A. 113, 7900–7905. doi: 10.1073/pnas.1602413113

PubMed Abstract | CrossRef Full Text | Google Scholar

Elam, J. (2019). HCP Subjects with Identified Quality Control Issues (QC_Issue measure codes explained). Available online at: https://wiki.humanconnectome.org/pages/viewpage.action?pageId=88901591. (Accessed January 14, 2021).

Google Scholar

Ernst, P. A., Shepp, L. A., and Wyner, A. J. (2017). Yule's “nonsense correlation” solved! Ann. Stat. 45, 1789–1809. doi: 10.1214/16-AOS1509

CrossRef Full Text | Google Scholar

Fischl, B., Sereno, M. I., and Dale, A. M. (1999). Cortical surface-based analysis: I. Segmentation and surface reconstruction. NeuroImage 9, 179–194. doi: 10.1006/nimg.1998.0396

PubMed Abstract | CrossRef Full Text | Google Scholar

Fornito, A., Zalesky, A., and Breakspear, M. (2013). Graph analysis of the human connectome: Promise, progress, and pitfalls. NeuroImage 80, 426–444. doi: 10.1016/j.neuroimage.2013.04.087

PubMed Abstract | CrossRef Full Text | Google Scholar

Fornito, A., Zalesky, A., and Bullmore, E. T. (2016). Fundamentals of Brain Network Analysis. Amsterdam: Elsevier.

Google Scholar

Fourier, J.-B. J. (1822). Theorie analytique de la chaleur. Paris: Firmin Didot.

Google Scholar

Fox, M. D, and Raichle, M. E. (2007). Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat. Rev. Neurosci. 8, 700–711. doi: 10.1038/nrn2201

PubMed Abstract | CrossRef Full Text | Google Scholar

Fries, P. (2005). A mechanism for cognitive dynamics: neuronal communication through neuronal coherence. Trends Cogn. Neurosci. 9, 474–480. doi: 10.1016/j.tics.2005.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J. (1994). Functional and effective connectivity in neuroimaging: a synthesis. Hum. Brain Mapp. 2, 56–78. doi: 10.1002/hbm.460020107

CrossRef Full Text | Google Scholar

Friston, K. J. (2011). Functional and effective connectivity: a review. Brain Conn. 1, 13–36. doi: 10.1089/brain.2011.0008

CrossRef Full Text | Google Scholar

Friston, K. J., Kahan, J., Razi, A., Stephan, K. E., and Sporns, O. (2014). On nodes and modes in resting state fMRI. NeuroImage 99, 533–547. doi: 10.1016/j.neuroimage.2014.05.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Gabay, N. C., Babaie-Janvier, T., and Robinson, P. A. (2018). Dynamics of cortical activity eigenmodes including standing, traveling, and rotating waves. Phys. Rev. E. 98:042413. doi: 10.1103/PhysRevE.98.042413

CrossRef Full Text | Google Scholar

Gabay, N. C., and Robinson, P. A. (2017). Cortical geometry as a determinant of brain activity eigenmodes: neural field analysis. Phys. Rev. E. 96:032413. doi: 10.1103/PhysRevE.96.032413

PubMed Abstract | CrossRef Full Text | Google Scholar

Galán, R. F. (2008). On How Network Architecture Determines the Dominant Patterns of Spontaneous Neural Activity. PLoS ONE 3:e2148. doi: 10.1371/journal.pone.0002148

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, X., and Robinson, P. A. (2020). Importance of self-connections for brain connectivity and spectral connectomics. Biol. Cybern. 114, 643–651. doi: 10.1007/s00422-020-00847-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghosh, A., Rho, Y., McIntosh, A. R., Kötter, R., and Jirsa, V. K. (2008). Noise during Rest Enables the Exploration of the Brain's Dynamic Repertoire. PLoS Comput. Biol. 4:e1000196. doi: 10.1371/journal.pcbi.1000196

PubMed Abstract | CrossRef Full Text | Google Scholar

Glasser, M. F., Coalson, T. S., Robinson, E. C., Hacker, C. D., Harwell, J., Yacoub, E., et al. (2016). A multi-modal parcellation of human cerebral cortex. Nature 536, 171–178. doi: 10.1038/nature18933

PubMed Abstract | CrossRef Full Text | Google Scholar

Glasser, M. F., Sotiropoulos, S. N., Wilson, J. A., Coalson, T. S., Fischl, B., Andersson, J. L., et al. (2013). The minimal preprocessing pipelines for the human connectome project. NeuroImage 80, 105–124. doi: 10.1016/j.neuroimage.2013.04.127

PubMed Abstract | CrossRef Full Text | Google Scholar

Gohel, S. R., and Biswal, B. B. (2015). Functional integration between brain regions at rest occurs in multiple-frequency bands. Brain Conn. 5, 23–34. doi: 10.1089/brain.2013.0210

PubMed Abstract | CrossRef Full Text | Google Scholar

Greicius, M. D., Supekar, K., Menon, V., and Dougherty, R. F. (2009). Resting-state functional connectivity reflects structural connectivity in the default mode network. Cereb. Cortex 19, 72–78. doi: 10.1093/cercor/bhn059

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, F. C. A., Battaglia, D., Spiegler, A., Deco, G., and Jirsa, V. K. (2015). Functional connectivity dynamics: Modeling the switching behavior of the resting state. NeuroImage 105, 525–535. doi: 10.1016/j.neuroimage.2014.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Henderson, J. A., Dhamala, M., and Robinson, P. A. (2021). Brain dynamics and structure-function relationships via spectral factorization and the transfer function. NeuroImage 235:117989. doi: 10.1016/j.neuroimage.2021.117989

PubMed Abstract | CrossRef Full Text | Google Scholar

Henderson, J. A., and Robinson, P. A. (2011). Geometric effects on complex network structure in the cortex. Phys. Rev. Lett. 107:018102. doi: 10.1103/PhysRevLett.107.018102

PubMed Abstract | CrossRef Full Text | Google Scholar

Hindriks, R., Adhikari, M. H., Murayama, Y., Ganzetti, M., Mantini, D., Logothetis, N. K., et al. (2016). Can sliding-window correlations reveal dynamic functional connectivity in resting-state fMRI? NeuroImage 127, 242–256. doi: 10.1016/j.neuroimage.2015.11.055

PubMed Abstract | CrossRef Full Text | Google Scholar

Hipp, J. F., Hawellek, D. J., Corbetta, M., Siegel, M., and Engel, A. K. (2012). Large-scale cortical correlation structure of spontaneous oscillatory activity. Nat. Neurosci. 15, 884–890. doi: 10.1038/nn.3101

PubMed Abstract | CrossRef Full Text | Google Scholar

Honey, C. J., Kötter, R., Breakspear, M., and Sporns, O. (2007). Network structure of cerebral cortex shapes functional connectivity on multiple time scales. Proc. Nat. Acad. Sci. U.S.A. 104, 10240–10245. doi: 10.1073/pnas.0701519104

PubMed Abstract | CrossRef Full Text | Google Scholar

Honey, C. J., Sporns, O., Cammoun, L., Gigandet, X., Thiran, J. P., Meuli, R., et al. (2009). Predicting human resting-state functional connectivity from structural connectivity. Proc. Nat. Acad. Sci. U.S.A. 106, 2035–2040. doi: 10.1073/pnas.0811168106

PubMed Abstract | CrossRef Full Text | Google Scholar

Honey, C. J., Thivierge, J.-P., and Sporns, O. (2010). Can structure predict function in the human brain? NeuroImage 52, 766–776. doi: 10.1016/j.neuroimage.2010.01.071

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunyadi, B., Woolrich, M. W., Quinn, A. J., Vidaurre, D., and De Vos, M. (2019). A dynamic system of brain networks revealed by fast transient EEG fluctuations and their fMRI correlates. NeuroImage 185, 72–82. doi: 10.1016/j.neuroimage.2018.09.082

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutchison, R. M., Womelsdorf, T., Allen, E. A., Bandettini, P. A., Calhoun, V. D., Corbetta, M., et al. (2013a). Dynamic functional connectivity: promise, issues, and interpretations. NeuroImage 80, 360–378. doi: 10.1016/j.neuroimage.2013.05.079

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutchison, R. M., Womelsdorf, T., Gati, J. S., Everling, S., and Menon, R. S. (2013b). Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques. Hum. Brain Mapp. 34, 2154–2177. doi: 10.1002/hbm.22058

PubMed Abstract | CrossRef Full Text | Google Scholar

Jezzard, P., Matthews, P. M., and Smith, S. M., (eds). (2001). Functional MRI: An Introduction to Methods. Oxford: Oxford University Press.

Google Scholar

Jirsa, V. K., and Haken, H. (1996). Field Theory of Electromagnetic Brain Activity. Phys. Rev. Lett. 77:960–963. doi: 10.1103/PhysRevLett.77.960

PubMed Abstract | CrossRef Full Text | Google Scholar

Jirsa V. K. Jantzen K. J. Fuchs A. Kelso J. A. S. (2001) “Neural field dynamics on the folded three-dimensional cortical sheet its forward EEG MEG,” in IPMI 2001: Information Processing in Medical Imaging, eds M. F. Insana R. M. Leahy (Berlin: Springer), 286–299. doi: 10.1007/3-540-45729-1_30

CrossRef Full Text | Google Scholar

Koch, C. (1999). Biophysics of Computation. New York, NY: Oxford.

Google Scholar

Kopell, N. J., Gritton, H. J., Whittington, M. A., and Kramer, M. A. (2014). Beyond the connectome: the dynome. Neuron 83, 1319–1328. doi: 10.1016/j.neuron.2014.08.016

CrossRef Full Text | Google Scholar

Leonardi, N., Richiardi, J., Gschwind, M., Simioni, S., Annoni, J.-M., Schluep, M., et al. (2013). Principal components of functional connectivity: A new approach to study dynamic brain connectivity during rest. NeuroImage 83, 937–950. doi: 10.1016/j.neuroimage.2013.07.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Leonardi, N., Shirer, W. R., Greicius, M. D., and Van De Ville, D. (2014). Disentangling dynamic networks: Separated and joint expressions of functional connectivity patterns in time. Hum. Brain Mapp. 35, 5984–5995. doi: 10.1002/hbm.22599

PubMed Abstract | CrossRef Full Text | Google Scholar

Leonardi, N., and Van De Ville, D. (2015). On spurious and real fluctuations of dynamic functional connectivity during rest. NeuroImage 104, 430–436. doi: 10.1016/j.neuroimage.2014.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Luck, S. J., and Kappenman, E. S. (Eds) (2012). The Oxford Handbook of Event-Related Potential Components. Oxford: Oxford University Press.

Google Scholar

MacLaurin, J. N., and Robinson, P. A. (2019). Determination of effective brain connectivity from activity correlations. Phys. Rev. E 99:042404. doi: 10.1103/PhysRevE.99.042404

PubMed Abstract | CrossRef Full Text | Google Scholar

Mantini, D., Perrucci, M. G., Del Gratta, C., Romani, G. L., and Corbetta, M. (2007). Electrophysiological signatures of resting state networks in the human brain. Proc. Nat. Acad. Sci. U.S.A. 104, 13170–13175. doi: 10.1073/pnas.0700668104

PubMed Abstract | CrossRef Full Text | Google Scholar

Marzetti, L., Della Penna, S., Snyder, A. Z., Pizzella, V., Nolte, G., de Pasquale, F., et al. (2013). Frequency specific interaction of MEG resting state activity within and across brain networks as revealed by the multivariate interaction measure. NeuroImage 79, 172–183. doi: 10.1016/j.neuroimage.2013.04.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Maximon, L. C. (2010). “3j, 6j, 9j symbols,” in NIST Handbook of Mathematical Functions, eds F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark (Cambridge: Cambridge University Press), 757.

Google Scholar

Mehta-Pandejee, G., Robinson, P. A., Henderson, J. A., Aquino, K. M., and Sarkar, S. (2017). Inference of direct and multistep effective connectivities from functional connectivity of the brain and of relationships to cortical geometry. J. Neurosci. Meth. 283, 42–54. doi: 10.1016/j.jneumeth.2017.03.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Melrose, D. B., and McPhedran, R. C. (1991). Electromagnetic Processes in Dispersive Media. Cambridge: Cambridge University Press.

Google Scholar

Messe, A., Hütt, M.-T., König, P., and Hilgetag, C. T. (2015). A closer look at the apparent correlation of structural and functional connectivity in excitable neural networks. Sci. Rep. 5:7870. doi: 10.1038/srep07870

PubMed Abstract | CrossRef Full Text | Google Scholar

Messe, A., Rudrauf, D., Benali, H., and Marrelec, G. (2014). Relating structure and function in the human brain: relative contributions of anatomy, stationary dynamics, and non-stationarities. PLoS Comput. Biol. 10:e1003530. doi: 10.1371/journal.pcbi.1003530

PubMed Abstract | CrossRef Full Text | Google Scholar

Michel, C. M., and Koenig, T. (2018). EEG microstates as a tool for studying the temporal dynamics of whole-brain neuronal networks: A review. NeuroImage 180:577–593. doi: 10.1016/j.neuroimage.2017.11.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitra, A., Snyder, A. Z., Hacker, C. D., and Raichle, M. E. (2014). Lag structure in resting-state fMRI. J. Neurophysiol. 111:2374–2391. doi: 10.1152/jn.00804.2013

CrossRef Full Text | Google Scholar

Mukta, K. N., Gao, X., and Robinson, P. A. (2019). Neural field theory of evoked response potentials in a spherical brain geometry. Phys. Rev. E 99:062304. doi: 10.1103/PhysRevE.99.062304

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukta, K. N., MacLaurin, J. N., and Robinson, P. A. (2017). Theory of corticothalamic brain activity in a spherical geometry: Spectra, coherence, and correlation. Phys. Rev. E 96:052410. doi: 10.1103/PhysRevE.96.052410

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukta, K. N., Robinson, P. A., Pagès, J. C., Gabay, N. C., and Gao, X. (2020). Evoked response activity eigenmode analysis in a convoluted cortex via neural field theory. Phys. Rev. E 102:062303.

PubMed Abstract | Google Scholar

Murphy, K., and Fox, M. D. (2017). Towards a consensus regarding global signal regression for resting state functional connectivity MRI. NeuroImage 154, 169–173. doi: 10.1016/j.neuroimage.2016.11.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Musso, F., Brinkmeyer, J, Mobascher, A, Warbrick, T, and Winterer, G. (2010). Spontaneous brain activity and EEG microstates. A novel EEG/fMRI analysis approach to explore resting-state networks. NeuroImage 52, 1149–1161. doi: 10.1016/j.neuroimage.2010.01.093

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunez, P. L. (1995). Neocortical Dynamics and Human EEG Rhythms. New York, NY: Oxford University Press.

Google Scholar

Nunez, P. L., Wingeier, B. M., and Silberstein, R. B. (2001). Spatial-temporal structures of human alpha rhythms: theory, microcurrent sources, multiscale measurements, and global binding of local networks. Hum. Brain Mapp. 13:125–164. doi: 10.1002/hbm.1030

PubMed Abstract | CrossRef Full Text | Google Scholar

Ogata, K., and Yang, Y. (1970). Modern Control Engineering. Englewood Cliffs, NJ: Prentice-Hall.

Google Scholar

Pang J. C. and Robinson, P. A. (2019). Power spectrum of resting-state blood-oxygen-level-dependent signal. Phys. Rev. E 100:022418. doi: 10.1103/PhysRevE.100.022418

CrossRef Full Text | Google Scholar

Peltier, S. J., Kerssens, C., Hamann, S. B., Sebel, P. S., Byas-Smith, M., and Hu, X. (2005). Functional connectivity changes with concentration of sevoflurane anesthesia. NeuroRep. 16:285–288. doi: 10.1097/00001756-200502280-00017

PubMed Abstract | CrossRef Full Text | Google Scholar

Pernice, V., Staude, B., Cardanobile, S., and Rotter, S. (2011). How structure determines correlations in neuronal networks. PLoS Comput. Biol. 7:e1002059. doi: 10.1371/journal.pcbi.1002059

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinotsis, D. A., Hansen, E., Friston, K. J., and Jirsa, V. K. (2013). Anatomical connectivity and the resting state activity of large cortical networks. NeuroImage 65:127–138. doi: 10.1016/j.neuroimage.2012.10.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Preti, M. G., Bolton, T. A. W., and Van De Ville, D. (2017). The dynamic functional connectome: State-of-the-art and perspectives. NeuroImage 160, 41–54. doi: 10.1016/j.neuroimage.2016.12.061

PubMed Abstract | CrossRef Full Text | Google Scholar

Preti, M. G., and Van De Ville, D. (2019). Decoupling of brain function from structure reveals regional behavioral specialization in humans, Nature Commun., 10:41–54. doi: 10.1038/s41467-019-12765-7

CrossRef Full Text | Google Scholar

Raichle, M. E. (2011). The restless brain. Brain Conn. 1, 3–12. doi: 10.1089/brain.2011.0019

CrossRef Full Text | Google Scholar

Raichle, M. E., MacLeod, A. M., Snyder, A. Z., Powers, W. J., Gusnard, D. A., and Shulman, G. L. (2001). A default mode of brain function. Proc. Nat. Acad. Sci. U.S.A. 98:676. doi: 10.1073/pnas.98.2.676

CrossRef Full Text | Google Scholar

Rennie, C. J., Robinson, P. A., and Wright, J. J. (1999). Effects of local feedback on dispersion of electrical waves in the cerebral cortex. Phys. Rev. E 59:3320–3329. doi: 10.1103/PhysRevE.59.3320

CrossRef Full Text | Google Scholar

Rennie, C. J., Robinson, P. A., and Wright, J. J. (2002). Unified neurophysical model of EEG spectra and evoked potentials. Biol. Cybern. 86:457–471.

PubMed Abstract | Google Scholar

Roberts, J. A., Gollo, L. L., Abeysuriya, R. G., Roberts, G., Mitchell, P. B., Woolrich, M. W., et al. (2019). Metastable brain waves Nat. Commun. 10:1056. doi: 10.1038/s41467-019-08999-0

CrossRef Full Text | Google Scholar

Robinson, P. A. (2003). Neurophysical theory of coherence and correlations of electroencephalographic and electrocorticographic signals. J. Theor. Biol. 222:163–175. doi: 10.1016/S0022-5193(03)00023-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A. (2012). Interrelating anatomical, effective, and functional brain connectivity using propagators and neural field theory. Phys. Rev. E 85:011912. doi: 10.1103/PhysRevE.85.011912

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A. (2013). Discrete-network versus modal representations of brain activity: Why a sparse regions-of-interest approach can work for analysis of continuous dynamics. Phys. Rev. E 88:054702. doi: 10.1103/PhysRevE.88.054702

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A. (2019). Physical Brain Connectomics. Phys. Rev. E 99:012421. doi: 10.1103/PhysRevE.99.012421

CrossRef Full Text | Google Scholar

Robinson, P. A., Drysdale, P. M., Van der Merwe, H., Kyriakou, E., Rigozzi, M. K., Germanoska, B., et al. (2006). BOLD responses to stimuli: dependence on frequency, stimulus form, amplitude, and repetition rate. NeuroImage 31:585–599. doi: 10.1016/j.neuroimage.2005.12.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A., Loxley, P. N., O'Connor, S. C., and Rennie, C. J. (2001). Modal analysis of corticothalamic dynamics, electroencephalographic spectra, and evoked potentials. Phys. Rev. E 63:041909. doi: 10.1103/PhysRevE.63.041909

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A., Pagès, J. C., Gabay, N. C., Babaie, T., and Mukta, K. N. (2018). Neural field theory of perceptual echo and implications for estimating brain connectivity. Phys. Rev. E 97:042418. doi: 10.1103/PhysRevE.97.042418

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A., Rennie, C. J., and Rowe, D. L. (2002). Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. Phys. Rev. E 65:041924. doi: 10.1103/PhysRevE.65.041924

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A., Rennie, C. J., and Wright, J. J. (1997). Propagation and stability of waves of electrical activity in the cerebral cortex. Phys. Rev. E 56:826. doi: 10.1103/PhysRevE.56.826

CrossRef Full Text | Google Scholar

Robinson, P. A., and Roy, N. (2015). Neural field theory of nonlinear wave-wave and wave-neuron processes. Phys. Rev. E 91:062719. doi: 10.1103/PhysRevE.91.062719

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A., Sarkar, S., Pandejee, G. M., and Henderson, J. A. (2014). Determination of effective brain connectivity from functional connectivity with application to resting state connectivities. Phys. Rev. E 90:012707. doi: 10.1103/PhysRevE.90.012707

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P. A., Zhao, X., Aquino, K. M., Griffiths, J. D., and Sarkar, S. Mehta-Pandejee, G. (2016). Eigenmodes of brain activity: neural field theory predictions and comparison with experiment. NeuroImage 142:79–98. doi: 10.1016/j.neuroimage.2016.04.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryle, G. (1949). The Concept of Mind. Chicago: University Chicago Press.

Google Scholar

Schaefer, A., Kong, R, Gordon, E. M, Laumann, T. O, Zuo, X. N, Holmes, A. J., et al. (2018). Global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI. Cereb. Cortex 29, 3095–3114. doi: 10.1093/cercor/bhx179

PubMed Abstract | CrossRef Full Text | Google Scholar

Schiff, L. I. (1968). Quantum Mechanics, 3rd Edn. New York, NY: McGraw-Hill.

Google Scholar

Smith, S. M., Fox, P. T, Miller, K. L., Glahn, D. C., Mickle Fox, P., Mackay, C. E., et al. (2009). Correspondence of the brain's functional architecture during activation and rest. Proc. Nat. Acad. Sci. U.S.A. 106:13040–13045. doi: 10.1073/pnas.0905267106

PubMed Abstract | CrossRef Full Text | Google Scholar

Sporns, O. (2011). Networks of the Brain. Cambridge: MIT Press.

Google Scholar

Tagliazucchi, E., and Laufs, H. (2014). Decoding wakefulness levels from typical fMRI resting-state data reveals reliable drifts between wakefulness and sleep. Neuron 82:695–708. doi: 10.1016/j.neuron.2014.03.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Tewarie, P., Prasse, B., Meier, J. M., Santos, F. A. N., Douw, L., Schoonheim, M. M., et al. (2020). Mapping functional brain networks from the structural connectome: Relating the series expansion and eigenmode approaches. NeuroImage 216:116805. doi: 10.1016/j.neuroimage.2020.116805

PubMed Abstract | CrossRef Full Text | Google Scholar

van Albada, S. J., Kerr, C. C., Chiang, A. K. I., Rennie, C. J., and Robinson, P. A. (2010). Neurophysiological changes with age probed by inverse modeling of EEG spectra. Clin. Neurophysiol. 121, 21–38. doi: 10.1016/j.clinph.2009.09.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Van De Ville, D., Britz, J., and Michel, C. M. (2010). EEG microstate sequences in healthy humans at rest reveal scale-free dynamics. Proc. Nat. Acad. Sci. U.S.A. 107:18179–18184. doi: 10.1073/pnas.1007841107

PubMed Abstract | CrossRef Full Text | Google Scholar

Vidaurre, D., Hunt, L. T., Quinn, A. J., Hunt, B. A. E., Brookes, M. J., Nobre, A. C., et al. (2018). Spontaneous cortical activity transiently organises into frequency specific phase-coupling networks. Nat. Commun. 9:2987. doi: 10.1038/s41467-018-05316-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaesoubi, M., Miller, R. L., and Calhoun, V. D. (2015). Mutually temporally independent connectivity patterns: a new framework to study resting state brain dynamics with application to explain group difference based on gender. NeuroImage 107, 85–94. doi: 10.1016/j.neuroimage.2014.11.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeo, B. T. T., Krienen, F. M., Sepulcre, J., Sabuncu, M. R., Lashkari, D., Hollinshead, M., et al. (2011). The organization of the human cerebral cortex estimated by intrinsic functional connectivity. J. Neurophysiol. 106:1125. doi: 10.1152/jn.00338.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Yule, G. U. (1926). Why do we sometimes get nonsense-correlations between time-series?–A study in sampling and the nature of time-series. J. Roy. Stat. Soc. 89, 1–63. doi: 10.2307/2341482

CrossRef Full Text | Google Scholar

Zalesky, A., and Breakspear, M. (2015). Towards a statistical test for functional connectivity dynamics. NeuroImage 114, 466-470. doi: 10.1016/j.neuroimage.2015.03.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Zalesky, A., Fornito, A., and Bullmore, E. T. (2012). On the use of correlation as a measure of network connectivity. NeuroImage 60, 2096–2106. doi: 10.1016/j.neuroimage.2012.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Zwillinger, D. (1989). Handbook of Differential Equations. San Diego, CA: Academic Press.

Google Scholar