1. Introduction

Brain dynamics implement brain functions by organizing brain activity; but we do not directly observe functions, activity, or dynamics in a functional neuroimaging experiment.  Although the successful performance of a task manifests brain functions, we observe only the sensory and motor task variables—not the underlying brain functions that coordinate, monitor, predict, plan, and control these variables.  We observe shadows of brain activity via the spatiotemporal signals obtained from a modality, such as BOLD (blood oxygenation level dependent) signals from fMRI, or fields/potentials from M/EEG.  However, in order to make inferences about the underlying brain activity, we must apply biophysical, physiological, and biostatistical theories to the given task variables and measured signals.  Finally, we clearly do not directly observe the underlying brain dynamics, which are the deterministic and/or stochastic laws that govern brain activity to realize functional requirements.  In all three cases, we must combine the available data with a body of theory in order to make inferences about the interesting, underlying unknowns.

Different modalities measure different signals that may reflect essentially the same underlying brain activity, dynamics, and functions.  Armed with complementary datasets such as fMRI and M/EEG that have different spatial and temporal characteristics, the general question is: Can we make stronger inferences about each of the interesting unknowns via spatiotemporally improved signal detection and parameter estimation?  Clearly, we will need some common theoretical machinery in order to make joint, multimodal inferences.

To begin, let us make an assumption.  When different-session fMRI and M/EEG datasets are available for the same subject who performed realizations of the same task, let us assume that the underlying brain dynamics were statistically equivalent across sessions.  Although the brain activities must have differed in details across sessions, we assume that the underlying dynamics were substantially equivalent in the sense that between-session differences were statistically small compared with other experimentally significant variables.  Based on this testable assumption, our goals are (A) to design cross-modal analyses of both datasets (fMRI and M/EEG) that improve, spatially and temporally, the estimation and detection of the common brain activity, and (B) to characterize the underlying dynamics.  Ultimately, we seek (C) to understand how the dynamics implement brain functions in the context of successful task performances.  While keeping in mind the larger context of problem C, this paper presents an approach for solving problem A that leads to the threshold of problem B.

1.1 Problem C: Spatiotemporal Integration Functions

In the long term, an ideal analysis of spatiotemporal imaging data should permit us to make inferences about brain processes as realizations of brain functions, i.e., as solutions to natural problems that engage the brain.  Thus, we ultimately seek methods that can bridge the gap between (i) raw neuroimaging measurements from various modalities and (ii) descriptions of underlying brain processes in terms of the problems they solve.  Such descriptions call for “functional insight” that is attuned to the “native interests” of the brain regions under study (e.g., Lettvin et al., 1959).  We can approach this objectively if we express a putative functional insight as a testable hypothesis that explicitly formulates the problem that a brain process solves.  Because the analysis we ultimately seek should permit us to test such functional hypotheses, it follows that the technical problem of designing imaging algorithms depends on the scientific problem of formulating functional theories about the natural problems solved by the brain.  In short, our methods must match the subject matter.  On the other hand, because the subject matter comprises unresolved functional issues that our methods should address, we must design the analysis to match a general theoretical framework that is relatively neutral, which can accommodate various specific functional hypotheses.  Some viable candidates for such a framework are: adaptive system theory (Ashby, 1960; Sommerhoff, 1974); schema theory (Arbib et al., 1998); and information theory used as a language to express functional hypotheses, e.g., as used by Phillips and Singer (1997) to formulate the coherent infomax hypothesis of cortical processing.  To our knowledge, these frameworks remain undeveloped along lines that would render them immediately useful for the design of imaging analysis algorithms.  However, in spite of this present lack of a suitable general framework for describing brain function, these basic considerations nevertheless can guide our current efforts toward the long-term goal.

Spatiotemporal integration functions play a central role in solving various problems that engage the brain, which we classify as extrinsic or intrinsic.  “Extrinsic problems” challenge the brain to dynamically organize skeletal and visceral behaviors that can achieve co-adaptations between the organism’s environment and the organism’s physiology.  Two extremes of such co-adaptive processes are environment-directed and physiology-directed brain processes.  Environment-directed processes adaptively control environment variables around the organism by coordinating, in space and time, two streams of brain processes: (a) motor planning and control of skeletal effector variables, and (b) sensory monitoring and prediction of the environment variables.  (We use the terms “coordination”, “planning”, “control”, “monitoring”, and “prediction” to describe objective forms of brain function without ascribing subjective content.)  To illustrate, consider Sommerhoff’s (1974) example of a bird pecking at a grain of seed, the position of which may be subject to wind perturbations: The planned/controlled trajectory of the beak must be coordinated with the monitored/predicted trajectory of the grain.  In a formally similar way, physiology-directed processes modulate the autonomic regulation of physiological variables in the organism by coordinating, in space and time, two streams of brain processes: (a) planning and control of visceral effector variables, and (b) monitoring and prediction of the physiological variables.  For example, efforts of biomedical engineers to externally control hemodynamic variables (Isaka and Sebald, 1993; Palerm et al., 1999) underscore the highly nonlinear nature of these dynamic physiological control problems.  More generally, co-adaptive processes are bi-directional and higher-order in the sense that they coordinate environment-directed (somatic) and physiology-directed (visceral) processes (Yates and Stocker, 1998) as, for example, when a runner’s physiology responds to meet the increased demand of climbing a steep grade. In summary, the requirement for coordination of multiple brain processes in space and time implies that spatiotemporal integration functions play a key role in solving extrinsic problems.

 “Intrinsic problems” are self-referential problems that the brain solves (in the context of solving extrinsic problems) by coordinating brain architectural development processes and brain physiological regulation processes.  Thus, two extreme types of processes that solve intrinsic problems are neurophysiological self-regulation and neuroarchitectural self-organization.  Processes of the former type modulate the regulation of physiological variables in the brain itself, and so are special cases of the physiology-directed processes discussed above.  An example is the regulation of the coupling between brain hemodynamic variables (such as regional cerebral blood flow) and neuronal activity variables (such as mean firing rate).  In order to connect the neuronal and the hemodynamic variables in this example, spatiotemporal integration must span different temporal and spatial scales.  At the other extreme, architectural self-organization (which descends from genetic and embryological processes) implies that the brain “develops itself”, with each developmental stage setting new constraints and providing new possibilities for the next stage (Kelso, 1995; Arbib et al., 1998).  For example, functional segregation constrains the dynamics of the adult brain while simultaneously providing combinatorial possibilities for functional integration.  This modular architecture is a general solution to an intrinsic “multimodal integration” problem (Sporns et al., 2000), which conditions the ongoing solving of particular “dynamic link” architectural problems that arise in the context of solving extrinsic problems (von der Malsberg and Singer, 1988).

The concept of neuronal transients (Friston, 1997; 2000) arises in the context of this intrinsic multimodal integration problem.  Although each specialized area operates on its own characteristic spatial and temporal scale, brain areas collectively coordinate their functions, as necessary.  In order to provide coherent output signals to other brain areas and/or motor systems, each area must therefore resolve spatial ambiguities and temporal asynchronies in its input signals as they arrive from other brain areas and/or sensory modalities.  In particular, finite intervals of integration are required in order to resolve uncertainties that arise from temporal asynchrony (Pöppel, 1997).  This necessitates asynchronous coupling between neuronal transients in different brain regions, which in turn requires nonlinear interactions (Friston, 2000).  Others (e.g., Poggio and Reichardt, 1981; Koch, 1999) have emphasized the importance of nonlinearities, such as multiplication of two or more variables, for realizing computational requirements for basic functions, such as pattern recognition, coordination, and control.  Thus at minimum, a model of the nonlinear, asynchronous dynamics of neuronal transients must be multilinear (permitting multiplications) and have transient memory (permitting temporal interactions).  The satisfaction of these conditions by the deterministic Volterra model (Poggio and Reichardt, 1981; Friston, 2000) sets the stage for using nonlinear system identification methods to approach problem A.

1.2 Problem A: Solution Strategy

Our strategy for solving problem A is as follows.  Task-related M/EEG or fMRI analyses treat the brain as a dynamic system, where the input variables are sensory and/or motor task variables, the output variables are measured M/EEG or fMRI signals, and the analysis yields estimates of system parameters that maximize the dynamic coupling between the input and the output variables.  Let us conduct a thought experiment that suspends the physical impossibility of simultaneous MEG and fMRI:  Suppose that we have an ideal, simultaneous M/EEG and fMRI dataset.  We could perform cross-modal analyses with such a dataset—call these “fMRI-related M/EEG” and “M/EEG-related fMRI”—by using the measurements of one modality as the input variables, in place of the sensory and/or motor task variables.  Here we suppose that the analysis method can handle continuous as well as discrete (“event-related”) input variables.  An estimator of M/EEG source current activity at a given location in the brain, for a given orientation, has the form of a linear combination of channels.  Using an fMRI BOLD signal at a voxel (or a local average of voxels), we can design an fMRI-related M/EEG analysis that “tunes” the linear estimator of source current activity to maximize the dynamic coupling between the estimated activity and the BOLD signal.  Then we can use this M/EEG-derived estimate of neuronal activity to make an M/EEG-related fMRI analysis that estimates the hemodynamic response to neuronal activity at the voxel (or local average of voxels).  Because task variables are remote from the hemodynamic system, the “neuronal-derived” hemodynamic response will differ from the task-derived hemodynamic response.  Consequently, matched filter signal detection using the neuronal-derived hemodynamic response may be more powerful for constructing fMRI statistic images of activated voxels.  In summary, fMRI spatially focuses M/EEG source time series estimates, which we use to produce an improved estimate of the fMRI hemodynamic response to neuronal activity and to derive a new fMRI statistic image.  Thus, the strategy is to use the temporal strengths of M/EEG to enhance the spatial strengths of fMRI, and vice versa.

Figure 1.  Schematic of the overall strategy for solving problem A.  System models that maximize the dynamic coupling between task variables and measured signals are independently identified for same-subject, same-task, different-session fMRI and M/EEG datasets.  These models are used to simulate a simultaneous fMRI and M/EEG dataset.  Then cross-modal system models are identified by maximizing the dynamic coupling between the simulated fMRI and M/EEG signals.

To enable this strategy, we must handle the typical situation where we do not have simultaneous datasets.  Our auxiliary strategy is to use the same dynamic system modeling methods (that enable the fMRI-related M/EEG and M/EEG-related fMRI analyses) to make task-related dynamic models, separately for each modality.  Then we can use these dynamic models to simulate “virtual simultaneous” fMRI and M/EEG datasets.  Figure 1 summarizes the overall strategy.

It remains to specify the nature of the system model.  Based on the previous section’s considerations, we require a nonlinear model with transient memory, such as the Volterra model.  Because problem B is concerned with the underlying dynamics of brain activity, we approach the question of selecting a dynamic system model in the context of a formal statement of problem B.

1.3 Problem B: Dynamics and Observables

At a hypothetical and formal level, we can state problem B in terms of a system of equations that comprise one state equation (1), which governs the spatiotemporal dynamics of brain activity, and one observation equation (2) per modality (Ozaki, 1985; Valdes et al., 1999):

(1)

(2)

where  is the (unknown) state of neuronal activity at location and time ;  is an (unknown) differential operator; is a (known) input vector of control variables; is an (unknown) deterministic “driving” function; and  is an (unknown) stochastic term.  For the observation equations,  is an observed measurement we obtain from neuroimaging modality ;  is a modality-dependent transfer function that maps the neuronal state onto the observables via biophysics, physiology, and instrumentation; and  is a noise term.  Then given a set of measurements , control variables , and knowledge of the transfer functions , we seek to infer first  in both space and time, then ultimately the dynamics  and  (i.e. the specific form of the state equation itself).  When multimodal information is available, this system of equations implies that our ultimate goal is to obtain a simultaneous solution, i.e., a joint inference about the interesting unknowns.

Simplifications of equation (1) lead to three working models of the dynamics of task-related neuronal transients:

(3)

(4)

(5)

In all three cases, we expand the driving function  as a Volterra series (see equation (6)).  Thus, the pure Volterra model of equation (3) directly expands the brain activity as a deterministic Volterra series of the task variables  (plus the stochastic activity).  In the differential-Volterra model of equation (4), a differential operator  characterizes the driven process, which is a derivation of the neuronal transient activity.  In the self-history Volterra model of equation (5), we include the brain activity as a variable that contributes to its own driving function.  Treatment of the stochastic term involves modeling (e.g., autoregressive modeling) of the task-unrelated brain (and other noise) processes to account for serially correlated background activity (Zarahn et al., 1997; Friston, Josephs, et al., 2000).

Generalizations (4) and (5) may contribute important theoretical and practical advances over the pure Volterra model.  For example, self-history and nonlinear differential operators are two ways of modeling deterministic-stochastic interactions between  and , which may provide a framework for explaining the observations of Brandt and Jansen (1991) that “deterministic” visual evoked responses interact with their “stochastic” prestimulus baselines.  Deterministic-stochastic interactions may also provide an explanatory framework for well-known (but poorly understood) response variability phenomena, such as statistically significant single-trial ERP latency jitter (e.g., Pflieger, 1991).  In practice, both approaches have the potential for more efficient modeling of the system: Recent self-history may reduce the history required for task variables; and a few differential parameters may “leverage” expensive Volterra kernels.

1.4 Possibilities for Combined Analysis of fMRI and M/EEG Datasets

Now we turn specifically to fMRI and M/EEG modalities.  Figure 2 diagrams the relationship between the measured M/EEG and fMRI signals and the underlying neuronal transient activity.

Figure 2.  Common underlying sources of M/EEG and fMRI signals.

Neuronal transient activity  (a function of space  and time ) underlies observed M/EEG signals across channels  (functions of time) and fMRI BOLD signals  (a function of space and time).  We assume that the task-related dynamics of neuronal transients have the form of a differential equation , where  is a differential operator,  is a deterministic driving function of the task variables  (functions of time) and self-history, and  represents stochastic influences (a function of space and time).  The transfer functions  and  describe the net effect of biophysical, physiological, and instrumental factors that convert the underlying neuronal transient activity to the measured M/EEG and fMRI signals, respectively.  The dashed ellipse represents a subspace of neuronal transients that can be sources of M/EEG signals; this subspace is the domain of the M/EEG transfer function .  Likewise, the dotted ellipse represents the source domain of the fMRI signals.  We assume that the M/EEG and fMRI source domains overlap, i.e., that there is a neuronal transient subspace that is the source of both M/EEG and fMRI BOLD signals.  The diamond represents the neuronal transient subspace sampled by an experiment.  Subspace  represents those neuronal transients in the experiment that affect M/EEG, but not fMRI, signals.  Subspace  represents neuronal transients that affect fMRI, but not M/EEG, signals.  The common subspace  represents neuronal transients that affect both M/EEG and fMRI signals (whereas subspace  represents those that affect neither).  The image of the common subspace  is  in the M/EEG signal space, and  in the fMRI signal space.  The signal subspaces  (in M/EEG signal space) and   (in fMRI signal space) are images of the disjoint neuronal activity subspaces  and . 

It is generally believed that postsynaptic neuronal currents in an open field dendritic arrangement, such as the apical dendrites of pyramidal cells, form the basis for the scalp potentials and fields measured as M/EEG (Lorente de Nó, 1947).  We obtain the observation equation that relates these currents and fields from Maxwell's equations, solved numerically for the volume conductor that represents the head. This is often referred to as the forward problem.  Volume conduction is linear and, for signal frequencies less than ~1000 Hz, instantaneous (quasi-static approximation; Malmivuo and Plonsey, 1995).  Forward solution methods include analytic/spherical (e.g., Arthur and Geselowitz, 1970), boundary element (e.g., Hämäläinen and Sarvas, 1987), and finite element (e.g., Awada and Greenblatt, 1999).  However, Helmholtz (1853) showed that there is no unique solution to the inverse problem of estimating these neuronal currents from measured extracranial fields.  All such solutions are inherently model-dependent.

The most commonly employed signal in fMRI is blood oxygenation level dependent (BOLD) contrast: a change in water proton signals caused by changes in the relative concentration of venous blood deoxyhemoglobin, which serves as an endogenous paramagnetic contrast agent (Kwong et al., 1992; Ogawa et al., 1992).  The fMRI observation equation must incorporate physiological mechanisms that relate neuronal activity to BOLD signals via local blood oxygenation, regional blood flow, and regional blood volume—the theoretical details of which are current topics of research (Buxton et al., 1998).  There are three main phenomenological components of the (visual) evoked BOLD response.  A small signal decrease (initial dip) that peaks about 2 s after stimulus onset (not routinely seen in low fields) results from a transient decrease of blood oxygenation due to local tissue demands (Menon et al., 1995).  The largest BOLD component, a positive overshoot that peaks about 5 s, reflects a decrease in relative deoxyhemoglobin concentration caused by increased local blood flow in excess of tissue demands.  Finally, a negative undershoot requires up to 30 s to return to baseline.  Timing details of the BOLD signal are task dependent, and vary across brain regions.  Because the BOLD response to task events is generally not identical with the response to neuronal activity, recent efforts (Ogawa et al., 2000) have sought to relate the BOLD response empirically to neuronal activity as reflected by EEG.  It is unknown whether neural-hemodynamic decoupling that occurs during normal non-REM sleep (Klingelhofer et al., 1995) and normal aging (D’Esposito et al., 1999) may play a wider role during other normal functions, or perhaps in some brain disorders.

What, then, can we hope to learn about the underlying neuronal transients by combining information from M/EEG and fMRI?  Of course, we cannot expect to gain information about subspace  of Figure 2, which represents brain activity that is invisible to both modalities. Similarly, because all information we have about subspace  is contained in the M/EEG data, and all information we have about subspace  is contained in the fMRI data, we cannot expect integrated analyses to cast new light here either: The most we can hope for is to delineate the boundaries of these disjoint subspaces.  Possibilities for , , and  include the following: MEG is insensitive to radial current sources; both EEG and MEG are insensitive to currents generated in closed field dendritic arrangements (e.g., basal dendrites of pyramidal cells), and more generally to unsynchronized local activity (Nunez, 1995); fMRI is insensitive to small or brief changes in neuronal activity that do not produce detectable BOLD responses; and neural-hemodynamic decoupling may occur to some extent in normal and abnormal brain functioning.  For our purposes, it is the common subspace , which is visible to both modalities, that provides the most promising territory for joint analyses.  Here we can reasonably hope that M/EEG and fMRI will provide complementary information in space and time about the same underlying brain activity.  Thus, two fundamental objectives of integrated analysis are (a) to identify the image subspaces  and  in the combined M/EEG and fMRI datasets that refer to the same underlying brain activity subspace c, and (b) to exploit complementary information in  and  to gain new knowledge about c and its underlying dynamics.

Methods for identifying the common subspace and exploiting complementarity can utilize overlapping spatial, temporal, or spatiotemporal information.  Given that M/EEG excels in the temporal domain and fMRI excels in the spatial domain, we expect that an optimal approach should incorporate both spatial and temporal information.  All of the spatial information in M/EEG comes from the locations and geometry of the sensors, which we can co-register with structural MRI.  Methods for solving the M/EEG inverse problem have this anatomical information at their disposal, and can map source estimates to an anatomical space in common with fMRI.  Some studies have reported convergent spatial results of independent M/EEG and fMRI analyses (e.g., Grimm et al., 1998; Ball et al., 1999).  Other studies have sought to utilize fMRI statistic images to constrain M/EEG inverse solutions by dipole seeding (Ahlfors et al., 1999; Opitz et al., 1999; Scherg et al., 1999) or regularization of a linear estimator (Liu et al., 1998; Babiloni et al., 1999).  Wagner et al. (1999) conducted a simulation study of both approaches.  An important goal of these methods is to use fMRI to bolster confidence in the spatial aspect of M/EEG solutions so that the latter can be used to extract source time series on a millisecond time scale, which in turn can serve as the basis for dynamic analysis of large-scale neural networks (Thatcher et al., 1994).  Ioannides (1999) cautioned that statistic images obtained from fMRI (or PET) block designs might be unsuitable for constraining the analysis of MEG transients due to disparate temporal scales.  This criticism does not necessarily apply to event-related fMRI designs.  However, in relation to the method presented in this paper, we note that (to our knowledge) no approach that uses fMRI statistic images to improve M/EEG source estimates has yet made explicit use of the temporal information contained in the BOLD signal: Although the datasets have been co-registered in space, they have not been co-registered in time, i.e., temporally aligned.  RS Menon et al. (1998) have shown that relative BOLD signal onsets in different regions may resolve timing sequences within tens of milliseconds.  Similarly, Wildgruber et al. (1997) associated two phases of movement-related potentials (BP and motor potential) with sequential event-related average hemodynamic response onsets (in SMA and M1, respectively).  V Menon et al. (1997) used same-subject, same-task event-related EEG and fMRI responses to gain information about generators of the P300 ERP component.  This study demonstrated that temporal alignment of the linear event-related EEG and fMRI responses is both possible and useful.  Guy et al. (1999) temporally aligned EEG and fMRI responses with respect to periodic visual stimuli, and found quasi-harmonic components of the stimulus driving frequency in both modalities, which demonstrate nonlinearities in the brain response.

The temporal alignment problem would vanish if fMRI acquisition were truly simultaneous with M/EEG.  Of course, this is physically impossible with MEG.  Ives et al. (1993) demonstrated this possibility with ongoing EEG, and Bonmasser et al. (1999) measured visual evoked potentials concurrently with fMRI.  Suitable controls can effectively eliminate image artifacts caused by the EEG apparatus (Krakow et al., 2000).  Artifacts in the EEG signal occur in the MR environment due to the following sources: the static magnetic field interacts with the EEG instrument; subject movements in the static field induce currents; the ballistocardiogram (pulse) artifact; mechanical vibrations; and especially inductive effects of gradient switching during MR sequences (Lemieux et al., 1999).  Various approaches, some hardware and others software, can reportedly deal with these issues.  Another unresolved issue is the possibility that the MR environment may alter somewhat the EEG signal itself at a physiological level.  With reference to the right half of Figure 1, we note that the cross-modal analyses described here are pertinent to actual simultaneous datasets.

However, because we do not assume simultaneous M/EEG and fMRI acquisitions, it is generally necessary to perform temporal alignment.  This, in turn, relies on inter-session reliability of the task-related brain responses, which is a testable assumption (Noll et al., 1997; Mangun et al., 1998; Aguirre et al., 1998).  The remaining temporal alignment difficulties would vanish if brain responses were strictly linear functions of event-related task variables.  Then same-type events would naturally align their respective event-related averages (Menon et al., 1997), which approximate linear impulse response functions.  However, this direct alignment approach cannot handle dynamic nonlinearities, which are theoretically important (Poggio and Reichardt, 1981; Friston, 2000) and empirically observed (Vazquez and Noll, 1998; Friston et al., 1998; Guy et al., 1999; Ances et al., 2000).  Furthermore, our strategy for cross-modal analysis requires that we handle continuous, as well as event-related, task variables.  A viable tactic for solving both the temporal alignment problem and the problem of continuous task variables is to use nonlinear system identification techniques (Schetzen, 1980; Bendat, 1998) (a) to model M/EEG and fMRI responses and interactions separately with respect to the same class of task realizations (which may be either event-related or continuous) and (b) to conduct virtual simultaneous experiments by driving both nonlinear models in parallel with the same task realizations.

The approach detailed below is particularly indebted to Friston et al. (1998), who (a) demonstrated the viability of the pure Volterra approach for modeling fMRI time series using second order kernels, (b) conducted a virtual experiment with the model, and (c) placed the Volterra model in the familiar context of the general linear model.  As an example of (c), the authors constructed SPM{F} images of voxels with significant nonlinear interactions.  Furthermore, the Volterra approach naturally extends to the analysis of effective connectivity and modulatory interactions in large-scale neural networks (Friston and Büchel, 2000) by replacing the task variables with inputs from other brain regions.

2. Identification of Volterra-Based System Models

In this section, we describe methods for identifying the parameters of the three nonlinear, time-invariant system models represented by equations (3), (4), and (5).  We start with the pure Volterra model (equation (6)), and subsequently extend to the differential-Volterra model (equation (22)) and the self-history Volterra model (equation (25)).  These identification methods utilize “raw” (not averaged or statistically processed) continuous signals.  To simplify the description below, we assume that there are no temporal gaps in the data.  However, in practice, the data may consist of multiple runs, and each run may contain deleted segments due to artifact rejection.  Suitable bookkeeping of the time variable can handle these discontinuities.

3. The Pure Volterra Model

A Volterra series (Schetzen, 1980; Bendat, 1998) of  task variables  can model the task-related responses and interactions of a measured time series , such as the BOLD fMRI signal at a voxel (Friston et al., 1998) or a linearly derived M/EEG signal, as follows:

(6)

Here  arises from a random background process, which is primarily task-unrelated brain activity.  We have included low-order powers of time with coefficients  in order to account for  slow drifts plus a constant, which we regard as uninteresting.  The experimentally driven responses and interactions of interest are modeled by the Volterra kernels  (order 1) and  (order 2).  Task variables linearly convolve with their first order kernels, and product pairs of task variables bilinearly convolve with their second order kernels.  Higher order Volterra kernels n-linearly convolve with product n-tuples of task variables; however, the number of unknown parameters that we must estimate increases sharply with the order (see equation (12) ).

A first order kernel  represents the linear response associated with task variable  at time lag .  We assume that responses are transient, i.e., the kernel is zero except for a finite interval of time lags.  If the task variable corresponds with a stimulus, then we assume that the brain response follows the stimulus, which implies that a sensory kernel is identically zero for all non-positive lags.  If the task variable corresponds with a behavior, then we assume that the brain “response” precedes the behavior, which implies that a motor kernel is identically zero for all non-negative lags.  For tasks that involve tight coupling between a stimulus and a behavior, we can allow sensorimotor kernels that are nonzero for both positive and negative lags.  Whenever possible, however, we prefer to preserve causality by teasing apart the sensory and motor responses.  In general, we assume that each task variable has an associated effective range of time lags, such that a kernel is identically zero outside the effective range.

A second order kernel  represents the bilinear interaction between task variables  and  at respective time lags  and .  Kernels with  represent quadratic self-interactions.  We assume that a second order kernel can be nonzero only when both time lags are within the effective ranges of their associated task variables.

For each task variable , we can specify a series of basis functions  that enable us to approximate the first order kernel as a weighted sum (Schetzen, 1980; Friston et al., 1998):

(7)

(We discuss below a method for optimizing these basis functions for particular task variables.)  Using the same basis functions, we similarly expand the second order kernels as

(8)

Noting that the b and g coefficients are not functions of time lags, we next substitute approximations (7) and (8) into equation (6) with an eye towards joint estimation of the a (uninteresting drift), b (interesting response), and g (interesting interaction) parameters:

(9)

By convolving task variables with basis functions, we derive new time series

(10)

Finally, we simplify equation (9) by substituting definition (10):

(11)

Because we obtain the  functions by convolving known task variables with predetermined basis functions, we thus can jointly estimate, via generalized least squares, the a (drift), b (linear response), and g (nonlinear interaction) parameters.  Equations (7) and (8) reconstruct the first and second order Volterra kernels from these parameters.  We note that the self-interaction kernels are symmetric.

4. Gauss-Markov Estimation of Volterra Kernel Parameters

The total number  of estimated parameters for the second order Volterra model is

(12)

This total includes a constant, the a trend coefficients, the b response parameters, the symmetric self-interaction g parameters, and the cross-interaction g parameters, from which we form a vector  of length

(13)

Corresponding with the column vector , we form a row vector of length  that reflects the experimental situation at time  (in terms of basis-convolved task variables):

(14)

For an experimental run consisting of  time samples, the full  design matrix  is

(15)

Collecting the measured signals for all time samples yields a data vector of length :

(16)

Likewise, we form an error vector  of length :

(17)

We may rewrite equation (11) in matrix form as

(18)

The time series  has a serial correlation structure that derives from “background” neural activity, other physiological activity, motion effects, and/or scanner/system noise (Zarahn et al., 1997).  This, in turn, implies that the frequency spectrum is not flat (i.e.,  is not white noise).  If the error process is approximately gaussian and stationary, we can solve (18) via ordinary least squares estimation if we first prewhiten the data.  In the frequency domain, we divide by the amplitude at each frequency, which flattens (“whitens”) the spectrum.  The time domain equivalent of prewhitening can be formally written as equation (19), where is the  variance-covariance matrix for the error time series.  Stationarity of the background activity implies that each autocovariance depends only on the lag .  If we extrapolate autocovariances beyond a specified lag  using the principle of maximum entropy, then we can use autoregressive (AR) methods to design a prewhitening filter in the time domain by inverting an AR model of the error process, which involves only autocovariances with lags less than  (Friston, Josephs et al., 2000; Appendix A).  Applying the prewhitening filter  to each term of equation (18), we obtain

(19)

Under the assumptions,  is identically and independently distributed, in which case the Gauss-Markov estimate for the parameters is

(20)

Here the “+” superscript represents the pseudoinverse operation.

In practice, we typically must use the same data to estimate the autocovariances of the background brain activity with lags less than .  A feasible generalized least squares approach (Dale et al., 2000) obtains a series of estimates, , where  is computed using  (or another initial estimate);  is computed using the approximation , where we form  by calculating lag covariances from the residual vector ; and so on.  This series converges to mutually consistent estimates of  and , given the data and the modeling assumptions.  Friston, Josephs et al. (2000) have recently discussed caveats regarding the effect of prewhitening filters on statistical bias, which we discuss further below.

Although not the only viable approach for solving equation (20), singular value decomposition (SVD) of  conveniently provides (a) a solution that is stable, and (b) standard errors in the fitted parameters, which reflect goodness-of-fit of the model to the data (Press et al., 1992, §15.4).  Clearly, this analysis fits the general linear model (GLM).  As Friston et al. (1998) have demonstrated, the statistical machinery of the GLM can be combined with Volterra models to construct SPM{F} maps (Worsley 1994) of (for example) brain regions that have significant nonlinear kernels.

4.1 Tailoring of Basis Functions

The total number of estimated parameters  of equation (12) depends primarily on the number of task variables  and the number of basis functions utilized per task variable.  We wish to minimize  in order to obtain stable estimates from a limited number of data samples.  However, in opposition to this goal, we also seek to minimize the approximation error of equations (7) and (8).  For fixed  (an index of task complexity), decreasing the number of basis functions decreases , whereas increasing the number of basis functions decreases approximation error.  To optimize this tradeoff, we strive to maximize the ability of a given number of basis functions to represent Volterra kernels.  Since each set of basis functions is associated with a task variable and its first order kernel, we take an empirical approach that tailors the representation of the first order kernel.  We base this approach on the observation that an average event-related response approximates the impulse response (which is the first order kernel) so long as the task events occur in temporal isolation.  Thus, if we have “calibration data” for each task variable, consisting of isolated task events, then we can make an event-related average to which we tailor a set of parameterized basis functions.  In particular, we propose to use  basis functions that have the form of a Poisson distribution:

(21)

where ;  is a “characteristic latency” that roughly corresponds to the peak latency of the event-related average response; and  is a “characteristic delay” that roughly corresponds with the onset of the response.  The fundamental basis function () can fit a “main peak”; the  functions can fit early oscillations that occur before ; and the  functions can fit slow waves that occur after .  For each number  there are optimal values of  and  for approximating the average response.  Thus, for each task variable, we seek the smallest  that produces an adequate approximation.

4.2 The Differential-Volterra Model

The differential-Volterra model extends the Volterra model (equation (6)) by applying a differential operator  to the data term on the left-hand side (LHS) of equation (6):

(22)

where

(23)

That is,

(24)

Here  is the highest derivative (order) of the differential operator , which has constant coefficients .  The  parameters are the first-degree (linear) coefficients; the  parameters are the second-degree (quadratic-bilinear) coefficients; etc.  The revised LHS is a derived data term, obtained by forming linear and higher degree combinations of the measured time series with its own temporal derivatives.  To prevent the derived data term from vanishing, we require  while considering order .  Thus, the pure Volterra model of equation (6) is the special case of equation (22) when the degree is 1 and the order  is 0  (noting that ).  Whereas the error term  for equation (6) is the error of the original measured data, we note that the error term for equation (22) is the error of the derived data.  For a first-degree operator, there are  undetermined linear coefficients.  For a second-degree operator, there are, in addition,  quadratic coefficients, and  bilinear coefficients.  If the order , for example, then the number of unknown  parameters is 2 for a first-degree operator, and 8 for a second-degree operator.

A differential operator applied to the data permits us, using relatively few additional parameters, to complement the model of the driving process (which is given on the RHS) by specifying a model of the driven process (which is given on the LHS).  In short, the RHS specifies how to drive some process, the nature of which the LHS specifies.  For example, the pure Volterra model of equation (6) assumes that the brain process driven by the task variables is identical to our measurement process.  However, it is equally plausible to conjecture that the task variables drive changes in the brain, as reflected by first derivatives of the measurement process.  If the original measurements and their derivatives are both plausible candidates for the driven brain process, then any linear combination of the two is yet another plausible candidate.  We can extend this to higher order derivatives.  Higher degree combinations give the differential operator a nonlinear character, which (as we noted above) may have theoretical significance by permitting the deterministic and stochastic terms on the RHS to interact. These considerations motivate the rationale for using the differential operator of equation (23) to specify the nature of the driven process.

Although Gauss-Markov methods cannot directly estimate the l parameters, i.e., we cannot simply append them to the  vector of equation (13), we can use the identical methods detailed above as part of an iterative approach.  The steps of this procedure are as follows (for a given degree and order , starting with degree 1 and order 0):

1.                  Numerically evaluate the time derivatives  from the data .  Initially set , and set all other coefficients to zero.

2.                  Use the Gauss-Markov methods detailed above, substituting  for , to obtain an estimate  for the RHS parameters.

3.                  Compute the modeled time series , which we obtain by applying the design matrix to the RHS parameter estimate and evaluating at each time .

4.                  Use numerical methods to solve the ordinary differential equation  for .

5.                  Compute the mean square error of the residual, .

6.                  Employ a suitable, general minimization method that iterates steps 2-5 to find  with minimum , subject to .  We note that with low degree and order,  consists of relatively few parameters compared with the number  of Volterra kernel parameters.

By design, the above steps start with the basic Volterra model (equation (6)), and generally will find a differential-Volterra model (equation (22)) with smaller residual error.  For low degree and order, we can attribute any improvement gained to the use of just a few additional parameters (relative to the number  of Volterra kernel parameters), which effectively change what the Volterra kernels are modeling.  Thus with moderation, we expect that this procedure should not risk overmodeling the data.  This point is important in regard to the objective of finding improved models of the data that generalize across acquisition sessions, i.e., we must avoid models that are overly specific to one session.  Thus, until we have characterized these methods, it will be necessary to acquire multiple session datasets for each modality in order to determine model order-degree that optimizes generalization across sessions.

Steps 4 and 6 depend on unspecified numerical methods for solving ordinary differential equations and minimization of multidimensional functions, respectively.  Here we are contemplating a small number of differential parameters (1-8), so we do not expect that the choice of minimization method is critical.  We can use analytic methods to solve second order linear (first-degree) differential equations (two unknown coefficients), and Laplace transform methods for higher order linear equations.  The quadratic-bilinear (second-degree) equations call for different methods, which we do not discuss here.

4.3 The Self-History Volterra Model

The self-history Volterra model extends the pure Volterra model (equation (6)) by including the past of the response variable  in the Volterra expansion for the present:

(25)

Here we have added a linear self-response kernel, a quadratic self-interaction kernel, and bilinear task-response interaction kernels.  Unlike the differential-Volterra model of equation (22), we can append the new parameters of the self-history Volterra model to the  vector of equation (13).  Although we cannot empirically “calibrate” basis functions as described above, the Laguerre functions are a natural basis for expanding these new kernels (Schetzen, 1980).  One caveat of the self-history approach is the introduction of measurement errors into the input variables (Press et al., 1992, §15.3).  Equation (12) gives the number of parameters to estimate, with the understanding that there is the equivalent of one extra task variable. 

Self-history Volterra model identification is computationally simpler than differential-Volterra model identification, although the latter typically requires fewer extra parameters than the former.  We can combine both models to make a self-history differential-Volterra model, which has the form of equation (1).  In practice, the key consideration for choosing from among the four model types (pure Volterra, differential-Volterra, self-history Volterra, and self-history differential-Volterra) is the ability to efficiently generalize across sessions, i.e., the ability to obtain a good fit of data acquired in one session using model parameters identified with a different session’s data.  All else being equal, we prefer the simplest model.  (We do not address these empirical issues here.)

Task-Related and Cross-Modal Analyses

At this point, we can implement the strategy summarized in Figure 1, where we substitute one of the four model types for “system model”.  If the task-related fMRI and M/EEG models (on the left side of Figure 1) efficiently generalize across sessions, then we can drive them simultaneously with the same simulated task variables to generate a virtual simultaneous dataset.  Although we used the stochastic component to identify the system models, we can omit this during the simulation.  In this case, the simulated data will reflect only the deterministic response, which is free from the influences of “brain noise”.

The goal of fMRI-related M/EEG analysis at a particular location is to determine a set of channel weights such that the derived M/EEG signal has maximum dynamic coupling with the fMRI BOLD signal.  We accomplish this in the context of identifying a system model, where the input is a BOLD signal at the location (perhaps averaged in a local neighborhood).  The measured signal in equations (6), (22), or (25) is replaced by

(26)

where  is a M/EEG channel vector at time , and  is a row of channel weights that we aim to determine.  More generally,  may contain up to three rows that permit up to three source current orientations at the location.  The norm operation converts the resultant one-to-three dimensional vector to a nonnegative number, which implicitly assumes that the BOLD signal is not orientation sensitive.

We can estimate  iteratively, as follows.  On the first iteration, we set  to a theoretical estimator of local activity, such as a beamformer (e.g., Sekihara et al., 2000).  Having fixed , we calculate , identify a system model with respect to the BOLD signal input, and measure goodness-of-fit.  On subsequent iterations, a multidimensional minimization method searches for  parameters that optimize goodness-of-fit.  The associated best-fitting model maximizes the dynamic coupling between the BOLD signal and the source activity estimated by (26).  In practice, we can substantially reduce the number of unknown parameters, i.e., the elements of , by representing  as the sum of a few principal components.

The goal of M/EEG-related fMRI analysis is to provide a more precise temporal context for statistical detection of activations in space by modeling the local BOLD response system using input obtained from the M/EEG-derived signal of equation (26).  The local BOLD signal depends on three main physiological variables: local oxygen metabolic rate, which is a consequence of neuronal activity; regional cerebral blood flow, which is regulated by neuronal activity; and regional cerebral blood volume, which depends on the blood flow (Buxton, 1998).  Thus, the local hemodynamic/metabolic system receives neuronal activity as input and produces BOLD signal as output.  “Activation” means that this system is active, i.e., that there is significant dynamic coupling between the input neuronal activity and the output BOLD signal.  We can test the hypothesis of significant activation in the context of dynamic modeling by forming an F-statistic as the ratio of (a) the magnitude of the projection of the data into the model space over (b) the magnitude of the residual error (Worsley and Friston, 1995; Buxton et al., 2000).  However, we cannot directly model the hemodynamic/metabolic system with respect to the unobserved neuronal activity.  Because the known task variables drive the unknown neuronal activities of interest, we typically must replace the latter with the former as input to the dynamic model.  On the other hand, these sensory and motor task variables are at least one stage removed from the actual neuronal activity of a brain region.  More typically, the neuronal activities of other brain regions intervene.  Moreover, whereas the task variables are uniform across voxels, actual neuronal activity varies from region to region.  Therefore, it may be highly beneficial to replace the task variables with region-by-region estimates of neuronal activity, such as those derived from M/EEG.  Although imperfect, these temporally improved local estimates of the task-related neuronal activity may yield increased statistical power for detecting activated voxels, especially for brain regions that are most removed from the sensory and motor variables.

5. Discussion

In brief, we have outlined a solution method for problem A that uses the framework of nonlinear system identification to simulate virtual simultaneous data, and then to conduct cross-modal analyses.  The fMRI-related M/EEG analysis permits us to “tune” an inverse estimator for source activity at a location.  Although BOLD responses lag behind the neuronal responses that M/EEG detects, we nevertheless can formally model the M/EEG signal as an output to a dynamic system with a BOLD signal input.  This is not a causal model, but the machinery permits us to empirically derive an estimator of source activity that maximizes the local dynamic coupling between BOLD and M/EEG signals.  M/EEG-related fMRI analysis models the BOLD signal as the output of a local hemodynamic/metabolic system that has M/EEG-derived signals as inputs.  This system is presumably causal, and it permits us to estimate the possibly nonlinear hemodynamic responsiveness of fMRI to neuronal activity.

We have not here directly addressed problem B, which is to use both sets of measurements to jointly estimate the underlying dynamics, which, in turn, requires a joint spatiotemporal estimate of the neuronal activity.  However, we have made some progress toward problem B by virtue of attending to problem A.  In particular, the state dynamics equation (1) has motivated our approach, and we have made some theoretical progress toward characterizing the observation equations (2): The cross-modal analysis produces an fMRI-restricted version of the inverse M/EEG transfer function  (embodied in the matrix ) and an M/EEG-restricted version of the forward fMRI BOLD transfer function  (embodied in the nonlinear system model of the M/EEG-related fMRI analysis).  These transfer functions are mutually restricted in the sense that they are relative to each other, having both been “tuned” to the common subspace  of neuronal activity that is strongly coupled with the measured signals of both modalities.  Nevertheless, we are positioned to contemplate a “simultaneous solution” to equations (1) and (2).

The question of optimal experimental design for model identification (parameter estimation) and associated statistics (signal detection) has been a recent topic of discussion in the fMRI literature.  Efficient parameter estimation utilizes randomized, short-latency inter-event intervals (Dale, 1999); the conventional block design with fixed intervals is optimal for statistical signal detection (Friston et al., 1999); and a hybrid block-probability approach represents a compromise (Friston et al., 1999; Buxton et al., 2000).  Discussions of the estimation-detection tradeoff have considered discrete designs that consist of binary point processes of events.  Identification of continuous-input physiological systems via white noise analysis (e.g., Marmarelis and Naka, 1974) is a purely random approach, which parallels the efficient estimators for discrete, binary inputs.  Experimental design via white noise inputs clearly represents the “least designed design” that attempts to sample the entire input space without bias.  However, actual neurons and organisms have extremely biased interests in “remote corners” of input space (Lettvin et al., 1959; Hubel and Wiesel, 1959), and tend to “tune out” random processes.  Because statistical signal detection relies on robust responses, these considerations suggest that we ultimately may desire an adaptive, real-time approach to task design that combines random “importance sampling” with functional hypothesis testing (problem C).  That is, in the context of an explicit hypothesis about an objective function that a brain region putatively optimizes, we might use online statistical signal detection to identify a subspace of task variables that we infer are interesting to a brain region after observing its robust response.  Such an adaptive importance sampling design should dynamically resolve the discrepancy between optimal estimation and optimal detection (which discrepancy signals that the current design is not functionally optimal).

A parallel issue concerns the effect of using (or not using) estimated serial correlations of the background noise on (i) the efficiency of parameter estimation versus (ii) robustness of signal detection.  The Gauss-Markov estimator described above makes use of estimated temporal covariances for more efficient parameter estimation; however, the associated signal detection statistics are not robust with respect to errors in the temporal covariance estimates (Friston, Josephs et al., 2000).  A key reason for this estimation-detection conflict is the requirement of gaussian field theory to use the same statistical model across all voxels (Worsley, 1994; Worsley and Friston, 1995), whereas temporal covariances of the background activity are spatially nonstationary (Zarahn et al., 1997).  However, except for the final step of forming a statistic image, our primary concern here is with a third criterion, which we can evaluate empirically: (iii) within-subject, between-session test-retest reliability of the dynamic models (Genovese et al., 1997; Noll et al., 1997).  Within-session efficiency of parameter estimation does not guarantee that a model will generalize across sessions.  Only inter-session stability can warrant the use of specific dynamic modeling features (such as voxel-specific noise estimates, differential operators, or self-history kernels) for conducting virtual simultaneous experiments and cross-modal analyses.

References