This paper has presented an extension of the DCM framework to the analysis of phase-coupled data in which a weakly coupled oscillator approach is used to describe dynamic phase changes in a network of oscillators. The WCO approach accommodates signal nonstationarity by using differential equations which describe how changes in phase are driven by pairwise differences in instantaneous phase. This is to be compared with data analysis approaches, such as the PLV or measures derived from autoregressive modelling, which assume signal stationarity.

Previous work on fitting WCO models to experimental data has relied on an evolution map approach (EMA) which is restricted to bivariate data

(Rosenblum and Pikovsky, 2001). A more recent approach considers multivariate data but is restricted to isotropic coupling

(Tokuda et al., 2007). We have shown using simulations that, for bivariate data, the DCM approach is more accurate than EMA. More importantly, DCM can be applied to networks with more than two regions. Inferences can then be made about network structures giving rise to the observed dynamics. For example, whether synchronization processes depend on master-slave (unidirectional) or mutual entrainment (bidirectional) mechanisms.

We have proposed a DCM approach in which each experimental condition of interest is represented using multiple trials of data. This is to be contrasted with DCM for ERPs

(David et al., 2006), for example, in which multiple trials of data are first averaged and then presented to DCM. In the context of phase coupling, the use of multiple trials is necessary so that a range of initial conditions, and therefore different domains of state space, are explored. We have shown using simulations that this leads to improved model inference. Such an approach is not of use in DCM for ERPs, as the underlying dynamical model operates around a known fixed point (zero activation), which is equal to the initial condition.

In the application of DCM to the MEG data, a first order series was used to approximate the PIFs, in which both even and odd terms were allowed to vary between conditions. This allowed DCM to identify both the Fixed Points (FPs; see e.g. ) and dynamics that led to them. In other applications, the FPs may be know a priori, in which case these values can be subtracted from observed data, and DCMs employed using simpler sine interaction PIFs. It may be that these FPs can be identified using multivariate phase clustering methods

(Hutt et al., 2003). A third option here would be to include an additive term in the observation model (Eq.

(14)), to accommodate any relative-phase offsets, and to then always use sine interaction PIFs. This would be suitable if we wished to model phase dynamics with a single, but unknown, stable fixed point.

Although the use of MEG data in this paper is only aimed at a proof of principle for the methodology, the analysis we have performed does appear consistent with findings from rodent neurophysiology in which fronto–cortical neurons were found to spike at consistent phases of the hippocampal theta rhythm, presumably for the transfer of phase-coded information

(Jones and Wilson, 2005). In this paper, in a working memory task, we found that phase dynamics in the memory condition led to near zero-lag correlation between theta activity in fronto–temporal regions. This was not the case for the dynamics in the control condition. Moreover, we can infer that this synchronization was brought about by a master-slave structure in which phase changes in both temporal lobe and frontal cortex were driven by activity in visual cortex.

The following subsections address a number of shortcomings of the current modelling approach, which we hope to address in future publications. We first note that an obvious drawback of the WCO approach is that it ignores amplitude variations, and so provides an incomplete description of brain dynamics. If, however, one is explicitly interested in amplitude changes that give rise to experimentally induced changes of spectral energy then a recent alternative DCM approach may be of use

(Chen et al., 2008).

4.1. Neurophysiological models

The PIFs which parameterise WCOs can be related to neurophysiologically realistic neural network models in a number of ways. In this paper we have considered categorisations of models from dynamical systems theory, in which models are classified according to the type of bifurcation underlying oscillatory behaviour. This leads to specific forms of PRC and PIF, and motivated us to use Fourier expansions for the PIF.

The connection between WCOs and underlying neurophysiological models can, however, be made more explicit. Indeed, given any differential equation model of neuronal or network activity, in which the system operates around a stable limit cycle, PRCs can be numerically evaluated using perturbation or adjoint methods. These are implemented in the XPP or MATCONT software packages

(Govaerts and Sautois, 2006; Ermentrout and Kleinfeld, 2001). Additionally, a simple way of seeing how changes in biological parameters,

, affect network synchronization is then via the derivatives

where

is the PIF

(Govaerts and Sautois, 2006). Use of such derivatives will be investigated in future work.

4.2. Stochastic dynamics

Another way to extend the biological validity of the WCO approach is to model the evolution of, not just a single phase variable in each region

, but a probability density over phases

. These densities can be considered as arising from multiple oscillators within a region and can be specified as solutions of weakly coupled oscillator dynamics based on stochastic differential equations (SDEs). Such an approach has been considered by

Brown et al. (2004) who derive analytic results for a population of oscillators with different initial phases, responding to transient inputs. Such stochastic dynamics can be characterized using Fokker–Planck equations

(Daffertshofer, 1998) or approximated using moment-closure methods

(Ly and Tranchina, 2007; Marreiros et al., 2009). Incorporating such behaviour would require extending the current DCM from deterministic differential equations (DDEs) to SDEs.

4.3. Phase resetting

Tass (2003) has also specified weakly coupled oscillator dynamics using SDEs in which experimental inputs give rise to responses which are transiently synchronized over trials. Here the population of responses is over trials, rather than over multiple oscillators within a single trial. Such dynamics can be accommodated by adding an extra term to Eq.

(11) describing effects of within-trial inputs

. Following

(Tass, 2003) this could take the form

which, provided the input parameters

were sufficiently large, would cause the phase (over trials) to lock at a certain peri-stimulus time point. This would, as suggested by

Makeig et al. (2002), provide a mechanism for the generation of ERF/P components in which system dynamics operate around limit cycles, rather than fixed points as in previous work

(David et al., 2006).

4.4. Conduction delays

In this paper, conduction delays have been absorbed into the representation of the phase interaction function, using a Fourier series approach. In future we will make use of independent sources of information about conduction delays, such as from diffusion imaging or from anatomical databases as in

Ghosh et al. (2008).