An electrophysiological validation of stochastic DCM for fMRI

Introduction

In addition to localizing brain regions that encode specific sensory, motor, or cognitive processes, contemporary neuroimaging research tries to understand how information is transmitted within brain networks (Sporns, 2010). Ultimately, the ambition is to understand the nature of the information that neuronal populations pass on to each other. This speaks to the notion of functional integration, which views cognitive function as resulting from information exchange within brain networks (Tononi et al., 1994). This means one has to understand the directed influence or effective connectivity that brain systems (e.g., cortical areas, neuronal populations, or single neurons) exert on each other. Analyzing effective connectivity rest on models that formalize assumptions about how neuronal systems are wired and how they respond to different stimuli. These models are then used to interpret brain responses measured using, e.g., functional magnetic resonance imaging (fMRI) or magneto-/electroencephalography (M/EEG). We refer to Valdes-Sosa et al. (2011) for a comprehensive review on effective connectivity, and its relationship with influence and causality.

In the context of M/EEG, detailed biophysical knowledge about the activity of neural ensembles has led the community to develop and validate dynamical models that can describe macro-scale dynamics in great detail (cf. Deco et al., 2008; Coombes, 2010; Bressloff, 2012). Most of these models are inspired by approaches in statistical physics based on the notion of a mean field, i.e., the idea that interactions among ensembles of neurons can be captured by summary statistics (i.e., moments of the relevant distribution). Dynamic causal modeling [DCM, see Daunizeau et al. (2011) for a recent review] for M/EEG embeds these models into a formal (Bayesian) statistical framework that allows for parameter estimation and model comparison when analyzing evoked or induced responses. In this context, the approach has proven successful in exploiting the realism of these biophysical models to capture experimental effects of cognitive manipulations in terms of network plasticity (e.g., Garrido et al., 2007 or Moran et al., 2011).

Even more established is the use of DCM for fMRI data, for which the approach was originally proposed by Friston et al. (2003), building on previous advances in the biophysics of hemodynamic processes (Buxton et al., 1998). Interestingly, fMRI paradigms elicit neural responses that unfold on a much slower time scale (seconds) than that typical of electrophysiological measurements (milliseconds). In contrast to hemodynamic processes, relatively little is known about the biophysical underpinning of these neuronal responses, which has precluded the use of detailed (realistic) models of neural dynamics in DCM for fMRI. Instead, the slow (and widespread) dynamics of interacting brain regions are captured by phenomenological models that are much simpler than those employed for M/EEG data. Even though these models serve as a reference when evaluating fMRI effective connectivity methods (Smith et al., 2011), their simplicity may impose limitations on the ensuing inferences (but see, for example, David et al., 2008 or Brodersen et al., 2011). This is why we have previously proposed variational Bayesian approaches that can deal with some of the (unavoidably simplifying) model assumptions by introducing random fluctuations of physiological states (Friston et al., 2008, 2010; Daunizeau et al., 2009b). Initial applications of this “stochastic” version of DCM [stochastic dynamic causal modeling (sDCM)] have established some face validity in the context of fMRI data analysis. For example, Friston et al. (2011) showed how sDCM can be used in combination with post-hoc model comparison (Friston and Penny, 2011) to explore large model spaces, and Li et al. (2011) examined the smoothness of the underlying neural fluctuations. In Daunizeau et al. (2012), we provide a comprehensive comparison of deterministic and stochastic DCM, in terms of parameter estimation and model comparison. In particular, we showed how accounting for random effects on the system's dynamics can improve network identification by exploiting the decorrelation of neural time series induced by neural noise.

However, sDCM for fMRI has its own methodological challenges and still lacks solid experimental validation. On the one hand, the ability of sDCM for fMRI to discriminate between the respective contributions of neuronal noise and measurement noise in the observed fMRI signal depends on modeling assumptions (Roebroeck et al., 2011; Daunizeau et al., 2012). On the other hand, it is yet unclear how to formally relate the neural states of DCM for fMRI to electrophysiological measures of activity. In this note, we use EEG data to provide empirical evidence for the predictive validity of sDCM for fMRI data. We therefore place this work in the context of the ongoing debate regarding the hemodynamic correlate of the EEG [see Rosa et al. (2010a) for a recent critical review of this literature]. Here, we appeal to neural field theory (Wilson and Cowan, 1972, 1973; Lopes da Silva et al., 1973; Nunez, 1974; Amari, 1975, 1977; Freeman, 1975; Jirsa and Haken, 1996) to revisit the heuristic hemodynamic correlate of EEG signals originally proposed in Kilner et al. (2005). In essence, this hypothesis suggests that the fMRI signal is correlated with the frequency modulation of the EEG (where larger fMRI signals reflect a shift of the EEG spectrum toward higher frequencies). This is based on the idea that slow dynamics underpinning fMRI responses reflect the instantaneous frequency of oscillating modes of activity (Cabral et al., 2011; Friston et al., 2011). Here, we show how slow macroscopic modes of activity that emerge from intrinsic (within region) connections may modulate the frequency spectrum of fast dynamical modes that contribute to the EEG signal. Guided by this notion, we obtain estimates of neural dynamics from fMRI data using sDCM and ask whether these estimates predict the spectral behavior of concurrent EEG data—in terms of the temporal changes in its center frequency (above and beyond physiological confounds). This enables us to establish the predictive validity of sDCM for fMRI data, in relation to electrophysiological responses.

This paper comprises three sections. In the first, we present the neural field model that motivates a separation of time scales and the ensuing heuristic hemodynamic correlate of the EEG. The second section is an empirical demonstration that serves to evaluate the predictive validity of stochastic DCM for fMRI for the ensuing electrophysiological correlates. We close with a discussion of the implications of our findings in the final section.

Model: Revisiting the Heuristic Hemodynamic Correlate of EEG

In this section, we first introduce a neural field model of a cortical region, and show how its spatiotemporal pattern of activity can be decomposed into canonical dynamic modes that have distinct time scales. In brief, slow macroscopic neural states control the rate of change of fast components. This motivates both the separation of time scales that underlies the generative model of sDCM for fMRI (Friston et al., 2011) and the heuristic hemodynamic correlate of EEG (Kilner et al., 2005). In brief, this section serves as a partial theoretical derivation of sDCM for fMRI, from the perspective of biologically-informed (neural field) models of macroscopic electrophysiological activity.

Neural Fields and Separation of Time Scales

We start with a description of the dynamics of single neurons within an ensemble, which we describe in terms of the probability density function p(ς) of neuronal post-synaptic membrane potentials (PSP) ς. We assume that neurons fire an action potential if their PSP surpasses a depolarization threshold χ. This means that 〈H(ς − χ)〉 represents the average (mean field) firing rate over the ensemble, where H is the Heaviside step function and 〈 〉 is the expectation under the ensemble density. Using the central limit theorem, one can show that the average firing rate over the ensemble is determined mainly by the first two moments of the PSP ensemble density p(ς) (Liley et al., 2002; Marreiros et al., 2008):

$\begin{array}{cc}\begin{array}{l}\langle H(\varsigma -\chi )\rangle ={\displaystyle {\int }_{\chi }^{\infty }p}(\varsigma )d\varsigma \approx S(\eta )=\frac{1}{1+\exp (-\rho (\eta -\chi ))}\\ \eta =\langle \varsigma \rangle \\ {\rho }^{-2}=\frac{3}{{\pi }^2}\langle {\left(\varsigma -\eta \right)}^2\rangle \end{array}& (1)\end{array}$

where S is a sigmoid function and ρ is proportional to the inverse of the standard deviation (dispersion) of depolarization within the ensemble. This is a statistical (mesoscopic) summary of neuronal ensemble activity, corresponding roughly to a cortical macrocolumn (about 100,000 neurons and 1 mm² of cortex). For our modeling purposes, this summary statistics is now taken to the macroscale, where each ensemble or column has a position r on the 2D-cortical manifold. At this macroscopic scale, states like depolarization can be regarded as a continuum or field (Jirsa and Haken, 1996; Liley et al., 2002; Coombes et al., 2007; Daunizeau et al., 2009a; Pinotsis et al., 2012), which is a function of space and time: η(t) → η(r, t). This means we can describe the local average input and output activity of neurons using the depolarization field η(r, t) and the sigmoid function. The spatiotemporal dynamics of the neural field essentially depends on how local ensembles influence each other—as described below.

We now consider a region of the cortical manifold that is defined by its continuous domain D, over which the neural field η(r, t) is defined. The field is considered to result from the temporal convolution of pre-synaptic input μ(r, t) as follows:

$\begin{array}{cc}\begin{array}{l}\eta (r,t)={\displaystyle {\int }_{-\infty }^{\infty }J}(t-t^{\prime })\mu (r,t)d{t}^{\prime }\\ =J\otimes \mu (r,t)\\ J(t)=\{\begin{array}{cc}\gamma \frac{t}{\tau }\exp \left(-\frac{t}{\tau }\right)& t\ge 0\\ 0& t<0\end{array},\end{array}& (2)\end{array}$

where J is the synaptic impulse response function, τ is a time-constant, and γ controls the maximum PSP following an action potential. The maximum PSP and time-constant summarize the microscopic properties of the ensemble (Jansen and Rit, 1995). Because action potentials propagate with a finite conduction velocity, the mean rate of arrival of pre-synaptic impulses can be expressed as a time-delayed integral (spatiotemporal convolution) of local firing rates:

$\begin{array}{cc}\begin{array}{l}\mu (r,t)={\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{D}G}}(|r-r^{\prime }|,t-t^{\prime })S(\eta (r^{\prime },t^{\prime }))d{r}^{\prime }d{t}^{\prime }\\ =G\otimes S°\eta (r,t),\end{array}& (3)\end{array}$

where D is the domain, over which the neural field is defined, and G is a homogeneous Green function that acts as a spatiotemporal convolution operator. Let r = |r − r′| be the distance between two points with positions r and r′. Below, we follow Bojak and Liley (2010) and parameterize the Green function as follows:

$\begin{array}{cc}G(r,t)=\frac{G^0}{4\pi }\frac{1}{{\sigma }^2}{t}^{-1}\exp \left(-\frac{r^2+{\overline{\nu }}^2{t}^2}{2\sigma \overline{\nu }t}\right)H(t),& (4)\end{array}$

where the Heaviside step function endows Equation (3) with temporal causality and G⁰ is the total number of synaptic connections. This dispersive propagator assumes a Gaussian falloff in the density of synaptic connections and accounts for a certain variability in velocities ($\overline{\nu }$ is the mode of the velocity distribution at “large distance,” i.e., when r/σ → ∞). The spatial scale of lateral connectivity is controlled by the spatial dispersion σ. Figure 1 depicts the dependence of the dispersive propagator with respect to space and time (see also Liley et al., 2012—Figure 3).

FIGURE 1

Figure 1. Neural field model: the dispersive propagator. This figure depicts the dispersive propagator of Equation (4) as a function of both time (x-axis) and distance (y-axis), in terms of the density 2πrG(r, t) of synaptic connections that are reached by an action potential emitted at r = t = 0 (on a logarithmic scale). Note that the 2πr scaling arises because of radial symmetry in 2D (Bojak and Liley, 2010). The white line shows the average distance, at which activity is propagated as a function of time. The black lines are contour lines of the connection density. These can be used to eyeball how fast the density attenuates (e.g., along the average distance path).

In addition to its plausible form, this homogeneous Green function has proven very useful, in that it has an analytical Fourier transform. Applying the Fourier transform back and forth to Equations (2) and (3) yields the following partial differential equation (see Bojak and Liley, 2010):

$\begin{array}{cc}\begin{array}{l}\left(\frac{\overline{\nu }}{2\sigma }+\frac{\partial }{\partial t}-\frac{\overline{\nu }\sigma }{2}{\nabla }^2\right)\mu (r,t)=\frac{G^0\overline{\nu }}{2\sigma }S(\eta (r,t))\\ \left({\tau }^2\frac{{\partial }^2}{\partial t^2}+2\tau \frac{\partial }{\partial t}+1\right)\eta (r,t)=\tau \gamma \mu (r,t).\end{array}& (5)\end{array}$

The first partial differential equation describes the diffusion of action potentials of spikes (through the Laplacian operator) and derives from the lateral connectivity kernel in the Green function. The second ordinary differential equation derives from Equation (2) and describes synaptic transmission; i.e., how propagated pre-synaptic firing is accumulated locally to produce depolarization.

The solution to Equation (5) can be decomposed into canonical dynamical modes, by projecting the neural field on the set of eigenvectors w^(k)(r) of the Laplacian operator defined over the domain D (see Daunizeau et al., 2009a). These spatial modes solve the Laplacian eigenvalue problem, i.e.,: ∇² w^(k) (r) = λ_k w^(k) (r), where the distribution of eigenvalues λ is such that: λ_k < λ₀ = 0. They form a complete and orthonormal basis function set for the neural field, allowing for an eigen-decomposition of the field, as follows:

$\begin{array}{cc}\begin{array}{l}\eta (r,t)={\displaystyle \sum _k{z}_1^{(k)}}(t)w^{(k)}(r)\\ \mu (r,t)={\displaystyle \sum _k{z}_2^{(k)}}(t)w^{(k)}(r),\end{array}& (6)\end{array}$

here, the dynamics of the spatial modes are defined as:

$\begin{array}{cc}\begin{array}{l}{z}_1^{(k)}(t)={\displaystyle {\int }_D{w}^{(k)}}(r)\eta (r,t)dr\\ z_2^{(k)}(t)={\displaystyle {\int }_D{w}^{(k)}}(r)\mu (r,t)dr.\end{array}& (7)\end{array}$

Furthermore, it can be shown that the dynamics of the spatial modes (Equation 5) can be approximated as third-order ordinary differential equations that are coupled through the fundamental mode k = 0 (see Appendix 1):

$\begin{array}{cc}\begin{array}{l}{\dot{z}}^{(k)}=\left[\begin{array}{c}{\dot{z}}_1^{(k)}\\ {\dot{z}}_2^{(k)}\\ {\dot{z}}_3^{(k)}\end{array}\right]\\ =\left[\begin{array}{ccc}0& 0& 1\\ c\left(z_1^{(0)}\right)& \frac{\overline{\nu }}{2\sigma }\left({\sigma }^2{\lambda }_k-1\right)& 0\\ -\frac{1}{{\tau }^2}& \frac{\gamma }{\tau }& -\frac{2}{\tau }\end{array}\right]z^{(k)}+O\left(z^2\right)\\ c\left(z_1^{(0)}\right)\equiv \frac{\rho G^0\overline{\nu }}{2\sigma }S\left(z_1^{(0)}\right)\left(1-S\left(z_1^{(0)}\right)\right)\end{array}& (8)\end{array}$

here, z₃ ≡ ż₁, and c(z⁽⁰⁾₁) is a non-linear function of the fundamental mode z⁽⁰⁾₁ that mediates the coupling between modes describing firing rate and depolarization. Roughly speaking, each mode decays at a rate that is proportional to $\sigma \overline{\nu }{\lambda }_k$. As a consequence, high propagation velocities in large cortical regions dissipate the spatial modes quickly, with the exception of the fundamental mode, which has a zero eigenvalue. More generally, smooth eigenmodes (low spatial scales) will be associated with slow dynamics (see, e.g., Schultze-Kraft et al., 2011). Note that the natural time scale of the fundamental mode can be very slow, if the conduction velocity $\overline{\nu }$ becomes small compared to the dispersion of lateral connections σ (i.e., if $\overline{\nu }/\sigma \to 0$). In this case, the fundamental mode z⁽⁰⁾ will decay so slowly that it will predominate over other stable (higher spatial frequency) modes z^(k), which disappear as soon as they are created. This means that Equation (8) effectively instantiates a separation of time scales, where the slowest mode “enslaves” the fast modes (Haken, 1983) through the coupling function c(z⁽⁰⁾₁). This slaving effect is such that the frequency profile of fast modes (k > 0) is controlled by the amplitude of slow modes (k = 0).

Modeling Local Activation in Terms of EEG Frequency Modulation

Equation (8) expresses the Jacobian of eigenmode dynamics as a function of the fundamental mode of activity z⁽⁰⁾₁ within a given brain region, over which the neural field is defined. This means that slow fluctuations in the fundamental mode z⁽⁰⁾₁ should be associated with changes in the EEG frequency spectrum P(ω), which derives from the Laplace transform of Equation (8) (see Appendix 2):

$\begin{array}{cc}\text{P}(\omega )\propto {\left|{\displaystyle \sum _{k\le k_c}\frac{\begin{array}{l}\left(i\omega +\frac{\overline{\nu }}{2\sigma }\left(1-{\sigma }^2{\lambda }_k\right)\right)\hfill \\ \left(i\omega +\frac{2}{\tau }+1\right)+\frac{\gamma }{\tau }\hfill \end{array}}{\begin{array}{l}\left(i\omega +\frac{\overline{\nu }}{2\sigma }\left(1-{\sigma }^2{\lambda }_k\right)\right)\hfill \\ \left(i\omega \left(i\omega +\frac{2}{\tau }\right)+\frac{1}{{\tau }^2}\right)-\frac{\gamma }{\tau }c\left(z_1^{(0)}\right)\hfill \end{array}}}\right|}^2& (9)\end{array}$

where k_c is the cut-off order induced by the skull's low electrical conductivity. Here, the summation is over the fast modes (up to k_c = 16), and we have assumed that EEG responses reflect postsynaptic potentials, i.e., z^(k)₁. Similarly to Kilner et al. (2005), let us now define the EEG center frequency $\overline{\omega }$ as the first moment of the (normalized) frequency power spectrum $\overline{\text{P}}(\omega )$:

$\begin{array}{cc}\begin{array}{l}\overline{\text{P}}(\omega )=\frac{\text{P}(\omega )}{{\displaystyle \underset{{\omega }_0}{\overset{\infty }{\int }}\text{P}}(\omega )d\omega }\\ \overline{\omega }\equiv {\displaystyle {\int }_{{\omega }_0}^{\infty }\omega }\overline{\text{P}}(\omega )d\omega \end{array}& (10)\end{array}$

By construction, normalized frequency power $\overline{\text{P}}(\omega )$ is the proportion of total power attributed to any given frequency of interest. Thus, one can think of it in terms of a probability density function over frequencies, whose first-order moment would be the center frequency. Intuitively, an increase in the center frequency $\overline{\omega }$ signals a shift of frequency power toward higher frequencies. Note that the center frequency is not the main peak in the frequency spectrum. Instead, variations of $\overline{\omega }$ over time quantify frequency modulations of the EEG signal. The key idea here is that the EEG center frequency $\overline{\omega }$ is a function of the slow eigenmode dynamics z⁽⁰⁾₁, and its modulation over time thus follows the slow eigenmode dynamics z⁽⁰⁾₁, i.e., $\overline{\omega }=\overline{\omega }\left(z_1^{(0)}\right)$.

Figure 2 depicts the frequency modulation of the EEG signal induced by changes in the slow eigenmode z⁽⁰⁾₁. The frequency power spectrum P (ω) given in Equation (9) was evaluated using the parameter values given in Table 1 below. One can see that when the slow eigenmode tends toward the action-potential firing threshold (z⁽⁰⁾₁ → χ), there is both a power increase in the low frequencies and a decrease in the high frequencies (cf. increase in the delta/theta/alpha band, and decrease in the beta/gamma band). This shift toward lower frequencies can be seen most clearly on the normalized frequency power $\overline{\text{P}}(\omega )$. As a consequence, the center frequency $\overline{\omega }$ increases as the slow eigenmode gets further away from the action potential firing threshold. In our example, this induces a change of about 10 Hz in the EEG center frequency. Note that the precise numerical value of the center frequency $\overline{\omega }$ (as well as its susceptibility to the eigenmode z⁽⁰⁾₁) depends upon the neural field's parameters (cf. Table 1). In fact, using more realistic neural fields models (accounting for, e.g., excitatory and inhibitory interneurons, see Pinotsis et al., 2012) would profoundly change the shape of the frequency power spectrum. However, the qualitative effect of neural activation upon the center frequency is qualitatively invariant to such changes. We will comment on this in the “Discussion” section.

FIGURE 2

Figure 2. Neural field model: EEG frequency modulation. This figure shows the effect of changes in the mean membrane depolarization fundamental eigenmode z⁽⁰⁾₁ on the frequency content of the EEG signal. Upper-left: EEG frequency spectrum P(ω) (z-axis) as a function of frequency ω (x-axis) and eigenmode amplitude z⁽⁰⁾₁ (y-axis). Upper-right: Normalized frequency spectrum $\overline{\text{P}}(\omega )$ (cf. Equation 10, same format). Lower-left: Magnitude of standard EEG frequency bands (y-axis; blue: delta, green: theta, red: alpha, magenta: beta, violet: gamma) as a function of the eigenmode amplitude z⁽⁰⁾₁ (x-axis). Lower-right: Centre frequency $\overline{\omega }$ (y-axis) as a function of the eigenmode amplitude z⁽⁰⁾₁ (x-axis). Note that here, the action potential firing threshold χ is 30 mV.

TABLE 1

Table 1. Neural field parameters.

The Link Between EEG Time-Frequency Response and fMRI Temporal Response

In this paper, we assume that the slow macroscopic dynamics x⁽ⁿ⁾ underlying fMRI data in a given region of interest follow the fundamental eigenmode z⁽⁰⁾₁ of its membrane depolarization cortical field. The motivations for this assumption include the following:

• Its relative simplicity.
• Sensitivity to input versus output spiking activity: there is a consensus now on the fact that mean membrane depolarization (as opposed to e.g., firing rate) drives hemodynamic changes [see Logothetis and Wandell (2004) for a review on the relationship between local field potentials and the BOLD signal).
• Energy budget: both the slow time scale and the low spatial frequency of the fundamental eigenmode make it likely to be the main driver of induced BOLD changes, when compared with faster (and higher spatial frequency) eigenmodes. This is because blood flow increases in local vascular trees are induced by coherent spatio-temporal summation of neural activity (see e.g., Daunizeau et al., 2010; Rosa et al., 2010a,b; Schmuel, 2010).
• It relates to theoretical accounts on critical slowing of high-dimensional dynamical systems, which has recently been used to motivate the rate of change of neural states in stochastic DCM for fMRI (Friston et al., 2010).

We will further comment on this point in the “Discussion” section.

In the previous section, we have described how a separation of time scales could arise from local interactions at the level of a cortical region. It turns out that this result can be generalized to a set of coupled local neural fields, in that the fundamental mode mediates the extrinsic (between-region) coupling. In other words, the influence different neural fields exert on each other expresses itself through the coupling of their respective fundamental modes, which locally drive higher order modes (as in Equation 8). Deriving the precise dynamical properties of coupled neural fields is beyond the scope of the present study. However, one can invoke the slaving effect when considering coupled (distal) brain regions, the network properties of which shape the temporal dynamics of local fundamental eigenmodes z⁽⁰⁾₁. This means that long-range connections control the frequency profile of brain regions responses (see e.g., Bojak et al., 2011), through their effect on their respective fundamental eigenmodes z⁽⁰⁾₁. Since we did not consider coupled neural fields in the above treatment, we replace their (unknown) evolution function with a first-order Taylor expansion on slow neural states, yielding the following Langevin equation for stochastic DCM (see Daunizeau et al., 2012):

$\begin{array}{cc}{\dot{x}}^{(n)}\approx A{x}^{(n)}+v& (11)\end{array}$

where θ⁽ⁿ⁾ = {A_ij} are neuronal parameters that measure extrinsic coupling strengths. Here, we have neglected non-linear interaction terms that act as gating factors on the network connections x⁽ⁿ⁾ (cf. Stephan et al., 2008). The influence of the fast modes on the motion of slow modes x⁽ⁿ⁾ is expressed via fluctuations v, which are comparatively much faster. This allows us, in the context of sDCM for fMRI, to treat v as stochastic (neural) noise and place priors p(v|m) on its generalized motion (see below).

In contradistinction, typical evoked EEG responses are driven by a mixture of slow and fast modes. Note that the relative contribution of the slow and fast modes is imbalanced (toward fast modes), simply because the summation span all eigenmodes up to cut-off order k_c. In fact, typical evoked EEG responses disappear within a second or so of peristimulus time; i.e., their limit frequency has an order of magnitude similar to the fMRI sampling rate (about 1 Hz).

Taken together, we expect to see concurrent changes in the slow neural states x⁽ⁿ⁾ that drive BOLD and in the EEG power spectrum P (ω). However, from the above section, we know that the underlying relationship is non-linear and depends upon intrinsic properties of local neural fields, such as the size of active brain regions and the spatial decay rate of synaptic connections (cf. Equation 9). In addition, the respective contribution of each region to the scalp EEG power spectrum is virtually unknown, as this would require solving the so-called EEG inverse problem (see Rosa et al., 2010a,b). These issues conspire with the unavoidable simplifications of our model to make the prediction of the full EEG frequency spectrum practically irrelevant. In the following, we will thus focus on the EEG center frequency $\overline{\omega }$, which captures global changes in the frequency spectrum. For the sake of simplicity, in the remainder of this work, we will invoke a first-order Taylor expansion of the following form:

$\begin{array}{cc}\overline{\omega }\left(x^{(n)}\right)=\overline{\omega }(0)+{\frac{\partial \overline{\omega }}{\partial x^{(n)}}|}_0{x}^{(n)}+O\left(x^2\right)& (12)\end{array}$

where the gradient ${\partial \overline{\omega }/\partial x^{(n)}|}_0$ measures the susceptibility of the center frequency to the slow neural states. This quantity will be directly estimated from experimental data (along with the intercept $\overline{\omega }(0)$, which acts as a confound). Testing for Equation (12) effectively bypasses the dependency of the EEG power spectrum to the respective contribution each local neural field included in the network, as well as their biophysical properties, and only retains the qualitative prediction of the model.

Summary

In this section, we have used the neural field formulation of local neural activity to demonstrate an emergent separation of temporal scales, in which the slow fluctuations of modes or patterns of depolarization enslave faster modes. Crucially, the amplitude of slow modes (which we associate with fMRI signals) and modulate the frequencies of fast modes (which we associate with EEG signals). This dynamical behavior is based upon a plausible mean field approximation to macroscopic neuronal activity on the cortical surface and, within the setting, provides a formal verification of the heuristic in Kilner et al. (2005), linking fMRI signals to changes in the center frequency of electrophysiological signals. In the next section, we will use stochastic DCM to obtain estimates of macroscopic (slow) neural fluctuations from fMRI data. We then test whether these estimates can predict fluctuations in concurrently measured EEG data—in terms of its center frequency, as suggested by the neural field treatment above (Equations 9–12).

Empirical Data: Resting State and Generalized Spike and Wave (GSW) Activity

In this section, we conduct an sDCM analysis of data from patients with epilepsy who exhibited generalized spike and wave (GSW) activity, while undergoing concurrent EEG/fMRI measurements in the scanner. We first describe the EEG/fMRI acquisition and conventional activation analysis. We then describe the sDCM analysis of fMRI data and ensuing validation using EEG data.

EEG/fMRI Acquisition and Activation Analysis

These data were part of a previous neuroimaging study of idiopathic generalized epilepsy (Hamandi et al., 2006; Vaudano et al., 2009), whose acquisition protocol and data pre-processing are briefly summarized here.

Ten-channel EEG (10–20 system) was recorded using MR-compatible equipment, along with bipolar electrocardiogram. After filtering and artifact correction (Krakow et al., 2000), the onsets and offsets of GSW activity were identified by two experts (see Vaudano et al., 2009 for details). Seven hundred T2*-weighted single-shot gradient echo echo-planar images (TE = 40 ms, TR = 3 s, 21 interleaved axial slices of 5 mm thickness, FOV = 24 × 24 cm², 64 × 64 matrix) were acquired over a 35 min session with a 1.5 Tesla MRI scanner (Horizon EchoSpeed, General Electric). Patients were asked to rest with their eyes closed and to keep still. FMRI data were pre-processed using SPM8 (http://www.fil.ion.ucl.ac.uk/spm/). EPI time series were realigned and spatially smoothed with an 8 mm FWHM isotropic Gaussian kernel and spatially normalized to the standard anatomical space. For each patient, a general linear model (GLM) was constructed to test for the presence of regional GSW-related BOLD changes. Periods of GSW activity were modeled as blocks, beginning at GSW onset and terminating at their offset. The GSW regressors were then convolved with the canonical hemodynamic response function (plus temporal and dispersion derivatives) before inclusion in the GLM. In addition, we included slow drifts (Fourier basis function set), as well as motion-related effects (head and eye movements), cardiac-related effects (see Liston et al., 2006) and scan-nulling regressors (modeling inter-scan motion events larger than 0.2 mm, cf. Lemieux et al., 2007) as confounding factors.

On average, the GLM design matrices contained about 100 confounds regressors. Significant positive and negative GSW-related BOLD responses were identified by means of an F-contrast on the GSW regressors for the nine patients included in this study. SPMs were thresholded at p < 0.05 (FWE whole-brain corrected) to define three regions of interest (ROIs), which were implicated in the initiation and termination of GSW discharges in all patients: thalamus, prefrontal cortex (PFC), and precuneus. A summary time series was derived for each ROI by computing the first eigenvariate of all suprathreshold voxel time series within a 10 mm of the ROI centers. The time series were adjusted for all confounding effects included in the GLM analysis.

Validating the sDCM Macroscopic Neural States Estimates with EEG: Methods

In addition to the neural evolution function given in Equation (11), DCM for fMRI requires the specification of an additional set of hemodynamic states that couple neural dynamics to observed BOLD signal changes:

$\begin{array}{cc}{\dot{x}}^{(h)}=\left[\begin{array}{c}{x}^{(n)}-{\kappa }_s{x}_1^{(h)}-{\kappa }_f\left(e^{x_2^{(h)}}-1\right)\\ x_1^{(h)}{e}^{-x_2^{(h)}}\\ \frac{1}{{\tau }_0}\left(e^{x_2^{(h)}}-e^{x_3^{(h)}/\alpha }\right)e^{-x_3^{(h)}}\\ \frac{1}{{\tau }_0}\left(\frac{e^{x_2^{(h)}-x_4^{(h)}}\left(1-{\left(1-E_0\right)}^{e^{-x_2^{(h)}}}\right)}{E_0}-e^{\left(1-\alpha \right)x_3^{(h)}/\alpha }\right)\end{array}\right]& (13)\end{array}$

here, x⁽ⁿ⁾ represents regional neural activity, whose dynamics is given in Equation (11). Equation (13) expresses changes in hemodynamic states x^(h), as a response to a neural perturbation x⁽ⁿ⁾ and as a function of hemodynamic parameters θ^(h) = {κ_s, κ_f, α, τ₀, E₀, ε₀} (for details, see the appendix in Stephan et al., 2008). Finally, one has to specify the observation mapping from hemodynamic states x^(h) to observed local BOLD changes y (Stephan et al., 2007):

$\begin{array}{cc}\begin{array}{l}y=V_0(4.3{\nu }_0{E}_{0}TE\left(1-e^{x_4^{(h)}}\right)+{\epsilon }_0{r}_0{E}_{0}TE\left(1-e^{x_4^{(h)}-x_3^{(h)}}\right)\\ +\left(1-{\epsilon }_0\right)\left(1-e^{x_3^{(h)}}\right))+\epsilon ,\end{array}& (14)\end{array}$

where ε is an additive measurement noise, and the parameters are defined in Table 2 below.

TABLE 2

Table 2. Prior means of hemodynamic parameters and acquisition-related constants for fMRI at 1.5 T.

In brief, the predicted data y depend non-linearly on the unknown model variables Θ = {x, θ}, through Equations (11), (13), and (14). This model (as well as statistical assumptions about measurement noise ε) is encoded in the likelihood function p(y|Θ, m). Priors p(Θ|m) specify our assumptions about the magnitude of state noise, evolution parameters and observation parameters—where the prior means of the hemodynamic parameters of Equations (13) and (14) are given in Table 2:

Inverting the generative model m means (1) approximating the conditional density p(Θ|y, m) of unknown variables Θ given the set of sampled measurements y and (2) quantifying the model evidence p(y|m). Non-linearities in the generative model eschew exact analytical solutions to this problem, which is finessed using variational approaches that rest on optimizing a free-energy lower bound F(q) to the model evidence, with respect to an approximate conditional density q(Θ):

$\begin{array}{cc}\begin{array}{l}F(q)={\langle \ln p\left(\Theta |m\right)+\ln p\left(y|\Theta ,m\right)-\ln q\left(\Theta \right)\rangle }_q\\ =\ln p\left(y|m\right)-D_{\text{KL}}\left(q(\Theta );p(\Theta |y,m)\right),\end{array}& (15)\end{array}$

where D_KL is the Kullback-Leibler divergence and the expectation 〈 〉_q is taken under q. From Equation (15), maximizing the functional F(q) with respect to q minimizes the Kullback-Leibler divergence between q(Θ) and the exact posterior p(Θ|y, m). This decomposition is complete in the sense that if q(Θ) = p(Θ|y, m), then F(q) = ln p(y|m). Typically, the iterative maximization of free energy proceeds under the Laplace approximation, where the approximate posterior q(Θ) ≈ p(Θ|y, m) is assumed to have a Gaussian form (see Friston et al., 2007). We refer to Daunizeau et al. (2009b) for details about the application of the variational Bayesian approach to stochastic DCM.

Following model inversion as described above, inference on hidden states is based on the conditional density q(x) and thus depends upon the generative model m. Usually, in DCM, inference on states or parameters is preceded by model selection. However, our focus was on the validity of sDCM for fMRI with regard to predicting the spectral behavior of concurrently acquired EEG data—based on the idea that slow macroscopic neural states x⁽ⁿ⁾ control the EEG frequency modulation. We therefore endowed the network with full reciprocal connectivity (no zero entry in matrix A of Equation 11) and simply applied the above variational Bayesian approach to derive the first moment $\widehat{x}={\langle x^{(n)}\rangle }_q$ of the conditional density q(x). In the following, $\widehat{x}$ thus refers to the sDCM estimate of the slow neural states, given fMRI time series y (and under model m).

To assess the relationship between EEG frequency modulation and neural states as estimated by sDCM, we constructed a GLM H₁, whose dependent variable was the trajectory $\overline{\omega }$ of the EEG center frequency sampled at each fMRI time sample and whose independent variables were the neural states estimated by sDCM plus confounds (cf. Equation 12):

$\begin{array}{cc}{H}_1:\overline{\omega }=\left[\widehat{x}{X}_0\right]\left[\begin{array}{c}\beta \\ {\beta }_0\end{array}\right]+e& (16)\end{array}$

where $\widehat{x}$ are the sDCM conditional estimates, X₀ contains the confounds (see below), β and β₀ are unknown regression coefficients, and e are i.i.d. residuals (with zero mean and unknown variance). Parameters β effectively capture the susceptibility of the EEG center frequency with respect to the slow neural states (cf. gradient ${\partial \overline{\omega }/\partial x^{(n)}|}_0$ in Equation 12). They lump together the contributions of the multiple regions that constitute the network as well as biophysical properties of the local neural fields. Parameters β₀ weight the contribution of the confounding factors to changes in the EEG center frequency. We then tested whether the sDCM neural states predicted EEG frequency modulation above and beyond what could be explained using the confounds. To this end, we used Bayesian model comparison at the group level (Stephan et al., 2009) to quantify the evidence in favor of the full model H₁, relative to a reduced model H₀ that only included the confounds:

$\begin{array}{cc}{H}_0:\overline{\omega }=X_0{\beta }_0+e& (17)\end{array}$

The assessment of sDCM predictive validity thus depends on the definition of the confounds, i.e., uncontrolled sources of variations in the EEG center frequency. In order to investigate the robustness of the statistical relationship between the slow eigenmode dynamics $\widehat{x}$ and the EEG center frequency trajectory $\overline{\omega }$, we have defined two different sets of confounds:

• Slow drifts: this assumes that the EEG center frequency undergoes slow and unspecific variations due to e.g., skin conductance changes, residual eye blinks artifacts, etc. …. Under this assumption, we set up X₀ similarly to the minimal confounds of a typical fMRI activation analysis, i.e., a Fourier basis function set (16 first harmonics; see Figure 6).
• Full confounds: this generalizes the intuition that confounds on the EEG center frequency are similar to those included in the fMRI activation analysis. For example, cardiac and respiratory effects, as well as body movements, could be modeled directly, as potential sources of (nuisance) variations. Under this assumption, we set up X₀ identically to the full confounding matrix used for each subject fMRI activation analysis (cf. section “EEG/fMRI Acquisition and Activation Analysis” above).

In brief, testing for the contribution of sDCM neural estimates $\widehat{x}$ on variations of the center frequency $\overline{\omega }$, above and beyond the full confounds is the most conservative of the two alternatives. However, it will be the least sensitive of both models as well; recall that β has only n_ROI entries, where n_ROI = 3 is the number of ROIs included in the sDCM analysis (to be compared with the dimension of β₀). In other words, these two sets of confounds can be motivated on principled grounds. We will thus perform a statistical comparison between H₀ and H₁ that marginalizes over the two sets of confounds.

Validating Neural States Estimate Obtained by sDCM Against EEG

Figure 3 shows the SPM of significant GSW-related BOLD changes (above and beyond physiological confounds), for the first patient of the study (F-test, p < 0.05, whole-brain corrected). This was used to define the three ROIs included in the sDCM analysis (Figure 4).

FIGURE 3

Figure 3. Absence seizure analysis: regions of interest. This figure summarizes the standard SPM activation results of a single epileptic patient with (petit mal) absence seizure (the same subject as in Figure 1). Significant (whole brain FWE-corrected, p < 0.05) positive and negative GSW-related BOLD responses were identified by means of an F-contrast on the GSW regressors. The color bar indicates the range of displayed F values.

FIGURE 4

Figure 4. Absence seizure analysis: sDCM analysis. This figure summarizes the sDCM analysis of a single subject fMRI data (same patient as in Figure 1). Upper-left: Time series of the observed and fitted BOLD signal in prefrontal cortex (PFC) as a function of time. The red shaded area indicates peri-GSW activity. Upper-right: Observed (y-axis) versus fitted (x-axis) BOLD signal in PFC. Lower-left: Macroscopic sDCM neural states in the three regions of interest (PF, prefrontal cortex; pC, precuneus; Th, thalamus) during peri-GSW activity, as a function of time. The red bar indicates strong head motion, which was modeled using scan-nulling regressors (see main text). Lower-right: sDCM conditional density on model parameters (left: first-order moment and ensuing network connectivity, right: posterior correlation matrix). Note: identifiability issues between pairs of parameters manifest as high posterior correlations.

One can see that most of the variance of the BOLD response is explained by the stochastic DCM. Focusing on peri-GSW periods, one can see that the conditional density q(x) shows high temporal variability, with increased uncertainty during head motion. This is expected, as scan-nulling confound regressors effectively prevent sDCM from deriving any information from fMRI data about hidden states during these periods. It is interesting to note that the network connectivity coefficients (elements of the A matrix in Equation 11) are small (about 0.02 Hz). This means that the average coupling of these regions (over the entire recording) is rather weak. Importantly, with the (partial) exception of the inhibitory intrinsic self-connections, the posterior correlation matrix (Figure 4) shows that most system parameters are identifiable and not confounded by hemodynamic parameter estimation. Generally, Figure 4 summarizes the motivation for this work—is the temporal variability of sDCM estimates of neuronal activity an artifact of measurement noise, or does it have an electrophysiological underpinning?

The cross-validation of neural state estimates by sDCM from fMRI data against EEG data requires the extraction of the time-dependent center frequency $\overline{\omega }$ for each channel of each subject, according to Equation (10). This is depicted in Figure 5, for the EEG channel C4 of the same patient in Figures 3, 4.

FIGURE 5

Figure 5. Absence seizure analysis: derivation of the EEG frequency modulation. This figure depicts the extraction of the EEG center frequency in a single subject (same patient as in Figure 1). Left: The EEG set-up used during the recording session is shown superimposed on the brain and skin surfaces (sensor C4 is highlighted). Right: Normalized square root power of windowed Fourier transform (z-axis) of the EEG traces of sensor C4 (x-axis: scanning time, y-axis: instantaneous frequency). The blue line shows the center frequency $\overline{\omega }$ as a function of scanning time t (cf. Equation 10 in the main text).

Recall that comparisons between the experimental and theoretical profiles of the frequency spectrum show differences that prevent meaningful quantitative fitting of Equation (9). For example, the model cannot account for peaks in the alpha band that appear in occipital channels (at least for the neural field parameter values of Table 1; results not shown). However, it is reassuring to see that the variations of the EEG frequency profile over time appear as small perturbations around a largely invariant average power spectrum. More precisely, the appearance of GSW activity does not dramatically distort its shape, but rather induces frequency modulations of small amplitude. This is important, since this affords face validity to the qualitative prediction of our neural field model. Interestingly, there is a clear change in the normalized frequency power beginning at the onset of the second GSW crisis. This leads to a marked decrease in the EEG center frequency (the first GSW event seems to be associated with a similar but weaker effect). We will briefly comment on GSW-related frequency modulations in the discussion.

Figure 6 summarizes the relationship between slow macroscopic sDCM states and the EEG frequency modulation of channel C4 for the same patients as in previous figures. Even though one can see that the center frequency $\overline{\omega }$ is very well fitted by the model (upper row of Figure 6: coefficient of determination: R² = 0.59), one might think that most of the frequency modulation could be explained by confounds. For example, eyeballing the trajectory of the observed EEG center frequency reveals clear low frequency tends that are likely to be captured by slow drifts confounds. To disclose the specific contribution of neural state estimates, we have plotted the data and model prediction after removing the confounds (adjusted data). One can see that, when considering slow drifts confounds, sDCM neuronal estimates explain EEG frequency modulations of about 10 Hz. This is a similar order of magnitude than center frequency variations captured by slow drifts (cf. parameter estimates on lower-right panel of Figure 6). Note that the corresponding correlation coefficient is about 0.49 (adjusted R² = 0.24), which indicates a reasonably good fit (above and beyond confounds).

FIGURE 6

Figure 6. Absence seizure analysis: EEG results (1). This figure summarizes the results of the analysis on the EEG frequency modulation in a single subject (same patient as in Figure 5) at channel C4. Upper-left: Observed (blue dashed line) and predicted (black plain line) frequency modulation (y-axis) of the best EEG channel (C4) as a function of scanning time (x-axis). For this patient, the two GSW episodes are indicated by the orange arrows. Upper-middle: Observed (y-axis) versus predicted (x-axis) frequency modulation for EEG channel C4. Upper-right: Matrix of confounds X₀ included in the analysis (slow drifts). Lower-left: Observed (blue dashed line) and predicted (black plain line) frequency modulation (y-axis) of the best EEG channel (C4) as a function of scanning time (x-axis), after adjustment for confounds. Lower-middle: Adjusted (y-axis) versus predicted (x-axis) frequency modulation for EEG channel C4. Lower-right: GLM parameter estimates, ± one standard deviation (orange: β, gray: β₀).

Figure 7 summarizes the result of Bayesian model comparison, in terms of the log-Bayes factor $\text{LBF}=\log p(\overline{\omega }|H_1)-\log p(\overline{\omega }|H_0)$, across all EEG channels for the same patient as before. One can see that all EEG channels show strong evidence in favor of the full model (H₁), whereas there is evidence against H₁ only for EOG and ECG channels. This suggests that the neuronal fluctuations estimated by stochastic DCM predict the EEG frequency modulation beyond slow drifts confounds. For completeness, we also show the results of classical inference (F-test with significance level p = 0.05, full model H₁ against the null H₀). As this test was repeated independently for each of the 30 channels, we applied a Bonferroni correction for multiple comparisons across channels. One can see that the profile of F-statistics and log-Bayes factors (LBFs) across channels is very similar.

FIGURE 7

Figure 7. Absence seizure analysis: EEG results (2). This figure summarizes the results of the analysis on the EEG frequency modulation in a single subject (same patient as in Figures 5, 6) across channels. Upper panel: Value of the F-statistic (y-axis) when testing for the significance of the contribution of the sDCM neural states dynamics in the EEG frequency modulation, as a function of EEG channels (x-axis). Red stars indicate channels that pass the 5% false positive rate threshold (with Bonferroni correction of the multiple comparisons across channels). Upper-right inset: Associated p-value across EEG channels (the red line indicates the corrected 5% threshold). Lower panel: Log Bayes factor log p(H₁|y) − log p(H₀|y) showing evidence in favor of H₁ versus H₀ across EEG channels (see main text). Red lines indicate posterior probabilities of model H₁ of 95% (upper line) and 5% (lower line), respectively.

Recall that the results summarized in Figures 6, 7 are only shown to illustrate the analysis. In fact, it is important to note that they (1) depend upon the definition of the confounds (here slow drifts) and (2) do not capture inter-individual variability. The question of whether the accuracy with which the neural states inferred by sDCM are able to predict the frequency modulations of the EEG signal (above and beyond confounds) generalizes at the group level is addressed below.

This Bayesian model comparison for the single patient shown in Figure 7 was then repeated for all subjects, for both sets of confounds (namely: slow drifts and full confounds). First, we pooled evidence over EEG channels for each subjects. This was simply done by summing the log- model evidences over EEG channels (within-subjects fixed effect. The ensuing distribution of LBFs across subjects—shown in Table 3 below—indicates that most subjects show very strong evidence in favor of H₁, irrespective of the particular set of confounds.

TABLE 3

Table 3. Log-Bayes factors (LBF) of H₁ versus H₀, for all subjects, and both sets of confounds.

To confirm this, we ran a random-effect analysis (Stephan et al., 2009) to assess the prevalence of H₁ and H₀ at the group level. Note that the two sets of confounds (“slow drifts” and “full confounds”) induce the four following models:

• H^(drifts)₁ : H₁ with slow drifts
• H^(full)₁ : H₁ with full confounds
• H^(drifts)₀ : H₀ with slow drifts
• H^(full)₀ : H₀ with full confounds.

The comparison of interest is in terms of whether or not adding the sDCM slow neural states adds anything to the prediction of the EEG center frequency, above and beyond confounds. Thus, we applied group-level family inference (Penny et al., 2010) to evaluate the evidence in favor or against H₁, irrespective of the definition of the confounds. We have thus partitioned the set of four models above into two families ${\overline{H}}_1$ and ${\overline{H}}_0$ that gather both sorts of confounds, as follows: ${\overline{H}}_1=\left\{H_1^{(\text{drifts})},H_1^{(\text{full})}\right\}$ and ${\overline{H}}_0=\left\{H_0^{(\text{drifts})},H_0^{(\text{full})}\right\}$. Deriving the relative evidence of ${\overline{H}}_1$ against ${\overline{H}}_0$ effectively marginalizes over the different confounds. First of all, we checked that there was a random effect at all, i.e.,: a difference between the respective prevalence of all models within the population. We found that the posterior probability of all models being equally frequent was p = 0.03. This Bayesian omnibus test indicates that the distribution of Bayes factors shown in Table 3 is unlikely to be due to chance. Figure 8 summarizes the group-level family inference results. One can see that the estimated prevalence of family ${\overline{H}}_1$ is about 0.82. In addition, the exceedance probability of ${\overline{H}}_1$ (versus ${\overline{H}}_0$) was 0.991. In other words, we can be 99% confident that sDCM estimates of neuronal fluctuations predict concurrent EEG frequency modulation above and beyond physiological confounds (in more than half the population).

FIGURE 8

Figure 8. Absence seizure analysis: EEG results (3). This figure summarizes the results of the analysis on the EEG frequency modulation at the group level (across subjects). Left: Posterior estimates of the expected frequency of families ${\overline{H}}_0$ and within the population (± one posterior standard deviation). Right: Posterior density (y-axis) over the expected frequency (x-axis) of family ${\overline{H}}_1$ within the population. The red shaded area indicates the mass of probability that is beyond 50% prevalence, yielding an exceedance probability of about 99%.

For completeness, we have also performed bayesian group-level model comparison, conditional on the definition of confounds. The results match those of the family-inference analysis above. In brief, the exceedance probability of H^(drifts)₁ (versus H^(drifts)₀) was 0.997, and the exceedance probability of H^(full)₁ (versus H^(full)₀) was 0.990. However, the bayesian omnibus test for these analyses was not as conclusive (slow drifts: p = 0.11, full confounds: p = 0.14). Taken together, this means that there is statistical evidence in favor of H₁ (versus H₀), irrespective of the definition of the confounds. This concludes our quantitative assessment of the predictive validity of sDCM for fMRI with regard to concurrently obtained EEG measures.

Discussion

Introducing random neural fluctuations in DCM for fMRI was originally motivated by the need to estimate hidden states in the absence of experimental perturbations. In this note, we have assessed the predictive validity of sDCM estimates of slow macroscopic neural states with respect to EEG frequency modulation. We have revisited the heuristic hemodynamic correlate of the EEG proposed in Kilner et al. (2005), using a neural field formulation to show how slow modes of local electrophysiological activity (which we assume drive changes in the fMRI signal) enslave the frequency spectrum of fast modes that contribute to the EEG signal. From this theoretical perspective, we tested whether slow neural estimates obtained by sDCM could explain concurrent EEG frequency modulation above and beyond physiological confounds.

To our knowledge, this work is the first attempt to provide empirical evidence that neural fluctuations inferred using stochastic DCM from fMRI time series have an electrophysiological underpinning. From an experimental perspective, this endorses the use of sDCM for studying spontaneous fluctuations that are beyond experimental control. Examples of this include—but are not limited to—multistable perception (Hesselmann et al., 2010) and epileptic events (e.g., absence seizures, cf. Vaudano et al., 2009). In this context, the added-value of sDCM is that it allows for comparing quantitative hypotheses about the propagation and/or the emergence of these phenomena within brain networks.

The acute reader might question the relevance of the modeling part of our paper, in the sense that quantitative predictions of the model were not used in the subsequent data analysis. However, we believe this model may contribute to the debate regarding the hemodynamic correlate of the EEG. First, it bears a number of interesting differences, with respect to the heuristic model of Kilner et al. (2005):

• Here, the main prediction, namely that fMRI activation is associated with a frequency modulation of the EEG signal, only relies on the assumption that BOLD changes are mostly driven by the fundamental (slow and smooth) eigenmode of the membrane potential field. This differential sensitivity of EEG and fMRI with respect to higher-order canonical eigenmodes bypasses most modeling assumptions of the heuristic model.
• Here, the relationship between activation and the EEG center frequency is not monotonic. More precisely, the center frequency decreases as one approaches the action potential threshold, either from hyperpolarized states (activation) of from depolarized states (de-activation). In brief, this form of critical slowing is qualitatively equivalent to the heuristic model, above the action potential threshold.

Second, we believe that the above qualitative theoretical prediction is very robust to model assumptions. First, we checked that the effect persisted in the context of our model, when changing the parameters that control the biophysical properties of the neural field (results not shown). Second, the effect is a byproduct of local propagation of activity over the field. In brief, smooth local patterns emerge in densely connected brain regions, the—slow—dynamics of which are shaped by effective connectivity between remote regions. Each local eigenmode is then expected to control the acceleration or slowing down of the response of its respective brain region. In contradistinction to e.g., quantitative predictions about frequency power spectra, we anticipate this qualitative prediction to hold, irrespective of models' specificities.

Let us now consider the (modest) theoretical contribution of this work in relation to the separation of time scales. We believe that—if properly extended—established biophysical models of fast electrophysiological responses may be a good starting point for understanding slow dynamical modes at macroscopic spatial scales that, presumably, drive fMRI responses. The neural field treatment we have developed in this note is a first step toward a mechanistic understanding of the macroscopic dynamics that emerge at the time scale of seconds or minutes. However, this model is limited in many ways. First, we did not account for extrinsic connectivity between regions (coupled neural fields). This is known to strongly impact on the stability of spatiotemporal brain dynamics (Knock et al., 2009), and will thus be critical for predicting the dynamical repertoire of slow macroscopic modes of activity (cf. point above). For example, Bojak et al. (2011) show that increasing the long-range coupling between remote regions strongly modulates the frequency response of the target region, in terms of the relative power in low- and high- frequency bands. Second, we did not separate the respective contributions of excitatory and inhibitory subpopulations (Pinotsis et al., 2012). Of particular importance here is the fact that inhibito-excitatory networks can possess stable limit cycles (Seung et al., 1995). This means that local neural fields can behave as oscillators or resonators. This relates to the work of Cabral et al. (2011) who examined the role of local network oscillations in the emergence of slow temporal synchrony between remote brain areas. The critical thing here is that the time it takes for one resonator to influence another resonator can be much slower than the local resonance frequencies. In addition, owing to the autonomous nature of local dynamics, the effect of between-regions coupling might be better understood in terms of changes in the resonance frequency of the target region. Taken together, these considerations suggest that challenging theoretical work still needs to be undertaken to establish the dynamical repertoire of coupled slow macroscopic modes of brain activity.

We take the fact that the sDCM neural states contribution was significant for all EEG channels (above and beyond slow drifts) as a compelling validation of the conditional estimates of neuronal activity from sDCM. Importantly, this also provides some empirical evidence in favor of the fMRI to EEG mapping itself. This supports the results of Rosa et al. (2010b), which is, to our knowledge, the only experimental study that tested predictions of the heuristic model of Kilner et al. (2005). In this study, the authors experimentally controlled the frequency content of the EEG signal, and showed that the BOLD response in a group subjects was better explained by frequency modulation than by amplitude modulation of the EEG. More specifically, Rosa et al. (2010b) showed that to correctly model BOLD variations it was important to use a normalized frequency spectrum, as in Equation (10). Here, we approach this issue from the other side and use fMRI to predict the EEG frequency spectrum. Note that this is only an indirect validation of the neurovascular coupling model, as we did not compare it with other quantitative scenarios. Among these is the idea that BOLD changes correlate with amplitude modulations within frequency bands of interest in the EEG spectrum. Typically, the BOLD signal is considered a filtered version of the EEG alpha power (see e.g., De Munck et al., 2008). More generally, this idea has been applied to each and every EEG frequency band with relative success [see Laufs et al. (2008) and references therein]. As an illustration, Lachaux et al. (2007) have found that EEG gamma band modulations co-localize with BOLD variations. However, although power fluctuations of different EEG frequency bands are known to be mutually highly correlated, multi frequency models are usually not used in such studies (De Munck et al., 2009). In addition to compromising the specificity of the results, this prevents a direct comparison with our and related work (cf. Rosa et al., 2010b).

Another limitation of our data analysis is that we neglected the unavoidable non-linearity implicit in the relationship $\overline{\omega }=\overline{\omega }\left(z_1^{(0)}\right)$ (cf. Equations 9–12 and Figure 2) when testing for the contribution of slow neural states to modulations of the EEG center frequency. This means we may have not taken full advantage of our neural field formulation. However, this does not detract from our quantitative assessment of the predictive validity of sDCM for fMRI data.

When examining the GSW-related EEG frequency modulations, we noted that the frequency spectrum did not change markedly during GSW activity. This might seem surprising, if one thinks of GSW as sudden bursts of activity within the network, associated with clearly recognizable 3 Hz transient waves of EEG activity (Destexhe and Sejnowski, 2001). Closer inspection of the EEG frequency profile confirmed the appearance of a peak at 3 Hz during GSW-activity (not shown). However, many other (faster) processes actually contribute to the overall frequency spectrum, which typically drags the center frequency toward 20 Hz. This is the main reason why we had to include as many confounding factors as possible, when assessing the specific contribution of the sDCM predictors. Importantly, physiological confounds (cardiac activity, head movements, etc.) span most of the slow frequencies in the EEG signal, which potentially masks interesting slow neural processes. This means that slow frequency EEG signals with a neural origin are confounded with physiological noise. Note that there are notorious exceptions to this detectability issue, for example readiness potentials (Jahanshahi and Hallett, 2003) that are thought to reflect the preparation of motor activity and have a typical time scale of a few seconds. Nevertheless, our results indicate that it may be possible to detect the contribution of slow neural processes to the scalp EEG signal, through their indirect effect on the frequency modulation of fast modes.

In conclusion, we have provided empirical evidence supporting the predictive validity of sDCM for fMRI data. In future work, we will examine the empirical and theoretical relationships that exist between slow and fast dynamical consequences of effective connectivity as can be accessed via fMRI and EEG, respectively.

Conflict of Interest Statement

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

The authors would like to thank Dr. Dimitris Pinotsis for useful discussions. K. J. Friston was funded by the Wellcome Trust. J. Daunizeau was funded by an ERC grant.

Footnotes

1. ^These are O(z²) terms that derive from the second term in the right-hand side of the second line of Equation (A11).

References

Amari, S.-I. (1975). Homogeneous nets of neuron-like elements. Biol. Cybern. 17, 211–220.

Pubmed Abstract | Pubmed Full Text

Amari, S.-I. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. Biol. Cybern. 27, 77–87.

Pubmed Abstract | Pubmed Full Text

Bojak, I., and Liley, D. T. (2010). Axonal velocity distributions in neural field equations. PLoS Comp. Biol. 6:e1000653. doi: 10.1371/journal.pcbi.1000653

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bojak, I., Oostendorf, T. F., Reid, A. T., and Kotter, R. (2011). Towards a model-based integration of co-registered electroencephalography/functional magnetic resonance imaging data with realistic neural population meshes. Philos. Transact. A Math. Phys. Eng. Sci. 369, 3785–3801.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bressloff, P. C. (2012). Spatiotemporal dynamics of continuum neural fields. J. Phys. A Math. Theor. 45:033001. doi: 10.1088/1751-8113/45/3/033001

CrossRef Full Text

Brodersen, K. H., Schofield, T. M., Leff, A. P., Ong, C. S., Lomakina, E. I., Buhmann, J. M., et al. (2011). Generative embedding for model-based classification of fMRI data. PLoS Comput. Biol. 7:e1002079. doi: 10.1371/journal.pcbi.1002079

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Buxton, R. B., Wong, E. C., and Franck, L. R. (1998). Dynamics of blood flow and oxygenation changes during brain activation: the Balloon model. Magn. Res. Med. 39, 855–864.

Pubmed Abstract | Pubmed Full Text

Cabral, J., Hugues, E., Sporns, O., and Deco, G. (2011). Role of local network oscillations in resting-state functional connectivity. Neuroimage 57, 130–139.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Coombes, S. (2010). Large-scale neural dynamics: simple and complex. Neuroimage 52, 731–739.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Coombes, S., Venkov, N. A., Shiau, L., Bojak, I., Liley, D. T. J., and Laing, C. R. (2007). Modeling electrocortical activity through local approximations of integral neural field equations. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 76, 1539–1547.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Daunizeau, J., David, O., and Stephan, K. E. (2011). Dynamic Causal Modelling: a critical review of the biophysical and statistical foundations. Neuroimage 58, 312–322.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Daunizeau, J., Kiebel, S. J., and Friston, K. J. (2009a). Dynamic causal modeling of distributed electromagnetic responses. Neuroimage 47, 590–601.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Daunizeau, J., Kiebel, S. J., and Friston, K. J. (2009b). Variational Bayesian identification and prediction of stochastic dynamic causal models. Physica D 289, 2089–2118.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Daunizeau, J., Laufs, H., and Friston, K. J. (2010). “EEG-fMRI information fusion: biophysics and data analysis,” in EEG-fMRI- Physiology, Technique and Applications, eds C. Mulert and L. Lemieux (Berlin, Heidelberg: Springer), 511–526. ISBN: 978-3-540-87918-3.

Daunizeau, J., Stephan, K. E., and Friston, K. J. (2012). Stochastic dynamic causal modelling of fMRI data: should we care about neural noise? Neuroimage 62, 464–481.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

David, O., Guillemain, I., Baillet, S., Reyt, S., Deransart, C., Segebarth, C., et al. (2008). Identifying neural drivers with functional MRI: an electrophysiological validation. PLoS Biol. 6:e315. doi: 10.1371/journal.pbio.0060315

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

De Munck, J. C. (1988). The potential distribution in a layered spheroidal volume conductor. J. Appl. Phys. 64, 464–470.

De Munck, J. C., Gonçalves, S. I., Faes, T. J., Kuijer, J. P., Pouwels, P. J., Heethaar, R. M., et al. (2008). A study of the brain's resting state based on alpha band power, heart rate and fMRI. Neuroimage 42, 112–121.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

De Munck, J. C., Goncalves, S. I., Mammoliti, R., Heethaar, R. M., and Lopes da Silva, F. H. (2009). Interactions between different EEG frequency bands and their effect on alpha-fMRI correlations. Neuroimage 47, 69–76.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Destexhe, A., and Sejnowski, T. J. (2001). Thalamocortical Assemblies. Oxford, UK: Oxford University Press.

Freeman, W. J. (1975). Mass Action in the Nervous System: Examination of the Neurophysiological Basis of Adaptive Behavior through the EEG, 1st Edn. New York, NY: Academic Press.

Friston, K. J., Harrison, L., and Penny, W. D. (2003). Dynamic causal modelling. Neuroimage 19, 1273–1302.

Pubmed Abstract | Pubmed Full Text

Friston, K. J., Li, B., Daunizeau, J., and Stephan, K. E. (2011). Network discovery with DCM. Neuroimage 56, 1202–1221.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Friston, K. J., Mattout, J., Trujillo-Barreto, N., Ashburner, J., and Peeny, W. (2007). Variational free energy and the Laplace approximation. Neuroimage 34, 220–234.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Friston, K. J., and Penny, W. (2011). Post-hoc Bayesian model comparison. Neuroimage 56, 2089–2099.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Friston, K. J., Stephan, K. E., Li, B., and Daunizeau, J. (2010). Generalised filtering. Math. Probl. Eng. 2010:621670. doi: 10.1155/2010/621670

CrossRef Full Text

Friston, K. J., Trujillo-Barreto, N. J., and Daunizeau, J. (2008). DEM: a variational treatment of dynamical systems. Neuroimage 42, 849–885.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Garrido, M. I., Kilner, J., Kiebel, S. J., and Friston, K. J. (2007). Evoked brain responses are generated by feedback loops. Proc. Natl. Acad. Sci. U.S.A. 104, 20961–20966.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Haken, H. (1983). Synergetics, an Introduction: Nonequilibrium Phase Transitions and Self-Organization in Physics, Chemistry, and Biology, 3rd Rev. enl. ed. New York, NY: Springer-Verlag.

Hamandi, K., Salek-Haddadi, A., Laufs, H., Liston, A., Friston, K., Fish, D. R., et al. (2006). EEG-fMRI of idiopathic and secondarily generalized epilepsy. Neuroimage 31, 1700–1710.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hesselmann, G., Sadaghiani, S., Friston, K. J., and Kleinschmidt, A. (2010). Predictive coding or evidence accumulation? False inference and neuronal fluctuations. PLoS ONE 5:e9926. doi: 10.1371/journal.pone.0009926

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jahanshahi, M., and Hallett, M. (eds.). (2003). The Bereitschaftspotential, Movement-Related Cortical Potentials. Springer Series. New York, NY: Kluwer Academic / Plenum Publishers. ISBN: 0306474077.

Jansen, B. H., and Rit, V. G. (1995). Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol. Cybern. 73, 357–366.

Pubmed Abstract | Pubmed Full Text

Jirsa, V., and Haken, H. (1996). Field theory of electromagnetic brain activity. Phys. Rev. Lett. 77, 960–963.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kass, R. E., and Raftery, A. E. (1995). Bayes factors. J. Am. Stat. Assoc. 90, 773–795.

Kilner, J. M., Mattout, J., Henson, R., and Friston, K. J. (2005). Hemodynamic correlates of EEG: a heuristic. Neuroimage 28, 280–286.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Knock, S. A., McIntosh, A. R., Sporns, O., Kötter, R., Hagmann, P., and Jirsa, V. K. (2009). The effects of physiologically plausible connectivity structure on local and global dynamics in large scale brain models. J. Neurosci. Methods 183, 86–94.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Krakow, K., Allen, P. J., Lemieux, L., and Fish, D. R. (2000). Methodology: EEG-correlated fMRI. Adv. Neurol. 83, 187–201.

Pubmed Abstract | Pubmed Full Text

Lachaux, J. P., Fonlupt, P., Kahane, P., Minotti, L., Ho mann, D., Bertrand, O., et al. (2007). Relationship between task-related gamma oscillations and BOLD signal: new insights from combined fMRI and intracranial EEG. Hum. Brain Mapp. 28, 1368–1375.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Laufs, H., Daunizeau, J., Carmichael, D. W., and Kleinschmidt, A. (2008). Recent advances in recording electrophysiological data simultaneously with magnetic resonance imaging. Neuroimage 40, 515–528.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lemieux, L., Salek-Haddadi, A., Lund, T. E., Laufs, H., and Carmichael, D. (2007). Modelling large motion events in fMRI studies of patients with epilepsy. Magn. Reson. Imaging 25, 894–901.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Li, B., Daunizeau, J., Stephan, K. E., Penny, W., and Friston, K. (2011). Stochastic DCM and generalized filtering. Neuroimage 58, 442–457.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liley, D. T. J., Cadusch, P. J., and Dafilis, M. P. (2002). A spatially continuous mean field theory of electrocortical activity. Network 13, 67–113.

Pubmed Abstract | Pubmed Full Text

Liley, D. T. J., Foster, B. L., and Bojak, I. (2012). “Cooperative populations of neurons: mean field models of mesoscopic brain activity,” in Computational Systems Neurobiology, ed N. Le Novère (Netherlands: Springer), 315–362.

Liston, A. D., Lund, T. E., Salek-Haddadi, A., Hamandi, K., Friston, K. J., and Lemieux, L. (2006). Modelling cardiac signal as a confound in EEG-fMRI and its application in focal epilepsy study. Neuroimage 30, 827–834.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Logothetis, N. K., and Wandell, B. A. (2004). Interpreting the BOLD signal. Annu. Rev. Physiol. 66, 735–769.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lopes da Silva, F. H., Hoeks, A., Smits, H., and Zetterberg, L. H. (1973). Model of brain rhythmic activity: the alpha-rhythm of the thalamus. Kybernetik 15, 27–37.

Pubmed Abstract | Pubmed Full Text

Marreiros, A. C., Kiebel, S. J., Daunizeau, J., Harrison, L. M., and Friston, K. J. (2008). Population dynamics under the Laplace assumption. Neuroimage 44, 701–714.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Moran, R. J., Jung, F., Kumagai, T., Endepols, H., Graf, R., Dolan, R. J., et al. (2011). Dynamic causal models and physiological inference: a validation study using isoflurane anaesthesia in rodents. PLoS ONE 6:e22790. doi: 10.1371/journal.pone.0022790

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nunez, P. L. (1974). The brain wave equation: a model for the EEG. Math. Biosci. 21, 279–297.

Penny, W., Joao, M., Flandin, G., Daunizeau, J., Stephan, K. E., Friston, K. J., et al. (2010). Comparing families of dynamic causal models. PLoS Comp. Biol. 6:e1000709. doi: 10.1371/journal.pcbi.1000709

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pinotsis, D. A., Moran, R., and Friston, K. J. (2012). Dynamic causal modelling with neural fields. Neuroimage 59, 1261–1274.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Roebroeck, A., Formisano, E., and Goebel, R. (2011). The identification of interacting networks in the brain using fMRI: model selection, causality and deconvolution. Neuroimage 58, 296–302.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rosa, M. J., Daunizeau, J., and Friston, K. J. (2010a). EEG-fMRI integration: a critical review of biophysical modelling and data analysis approaches. J. Integr. Neurosci. 9, 453–476.

Pubmed Abstract | Pubmed Full Text

Rosa, M. J., Kilner, J., Blankenburg, F., Josephs, O., and Penny, W. (2010b). Estimating the transfer function from neuronal activity to BOLD using simultaneous EEG-fMRI. Neuroimage 49, 1496–1509.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Schmuel, A. (2010). “Locally measured neuronal correlates of functional MRI signals,” in EEG-fMRI- Physiology, Technique and Applications, eds C. Mulert and L. Lemieux (Berlin, Heidelberg: Springer), 63–82. ISBN: 978-3-540-87918-3

Schultze-Kraft, M., Becker, R., Breakspear, M., and Ritter, P. (2011). Exploiting the potential of three dimensional spatial wavelet analysis to explore nesting of temporal oscillations and spatial variance in simultaneous EEG-fMRI data. Progr. Biophys. Mol. Biol. 105, 67–79.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Seung, H. S., Richardson, T. J., Lagarias, J. C., and Hopfield, J. J. (1995). “Minimax and Hamiltonian dynamics of excitatory-inhibitory networks,” in NIPS '97 Proceedings of the 1997 Conference on Advances in Neural Information Processing Systems, Vol. 10 (Cambridge, MA: MIT Press), 329–335.

Smith, S. M., Miller, K. L., Salimi-Khorshidi, G., Webster, M., Beckmann, C. F., Nichols, T. E., et al. (2011). Network modelling methods for fMRI. Neuroimage 54, 875–891.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sporns, O. (2010). Networks of the Brain. Cambridge, MA: MIT Press. ISBN: 9780262014694.

Srinivasan, R., Tucker, D. M., and Murias, M. (1998). Estimating the spatial Nyquist of the human EEG. Behav. Res. Methods Instrum. Comput. 30, 8–19.`

Stephan, K. E., Kasper, L., Harrison, L., Daunizeau, J., Den Ouden, H., Breakspear, M., et al. (2008). Nonlinear dynamic causal models for fMRI. Neuroimage 42, 649–662.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stephan, K. E., Penny, W. D., Daunizeau, J., Moran, R. J., and Friston, K. J. (2009). Bayesian model selection for group studies. Neuroimage 46, 1004–1017.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stephan, K. E., Weiskopf, N., Drysdale, P. M., Robinson, P. A., and Friston, K. J. (2007). Comparing hemodynamic models with DCM. Neuroimage 38, 387–401.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tononi, G., Sporns, O., and Edelman, G. M. (1994). A measure for brain complexity: relating functional segregation and integration in the nervous system. Proc. Natl. Acad. Sci. U.S.A. 91, 5033–5037.

Pubmed Abstract | Pubmed Full Text

Valdes-Sosa, P., Roebroeck, A., Daunizeau, J., and Friston, K. (2011). Effective connectivity: influence, causality and biophysical modelling. Neuroimage 58, 339–361.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Vaudano, A. E., Laufs, H., Kiebel, S. J., Carmichael, D. W., Hamandi, K., Guye, M., et al. (2009). Causal hierarchy within the thalamo-cortical network in spike and wave discharges. PLoS ONE 4:e6475. doi: 10.1371/journal.pone.0006475

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Von Holdt, R. E. (1965). Rational power of power series. Am. Math. Mon. 72, 740–743.

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neuron. Biophys. J. 12, 1–24.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wilson, H. R., and Cowan, J. D. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik 13, 55–80.

Pubmed Abstract | Pubmed Full Text

Appendix 1: Approximating the Local Neural Fields with a Spatial Eigenmodes Expansion

Let ∇² = ∂²/∂x² +∂²/∂y² be the 2D-Laplacian operator defined on a bounded (square) 2D domain of size $\Vert D\Vert$, where {x, y} = r are the two cartesian coordinates on the domain. Recall that eigenvectors w^{(k, l)} (r) of ∇² are such that:

$\begin{array}{cc}\begin{array}{l}{w}^{(k,l)}(r)=e^{i\frac{\pi }{\Vert D\Vert }(kx+ly)}\\ =w^{(k)}(x)w^{(l)}(y)\\ ={\left[w^{(1)}(x)\right]}^k{\left[w^{(1)}(y)\right]}^l\end{array}& \text{(A1)}\end{array}$

where we have used index pairs (k, l) for spatial frequencies along the x- and y-axis, respectively. The associated eigenvalues are given by: ${\lambda }_{k,l}=-\left(k^2+l^2\right){\pi }^2/\Vert D{\Vert }^2$. By convention, w^(k) refer to the eigenfunctions of the 1D-Laplacian, and its complex conjugate will be denoted as w^(−k). Now, projecting Equation (5) onto the (complex conjugate) spatial eigenmodes yields the eigendynamics of the neural field (Daunizeau et al., 2009a):

$\begin{array}{cc}\begin{array}{l}{\tau }^2{\ddot{z}}_1^{(k,l)}(t)=-2\tau {\dot{z}}_1^{(k,l)}(t)-z_1^{(k,l)}(t)+\tau \gamma z_2^{(k,l)}(t)\\ {\dot{z}}_2^{(k,l)}(t)=\frac{\overline{\nu }}{2\sigma }\left(\left({\sigma }^2{\lambda }_{kl}-1\right)z_2^{(k,l)}(t)+G^0\underset{{\varsigma }_{k,l}(t)}{\underbrace{{\displaystyle {\int }_D{w}^{(-k,-l)}}(r)S\left(\eta (r,t)\right)dr}}\right).\end{array}& (\text{A}2)\end{array}$

where the projections ς_{k, l} (t) of the firing rate form a complete basis function set.

Together with a truncated eigenmode expansion, Equation (A2) can be used to approximate the solution to Equation (5) of the main text. Such a truncation basically acts as a low-pass filter, neglecting the transient effects of higher spatial frequency modes. This can be motivated by noting that (1) the eigenmodes w^{(k, l)} (r) of the Laplacian operator are of increasing spatial frequency as the order (k, l) increases and (2) the associated eigenvalues λ_kl are negative (except for λ₀₀ = 0) with increasing magnitude. This implies that the dynamics of high-order spatial modes will be strongly damped and thus quickly disappear in the absence of external forcing.

We now focus on deriving an approximation to the projections ς_{k, l} (t) of the firing rate.

First, let us Taylor-expand the sigmoidal mapping of PSP activity around zero:

$\begin{array}{cc}S(\eta (r,t))={\displaystyle \sum _{m=0}^{\infty }\frac{1}{m!}}{S}_0^{[m]}{\eta }^m(r,t),& \text{(A3)}\end{array}$

where $S_0^{[m]}={\frac{{\partial }^{m}S}{\partial x^m}|}_0$ is the mth derivative of the sigmoidal mapping evaluated at η = 0 and η^m (r, t) is the mth power of mean PSP η(r, t). Assuming that Equation (A1) holds for the eigenvectors of the Laplacian operator defined onto the cortical domain D, Equation (6) becomes:

$\begin{array}{cc}\begin{array}{l}\eta (r,t)={\displaystyle \sum _k{\displaystyle \sum _l{z}_1^{(k,l)}}}(t)w^{(k,1)}(r)\\ ={\displaystyle \sum _k{\tilde{z}}_1^{(k)}}(t,y){\left[w^{(1)}(x)\right]}^k\\ {\tilde{z}}_1^{(k)}(t,y)\equiv {\displaystyle \sum _l{z}_1^{(k,l)}}(t){\left[w^{(1)}(y)\right]}^l\end{array}& \text{(A4)}\end{array}$

where ${\tilde{z}}_1^{(k)}(t,y)$ is a dummy eigenmode dynamics, that is still a function of space (more precisely, of the second coordinate y). Inserting Equation (A4) into (A3) yields the following power series of power series in w⁽¹⁾ :

$\begin{array}{cc}\begin{array}{l}S(\eta (r,t))={\displaystyle \sum _{m=0}^{\infty }\frac{1}{m!}}{S}_0^{[m]}{\left({\displaystyle \sum _k{\tilde{z}}_1^k}(t,y){\left[w^{(1)}(x)\right]}^k\right)}^m\\ ={\displaystyle \sum _{m=0}^{\infty }\frac{1}{m!}}{S}_0^{[m]}{\displaystyle \sum _k{A}_k}(t,y,m)w^{(k)}(x)\end{array}& \text{(A5)}\end{array}$

where the coefficients A_k (t, y, m) obey the following recursive equation (Von Holdt, 1965):

$\begin{array}{cc}\{\begin{array}{l}{A}_k(t,y,m)=\frac{1}{k{\tilde{z}}_1^{(0)}(t)}{\displaystyle {\sum }_{i=1}^k((m+1)i-k){\tilde{z}}_1^{(i)}(t)A_{k-i}(t,y,m)}\\ A_0(t,y,m)={\left[{\tilde{z}}_1^{(0)}(t)\right]}^m\end{array}& \text{(A6)}\end{array}$

The contribution ϒ_k(t, y) of the kth spatial mode to the expansion in Equation (A5) is given by:

$\begin{array}{cc}\begin{array}{l}S(\eta (r,t))={\displaystyle \sum _k{\Upsilon }_k}(t,y)w^{(k)}(x)\\ {\Upsilon }_k(t,y)={\displaystyle \sum _{m=0}^{\infty }\frac{1}{m!}}{S}_0{}^{[m]}{A}_k(t,y,m)\end{array}& \text{(A7)}\end{array}$

These can be analytically derived from the recursive relationship given in Equation (A6). They are given in Table A1 bellow, up to order k = 3.

TABLE A1

This table gives both the coefficients A_k (t, y, m) (left) and the contribution of the kth spatial mode to the expansion in Equation (A4), as given by ϒ_k(t, y) (right) up to k = 3.

It can be seen form Table A1 that, except for the fundamental mode, the spatial mode contribution ϒ_k is of the form:

$\begin{array}{cc}\begin{array}{l}{\Upsilon }_k\propto {\displaystyle \sum _{j=1}^k\text{O}}\left({\left[{\tilde{z}}_1\right]}^j\right)S^{[j]}\left({\tilde{z}}_1^{(0)}\right)\\ ={\displaystyle \sum _{j=1}^k\text{O}}\left({\left[{\tilde{z}}_1\right]}^j\right){\rho }^{j}S\left({\tilde{z}}_1^{(0)}\right){\displaystyle \prod _{m=1}^j\left(1-mS\left({\tilde{z}}_1^{(0)}\right)\right)},\\ \approx {\displaystyle \sum _{j=1}^k\text{O}}\left({\left[{\tilde{z}}_1\right]}^{2j+1}\right)\end{array}& \text{(A8)}\end{array}$

where O(x^j) is a polynomial of order j in x and ρ is the sigmoidal slope. The second line in Equation (A8) follows from the following property of the sigmoid function: S^[j] (x) = ρ^j S(x) ∏^j_{m = 1}(1 − mS(x)) and the last line would follow from a first-order Taylor expansion of the sigmoidal mapping. This means that the lower order term in each spatial mode contribution ϒ_k(t, y) is already a polynomial of order three in the dummy eigenmodes dynamics ${\tilde{z}}_1^{(k)}(t,y)$. We will therefore truncate the contribution to its first term, i.e.,: ${\Upsilon }_0(t,y)=S\left({\tilde{z}}_1^{(0)}(t,y)\right)$ for k = 0, and ${\Upsilon }_k(t,y)\approx {\tilde{z}}_1^{(k)}(t,y)S^{\prime }\left({\tilde{z}}_1^{(0)}(t,y)\right)$ for k > 0. This yields the following approximation to the firing rate field:

$\begin{array}{cc}S(\eta (r,t))\approx S\left({\tilde{z}}_1^{(0)}(t,y)\right)w^{(0)}(x)+{\displaystyle \sum _{k\ge 0}\rho }{S}^{\prime }\left({\tilde{z}}_1^{(0)}(t,y)\right){\tilde{z}}_1^{(k)}(t,y)w^{(k)}(x)& \text{(A9)}\end{array}$

Equation (A9) can now be used to approximate the integral with respect to x in Equation (A2):

$\begin{array}{cc}\begin{array}{l}{\varsigma }_{k,l}(t)={\displaystyle \int w^{(-l)}}(y)\left[{\displaystyle \int w^{(-k)}}(x)S(\eta (r,t))dx\right]\text{}dy\\ \approx {\delta }_k{\displaystyle \int w^{(-l)}}(y)S\left({\tilde{z}}_1^{(0)}(t,y)\right)dy+(1-{\delta }_k)\rho {\displaystyle \int w^{(-l)}}(y)S^{\prime }({\tilde{z}}_1^{(0)}(t,y)){\tilde{z}}_1^{(k)}(t,y)dy\\ ={\delta }_k{\displaystyle \int w^{(-l)}}(y)S\left({\tilde{z}}_1^{(0)}(t,y)\right)dy+(1-{\delta }_k)\rho {\displaystyle \sum _{l^{\prime }}{z}_1^{k,l}}(t){\displaystyle \int w^{(l^{\prime }-l)}}(y)S^{\prime }\left({\tilde{z}}_1^{(0)}(t,y)\right)\text{}dy\end{array}& \text{(A10)}\end{array}$

where δ_k is a Dirac delta, which is zero except when k = 0.

Similarly, the sigmoidal terms in the integrands of Equation (A10) can be approximated using power series of power series, as follows:

$\begin{array}{cc}\begin{array}{l}S\left({\tilde{z}}_1^{(0)}(t,y)\right)\approx S\left(z_1^{(0,0)}(t)\right)w^{(0)}(y)+{\displaystyle \sum _{l\ge 0}\rho }{S}^{\prime }\left(z_1^{(0,0)}(t)\right)z_1^{(0,l)}(t)w^{(l)}(y)\\ S^{\prime }\left({\tilde{z}}_1^{(0)}(t,y)\right)\approx S^{\prime }\left(z_1^{(0,0)}(t)\right)w^{(0)}(y)+{\displaystyle \sum _{l\ge 0}\rho }{S}^{″}\left(z_1^{(0,0)}(t)\right)z_1^{(0,l)}(t)w^{(l)}(y)\end{array}& \text{(A11)}\end{array}$

Again, inserting Equation (A11) into (A10), and ignoring high-order terms¹ yields:

$\begin{array}{cc}{\varsigma }_{k, l}(t)={\delta }_k{\delta }_{l}S\left(z_1^{(0,0)}(t)\right)+\left(1-{\delta }_k{\delta }_l\right)\rho S^{\prime }\left(z_1^{(0,0)}(t)\right)z_1^{(k,l)}(t)& \text{(A12)}\end{array}$

Inserting Equation (A12) in (A2) yields:

$\begin{array}{cc}{\dot{z}}_2^{(k,l)}(t)=\{\begin{array}{l}\frac{\overline{\nu }}{2\sigma }\left(G^{0}S\left(z_1^{(0,0)}(t)\right)-z_2^{(0,0)}(t)\right)\text{ if} k=l=0\\ \frac{\overline{\nu }}{2\sigma }\left(\left({\sigma }^2{\lambda }_{kl}-1\right)z_2^{(k,l)}(t)+\rho G^0{S}^{\prime }\left(z_1^{(0,0)}(t)\right)z_1^{(k,l)}(t)\right)\text{ otherwise}\end{array}& \text{(A13)}\end{array}$

The second line of Equation (A13) is Equation (8) of the main text, having re-ordered the eigenmodes index pairs (k, l) → k, for the sake of notational simplicity.

The critical properties of this derivation can be summarized as follows:

First, the fundamental eigenmode response will be slower than higher-order eigenmodes (because of smaller dampening terms). Second, the response of fast eigenmodes is enslaved by the slow eigenmode. Third, the fast eigenmodes do not feedback onto the slow eigenmodes. Taken together, we can think of the neural field as a feedforward system: inputs to the field activate the slow eigenmode z^{(0, 0)}₁ (t), which controls the fast modes z^{(k, l)}₁ (t). This is further discussed in the main text of this document.

Appendix 2: Deriving the Frequency Spectrum of Eigenmodes

Recall that the eigenmode dynamics obey the following ODE (cf. Equation 8 in the main text):

$\begin{array}{cc}\begin{array}{l}{\dot{z}}^{(k)}=A{z}^{(k)}+\nu \hfill \\ A=\left[\begin{array}{ccc}0& 0& 1\\ c(z_1^{(0)})& -\frac{\overline{\nu }}{2\sigma }(1-{\sigma }^2{\lambda }_k)& 0\\ -\frac{1}{{\tau }^2}& \frac{\gamma }{\tau }& -\frac{2}{\tau }\end{array}\right]\hfill \end{array}& (\text{A}14)\end{array}$

where ν = O(z²) are high-order error terms that act as a perturbation to the deterministic evolution of eigenmodes. As a simplifying assumption, we consider that ν behaves as i.i.d. white noise, the Fourier transform F[z^(k)] of Equation (A8) writes:

$\begin{array}{cc}\begin{array}{l}F\left[z^{(k)}\right]={\left(2i\pi \omega I_3-A\right)}^{-1}F[\nu ]\\ \propto {\left(2i\pi \omega I_3-A\right)}^{-1}{1}_3\end{array}& \text{(A15)}\end{array}$

where 1₃ is a 3 × 1 vector of ones, and I₃ is the 3 × 3 identity matrix.

Note that the inverse translated Jacobian writes:

$\begin{array}{cc}\begin{array}{l}{\left(2i\pi \omega I_3-A\right)}^{-1}=\frac{1}{\left|2i\pi \omega I_3-A\right|}\left[\begin{array}{ccc}{\xi }_k(\omega )\left(i\omega +\frac{2}{\tau }\right)& \frac{\gamma }{\tau }& {\xi }_k(\omega )\\ c\left(z_1^{(0)}\right)\left(i\omega +\frac{2}{\tau }\right)& i\omega \left(i\omega +\frac{2}{\tau }\right)+\frac{1}{{\tau }^2}& c\left(z_1^{(0)}\right)\\ \frac{\gamma }{\tau }c\left(z_1^{(0)}\right)-{\xi }_k(\omega )\frac{1}{{\tau }^2}& i\omega \frac{\gamma }{\tau }& i\omega {\xi }_k(\omega )\end{array}\right]\\ {\xi }_k(\omega )=i\omega +\frac{\overline{\nu }}{2\sigma }\left(1-{\sigma }^2{\lambda }_k\right)\end{array}& \text{(A16)}\end{array}$

where ξ_k (ω) is the only term that depends upon the mode order k, and the determinant of the translated Jacobian is given by:

$\begin{array}{cc}\begin{array}{l}\left|2i\pi \omega I_3-A\right|=i\omega {\xi }_k(\omega )\left(i\omega +\frac{2}{\tau }\right)-\frac{\gamma }{\tau }c\left(z_1^{(0)}\right)+{\xi }_k(\omega )\frac{1}{{\tau }^2}\\ ={\xi }_k(\omega )\left(i\omega \left(i\omega +\frac{2}{\tau }\right)+\frac{1}{{\tau }^2}\right)-\frac{\gamma }{\tau }c\left(z_1^{(0)}\right)\end{array}& \text{(A17)}\end{array}$

where ξ_k (ω) is given in Equation (A12) (with appropriate eigenmodes index reordering (k, l) → k). The explicit Fourier transform of eigenmode dynamics now derives from inserting Equation (A16) into (A15):

$\begin{array}{cc}F\text{}\left[z^{(k)}\right]\propto \frac{1}{\left|2i\pi \omega I_3-A\right|}\left[\begin{array}{c}{\xi }_k(\omega )\left(i\omega +\frac{2}{\tau }+1\right)+\frac{\gamma }{\tau }\\ \left(i\omega +\frac{2}{\tau }\right)\left(c\left(z_1^{(0)}\right)+i\omega \right)+\frac{1}{{\tau }^2}+c\left(z_1^{(0)}\right)\\ {\xi }_k(\omega )\left(i\omega -\frac{1}{{\tau }^2}\right)+\frac{\gamma }{\tau }\left(i\omega +c\left(z_1^{(0)}\right)\right)\end{array}\right]& (\text{A}18)\end{array}$

Now, we are interested in the total frequency power of EEG signal sampled at some position r' on the scalp that is induced by mean membrane depolarization:

$\begin{array}{cc}\begin{array}{l}\text{P}(\omega ,r^{\prime })\equiv {\left|{\displaystyle {\int }_{D}L}(r,r^{\prime })F[\eta (r,t)]dr\right|}^2\\ ={\left|{\displaystyle \sum _k{\alpha }^{(k)}}(r^{\prime })F\left[z_1^{(k)}\right]\right|}^2\\ {\alpha }^{(k)}(r^{\prime })={\displaystyle {\int }_{D}L}(r,r^{\prime })w^{(k)}(r)dr\end{array}& (\text{A}19)\end{array}$

where the second line derives from the eigendecomposition of the neural field (cf. Equation 6), L(r, r′) is a spatial convolution operator that describes the quasi-static volume propagation of electromagnetic fields through head tissues, and α^(k) (r′) captures the contribution of the kth eigenmode at position r′ on the scalp. It is known that the geometrical properties of electrical conductivity within head tissues eventually shape the EEG susceptibility to different cortical spatial scales. For the sake of simplicity, we will assume the skull acts as a perfect low-pass filter (see e.g., Srinivasan et al., 1998), which induces a cut-off order k_c, i.e.,:

$\begin{array}{cc}{\alpha }^{(k)}({\text{r}}^{\prime })=\{\begin{array}{cc}\begin{array}{l}{\alpha }^{(0)}(r^{\prime })\\ 0\end{array}& \begin{array}{l}\text{if\hspace{0.17em}}k\le k_c\\ \text{if\hspace{0.17em}}k>k_c\end{array}\end{array}& (\text{A}20)\end{array}$

Equation (A20) basically neglects the continuous loss of sensitivity as one increases the spatial frequency of (cortical) membrane depolarization fields. Inserting Equation (A20) into (A19) yields Equation (9) of the main text, i.e.,:

$\begin{array}{cc}\begin{array}{l}\text{P}\left(\omega ,r^{\prime }\right)={\alpha }^{(0)}{\left(r^{\prime }\right)}^2{\left|{\displaystyle \sum _{k\le k_c}F}\left[z_1^{(k)}\right]\right|}^2\\ ={\alpha }^{(0)}{\left(r^{\prime }\right)}^2{\left|{\displaystyle \sum _{k\le k_c}\frac{{\xi }_k(\omega )\left(i\omega +\frac{2}{\tau }+1\right)+\frac{\gamma }{\tau }}{{\xi }_k(\omega )\left(i\omega \left(i\omega +\frac{2}{\tau }\right)+{\frac{1}{\tau }}^2\right)-\frac{\gamma }{\tau }c\left(z_1^{(0)}\right)}}\right|}^2\end{array}& (\text{A}21)\end{array}$

where α⁽⁰⁾ (r′) is a scaling factor that will be different for different EEG sensors.

Equation (A21) is a valid approximation of the EEG power spectrum at position r′ on the scalp generated by one local brain region. The generalization to a network of remote brain regions requires the specification of the relative contribution of each region, which depends on their relative distance with EEG scalp sensors. In principle, this can be derived from knowledge of the forward operator L(r, r′) (see e.g., De Munck, 1988). This is beyond the scope of the present study.

Lastly, note that the numerical values of the Laplacian's eigenvalues are approximated from the 2D-Euclidean case, as follows: ${\lambda }_{k,l}=-\left(k^2+l^2\right){\pi }^2/\Vert D{\Vert }^2$. The effective frequency cut-off k_c is such that k² + l² ≤ k²_c. In practice, we use k²_c = 16 and $\Vert D\Vert =1$Cm.