An electrophysiological validation of stochastic DCM for fMRI

In this note, we assess the predictive validity of stochastic dynamic causal modeling (sDCM) of functional magnetic resonance imaging (fMRI) data, in terms of its ability to explain changes in the frequency spectrum of concurrently acquired electroencephalography (EEG) signal. We first revisit the heuristic model proposed in Kilner et al. (2005), which suggests that fMRI activation is associated with a frequency modulation of the EEG signal (rather than an amplitude modulation within frequency bands). We propose a quantitative derivation of the underlying idea, based upon a neural field formulation of cortical activity. In brief, dense lateral connections induce a separation of time scales, whereby fast (and high spatial frequency) modes are enslaved by slow (low spatial frequency) modes. This slaving effect is such that the frequency spectrum of fast modes (which dominate EEG signals) is controlled by the amplitude of slow modes (which dominate fMRI signals). We then use conjoint empirical EEG-fMRI data—acquired in epilepsy patients—to demonstrate the electrophysiological underpinning of neural fluctuations inferred from sDCM 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 macroscale 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., 2007or 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., 2008or 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 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,  showed how sDCM can be used in combination with post-hoc model comparison  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 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 Cowan, 1972, 1973;Lopes da Silva et al., 1973;Nunez, 1974;Amari, 1975Amari, , 1977Freeman, 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;. 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  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): 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 2 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 Frontiers in Computational Neuroscience www.frontiersin.org January 2013 | Volume 6 | Article 103 | 2 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: where J is the synaptic impulse response function, τ is a timeconstant, 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: 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: where the Heaviside step function endows Equation (3) with temporal causality and G 0 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 (ν 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).
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): 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 presynaptic firing is accumulated locally to produce depolarization.  (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). 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.,: ∇ 2 w (k) (r) = λ k w (k) (r), where the distribution of eigenvalues λ is such that: λ k < λ 0 = 0. They form a complete and orthonormal basis function set for the neural field, allowing for an eigen-decomposition of the field, as follows: 1 is a non-linear function of the fundamental mode z (0) 1 that mediates the coupling between modes describing firing rate and depolarization. Roughly speaking, each mode decays at a rate that is proportional to σνλ 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ν becomes small compared to the dispersion of lateral connections σ (i.e., ifν/σ → 0). In this case, the fundamental mode z (0) 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 (0) 1 . 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 (0) 1 within a given brain region, over which the neural field is defined. This means that slow fluctuations in the fundamental mode z (0) 1 should be associated with changes in the EEG frequency spectrum P(ω), which derives from the Laplace transform of Equation (8) (see Appendix 2): 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 1 . Similarly to Kilner et al. (2005), let us now define the EEG center frequencyω as the first moment of the (normalized) frequency power spectrumP(ω): By construction, normalized frequency powerP(ω) 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ω 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ω over time quantify frequency modulations of the EEG signal. The key idea here is that the EEG center frequencyω is a function of the slow eigenmode dynamics z (0) 1 , and its modulation over time thus follows the slow eigenmode dynamics z (0) 1 , i.e.,ω =ω z (0) 1 . Figure 2 depicts the frequency modulation of the EEG signal induced by changes in the slow eigenmode z (0) 1 . 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 (0) 1 → χ), 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 powerP(ω). As a consequence, the center frequencyω 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ω (as well as its susceptibility to the eigenmode z (0) 1 ) 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.

THE LINK BETWEEN EEG TIME-FREQUENCY RESPONSE AND fMRI TEMPORAL RESPONSE
In this paper, we assume that the slow macroscopic dynamics x (n) underlying fMRI data in a given region of interest follow the fundamental eigenmode z (0) 1 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 highdimensional dynamical systems, which has recently been used to motivate the rate of change of neural states in stochastic DCM for fMRI .
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 (0) 1 . 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 (0) 1 . 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):ẋ where θ (n) = {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 (n) (cf. Stephan et al., 2008). The influence of the fast modes on the motion of slow modes x (n) 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 (n) 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ω, 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:ω where the gradient ∂ω/∂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ω(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 Vaudano et al., 2009), whose acquisition protocol and data pre-processing are briefly summarized here.
Ten-channel EEG (10-20 system) was recorded using MRcompatible 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 2 , 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 GSWrelated 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 interscan 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 GSWrelated 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: here, x (n) 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 (n) and as a function of hemodynamic parameters θ (h) = κ s , κ f , α, τ 0 , E 0 , ε 0 (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): where ε is an additive measurement noise, and the parameters are defined in Table 2 below. 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 Frontiers in Computational Neuroscience www.frontiersin.org January 2013 | Volume 6 | Article 103 | 7 free-energy lower bound F(q) to the model evidence, with respect to an approximate conditional density q( ): 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 (n) 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 x = x (n) q of the conditional density q(x). In the following,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 1 , whose dependent variable was the trajectoryω 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): wherex are the sDCM conditional estimates, X 0 contains the confounds (see below), β and β 0 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 ∂ω/∂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 β 0 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 1 , relative to a reduced model H 0 that only included the confounds: 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 dynamicsx and the EEG center frequency trajectoryω, 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 0 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 0 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 estimateŝ x on variations of the center frequencyω, 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 β 0 ). In other words, these two sets of confounds can be motivated on principled grounds. We will thus perform a statistical comparison between H 0 and H 1 that marginalizes over the two sets of confounds. 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).

VALIDATING NEURAL STATES ESTIMATE OBTAINED BY sDCM AGAINST EEG
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ω for each channel of each subject, according to Equation (10). This is depicted in  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ω is very well fitted by the model (upper row of Figure 6: coefficient of determination: R 2 = 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 2 = 0.24), which indicates a reasonably good fit (above and beyond confounds). Figure 7 summarizes the result of Bayesian model comparison, in terms of the log-Bayes factor LBF = log p(ω |H 1 ) − log p(ω |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 1 ), whereas there is evidence against H 1 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 1 against the null H 0 ). 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.
Recall that the results summarized in Figures 6, 7 are only shown to illustrate the analysis. In fact, it is important to note 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.
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 1 , irrespective of the particular set of confounds. To confirm this, we ran a random-effect analysis (Stephan et al., 2009) to assess the prevalence of H 1 and H 0 at the group level. Note that the two sets of confounds ("slow drifts" and "full confounds") induce the four following models: 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 1 , irrespective of the definition of the confounds. We have thus partitioned the set of four models above into two families H 1 and H 0 that gather both sorts of confounds, as follows: 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 H 1 is about 0.82. In addition, the exceedance probability of H 1 (versus 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). 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 ) 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 1 (versus H 0 ), 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 includebut 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, theslow-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    Kass and Raftery, 1995). Highlighted cells depict the two subjects that show relative evidence in favor of H 0 . 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 betweenregions 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ω = ω z (0) 1 (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.
where the coefficients A k (t, y, m) obey the following recursive equation (Von Holdt, 1965): The contribution ϒ k (t, y) of the kth spatial mode to the expansion in Equation (A5) is given by: 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. It can be seen form Table A1 that, except for the fundamental mode, the spatial mode contribution ϒ k is of the form: 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 dynamicsz (k) 1 (t, y). We will therefore truncate the contribution to its first term, i.e.,: ϒ 0 (t, y) = S z (0) 1 (t, y) for k = 0, and ϒ k (t, y) ≈z  1 (t, y) for k > 0. This yields the following approximation to the firing rate field:  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: λ k, l = − k 2 + l 2 π 2 / D 2 . The effective frequency cut-off k c is such that k 2 + l 2 ≤ k 2 c . In practice, we use k 2 c = 16 and D = 1 cm.