\chapter{Dynamic Causal Modelling for M/EEG \label{Chap:eeg:DCM}}

\section{Introduction}
Dynamic Causal Modelling (DCM) is based on an idea initially developed for fMRI
data: The measured data are explained by a network model consisting of
a few sources, which are interacting dynamically. This network model
is inverted using a Bayesian approach, and one can make inferences
about connections between sources,
or the modulation of connections by task.

For M/EEG data, DCM
is a powerful technique for inferring about parameters that one
doesn't observe with M/EEG directly. Instead of asking 'How does the
 strength of the source in left superior temporal gyrus (STG) change
between condition A and B?', one can ask questions like 'How does the backward
connection from this left STG source to left primary auditory cortex
change between condition A and B?'. In other words, one isn't limited
to questions about source strength as estimated using a source
reconstruction approach, but can test hypotheses about what is
happening between sources, in a network.

As M/EEG data is highly resolved
in time, as compared to fMRI, the inferences are about more neurobiologically plausible
parameters. These relate more directly to the causes of the underlying
neuronal dynamics.

The key DCM for M/EEG methods paper appeared in 2006, and the
first DCM studies about mismatch negativity came out
in 2007/2008. At its heart DCM for M/EEG is a source reconstruction technique, and for
the spatial domain we use exactly the same leadfields as other
approaches. However, what makes DCM unique, is that is combines the
spatial forward model with a biologically informed temporal forward model, describing e.g. the
connectivity between sources. This critical ingredient not only makes
the source reconstruction more robust by implicitly constraining the
spatial parameters, but also allows one to infer about connectivity.
\\
\\
Our methods group is continuing to work on further
improvements and extensions to DCM. In the following, we will describe the usage of DCM for evoked
responses (both MEG and EEG), DCM for induced responses (i.e., based
on power data in the time-frequency domain), and DCM for local field
potentials (measured as steady-state responses). All three DCMs share the
same interface, as many of the parameters that need to be specified are the same for all three approaches. Therefore, we will first
describe DCM for evoked responses, and then point out where the differences to the other two DCMs lie.
\\
\\
This manual provides only a procedural guide for the practical use of DCM for
M/EEG. If you want to read more about the scientific background, the
algorithms used, or how one would typically use DCM in applications,
we recommend the following reading. The two key methods
contributions can be found in \cite{od_dcm_erp} and
\cite{sjk_dcm_erp}. Two other contributions using the model for
testing interesting hypotheses about neuronal dynamics are described
in \cite{sjk_dcm_intrinsic} and \cite{matthias_dcm_constraints}. At
the time of writing, there were also three application papers published which demonstrate
what kind of hypotheses can be tested with DCM
\cite{mg_dcm_repro,mg_feedback,marta_mmndcm}. Another good source of
background information is the recent SPM book
\cite{spm_book}, where Parts 6 and 7 cover not only DCM for M/EEG but
also related research from our group. The DCMs for induced responses
and steady-state responses are covered in \cite{cc_induced, cc_asymm} and
\cite{dcm_ssr,rm_spectralnmm,rm_massmodelspectral}. Also note that there is a DCM example file, which we put onto the webpage http://www.fil.ion.ucl.ac.uk/spm/data/eeg\_mmn/. After downloading DCMexample.mat, you can load (see below) this file using the DCM GUI, and have a look at the various options, or change some, after reading the description below.

\section{Overview}
In summary, the goal of DCM is to explain measured data (such as evoked
responses) as the output of an interacting network consisting of a several
areas, some of which receive input (i.e., the stimulus). The differences between
evoked responses, measured under different conditions, are modelled as
a modulation of selected DCM parameters, e.g. cortico-cortical
connections \cite{od_dcm_erp}. This interpretation of the evoked
response makes hypotheses about connectivity directly testable. For
example, one can ask, whether the difference between two
evoked responses can be explained by top-down modulation of early areas
\cite{mg_dcm_repro}. Importantly, because model inversion is implemented
using a Bayesian approach, one can also compute Bayesian model evidences. These
can be used to compare alternative, equally plausible, models and
decide which is the best \cite{stefan_neurodynamics}.

DCM for evoked responses takes the spatial forward model into
account. This makes DCM a spatiotemporal model of the full data set
(over channels and peri-stimulus time). Alternatively, one can
describe DCM also as a spatiotemporal source reconstruction algorithm which uses
additional temporal constraints given by neural mass dynamics and
long-range effective connectivity. This is achieved by parameterising
the lead-field, i.e., the spatial projection of source
activity to the sensors. In the current version, this can be done
using two different approaches. The first assumes that the leadfield
of each source is modelled by a single equivalent current dipole
(ECD) \cite{sjk_dcm_erp}. The second approach posits that each source
can be presented as a 'patch' of dipoles on the grey matter sheet
\cite{jean_distributed}. This spatial model is complemented
by a model of the temporal dynamics of each source. Importantly, these
dynamics not only describe how the intrinsic source dynamics evolve
over time, but also how a source reacts to external input, coming
either from subcortical areas (stimulus), or from other cortical
sources.

The GUI allows one to enter all the information necessary for specifying a
spatiotemporal model for a given data set. If you want to fit multiple
models, we recommend using a batch script. An example of such
a script (\textit{DCM\_ERP\_example}), which can be adapted to your
own data, can be found in the \textit{man/example\_scripts/} folder of
the distribution. You can run this script on example data provided by via the SPM webpage (http://www.fil.ion.ucl.ac.uk/spm/data/eeg\_mmn/). However, you first have to preprocess these data to produce an evoked response by going through the preprocessing tutorial (chapter \ref{Chap:data:mmn}) or by running the \texttt{history\_subject1.m} script in the \texttt{example\_scripts} folder.

\section{Calling DCM for ERP/ERF}
After calling \textit{spm eeg}, you see SPM's graphical user interface,
the top-left window. The button for calling the DCM-GUI is found
in the second partition from the top, on the right hand side. When
pressing the button, the GUI pops up. The GUI is partitioned into five
parts, going from the top to the bottom. The first part is about
loading and saving existing DCMs, and selecting the type of model. The second part is about selecting
data, the third is for specification of the spatial forward model, the
fourth is for specifying connectivity, and the last row of
buttons allows you to estimate parameters and view results.

You have to select the data first and specify the model in
a fixed order (data selection $>$ spatial model $>$
connectivity model). This order is necessary, because there are
dependencies among the three parts that would be hard to resolve
if the input could be entered in any order. At any time, you can switch back and forth from one part to
the next. Also, within each part, you can specify information in any
order you like.

\section{load, save, select model type}
At the top of the GUI, you can load an existing DCM or save the one
you are currently working on. In general, you can \textit{save} and
\textit{load} during model specification at any time. You can also
switch between different DCM analyses (the left menu). The default is 'ERP' which is DCM
for evoked responses described here. Currently, the other types are cross-spectral densities (CSD), induced responses (IND) and phase coupling (PHA) described later in this chapter. The menu on the right-hand side lets you choose the neuronal model.   'ERP'  is the standard model described in most of our older papers, e.g.~\cite{od_dcm_erp}.  'SEP' uses a variant of this model, however, the dynamics tend to be faster \cite{andre_sigmoid}. 'NMM' is a nonlinear neural mass model based on a first-order approximation, and  'MFM', is also nonlinear and is based on a second-order approximation. 'NMDA' is a variant of the 'NMM' model which also includes a model of NMDA receptor. 'CMC' and 'CMM'  are canonical microcircuit models ~\cite{Bastos2012} used in the more recent paper to link models of neurophysiological phenomena with canonical models of cortical processing based on the idea of predictive coding. 

\section{Data and design}
In this part, you select the data and model between-trial
effects. The data can be either event-related
potentials or fields. These data must be in the SPM-format. On the
right-hand side you can enter trial indices of the evoked responses in
this SPM-file. For example, if you want to model the
second and third evoked response contained within an SPM-file, specify
indices 2 and 3. The indices correspond to the order specified by the \textit{condlist} method (see \ref{Chap:eeg:preprocessing}). If the two evoked responses, for some reason, are in
different files, you have to merge these files first. You can do this
with the SPM preprocessing function \textit{merge}
(\textit{spm\_eeg\_merge}), see \ref{Chap:eeg:preprocessing}. You can also choose how you want to
model the experimental effects (i.e. the differences between conditions). For example, if trial \texttt{1} is the standard and trial \texttt{2} is the deviant response in an oddball paradigm, you can use the standard as the baseline and model the differences in the connections that are necessary to fit the deviant. To do that type \texttt{0 1} in the text box below trial indices. Alternatively, if you type \texttt{-1 1} then the baseline will be the average of the two conditions and the same factor will be subtracted from the baseline connection values to model the standard and added to model the deviant. The latter option is  perhaps not optimal for an oddball paradigm but might be suitable for other paradigms where there is no clear 'baseline condition'.  When you want to model three or more evoked responses, you can model the modulations of a connection strength of the second and third evoked responses as two separate experimental effects relative to the first evoked response. However, you can also choose to couple the connection strength of the first evoked response with the two gains by imposing a linear relationship on how this connection changes over trials. Then you can specify a single effect (e.g. \texttt{-1 0 1}).  This can be useful when one wants to add constraints on how connections (or other DCM parameters) change over trials. A compelling example of this can be found in \cite{marta_mmndcm}. For each experimental effect you specify, you will be able to select the connections in the model that are affected by it (see below).

Press the button \texttt{'data file'} to load the M/EEG dataset. Under 'time window (ms)' you have to enter the
peri-stimulus times which you want to model, e.g. 1 to 200 ms.

You can choose whether you want to model the mean or drifts of the data at sensor level. Select 1 for 'detrend' to just model the mean. Otherwise select the number of discrete cosine transform terms you want to use to model low-frequency drifts ($> 1$). In DCM, we use a projection of the data to a subspace to reduce the amount of data. The type of spatial projection is described in \cite{matthias_dcm_constraints}. You can select the number of modes you wish to keep. The default is 8.

You can also choose to window your data, along peri-stimulus time,
with a hanning window (radio button). This windowing will reduce the
influence of the beginning and end of the time-series.

If you are happy with your data selection, the projection and the detrending
terms, you can click on the $>$ (forward) button, which will bring you to
the next stage \textit{electromagnetic model}. From this part, you can
press the red $<$ button to get back to the data and design part.

\section{Electromagnetic model}
With the present version of DCM, you have three options for how to spatially
model your evoked responses. Either you use a single equivalent
current dipole (ECD) for each source, or you use a patch on the
cortical surface (IMG), or you don't use a spatial model at all (local field potentials (LFP)). In all three cases, you have to enter the source names (one
name in one row). For ECD and IMG, you have to specify the prior source locations (in mm in MNI
coordinates). Note that by default DCM uses  uninformative priors on
dipole orientations, but tight priors on locations. This is because tight priors on locations ensure that the posterior location will not deviate to much from its prior location. This means each dipole stays in its designated area and retains its meaning. The prior location for each
dipole can be found either by using available anatomical knowledge or
by relying on source reconstructions of comparable studies. Also note
that the prior location doesn't need to be overly exact, because the
spatial resolution of M/EEG is on a scale of several millimeters.
You can also load the prior locations from a file
('load'). You can visualize the locations of all sources when you
press 'dipoles'.

The onset-parameter determines when the stimulus,
presented at 0 ms peri-stimulus time, is assumed to activate the cortical area to which it is connected. In DCM, we usually do not model the rather small early responses, but start modelling at the first large
deflection. Because the propagation of the stimulus impulse through
the input nodes causes a delay, we found that the default value of 60
ms onset time is a good value for many evoked responses where the
first large deflection is seen around 100 ms. However, this value is a
prior, i.e., the inversion routine can adjust it. The prior mean should be chosen according to the specific responses of interest. This is because the time until the first large deflection is dependent on the paradigm or the modality you are working in, e.g. audition or vision.
You may also find that changing the onset prior has an effect on how your data are
fitted. This is because the onset time has strongly nonlinear effects
(a delay) on the data, which might cause differences in which maximum
was found at convergence, for different prior values. It is also possible to type several numbers in this box (identical or not) and then there will be several inputs whose timing can be optimized separately. These inputs can be connected to different model sources. This can be useful, for instance, for modelling a paradigm with combined auditory and visual stimulation.The 'duration (sd)' box makes it possible to vary the width of the input volley, separately for each of the inputs. This can be used to model more closely the actual input structure (e.g. a long tone or extended presentation of a visual input). By combining several inputs with different durations one can approximate an even more complex input waveform (e.g. speech).

When you want to proceed to the next model
specification stage, hit the $>$ (forward) button and proceed
to the \textit{neuronal model}.

\section{Neuronal model}
There are five (or more) matrices which you need to specify by button presses. The
first three are the connection strength parameters for the first
evoked response. There are three types of connections,
\textit{forward}, \textit{backward} and \textit{lateral}. In each of
these matrices you specify a connection \textit{from} a source area
\textit{to} a target area. For example, switching on the element
$(2,\;1)$ in the intrinsic forward connectivity matrix means that you
specify a forward connection from area 1 to 2. Some people find the
meaning of each element slightly counter-intuitive, because the
column index corresponds to the source area, and the row index to the
target area. This convention is motivated by direct correspondence between the matrices of buttons in the GUI and connectivity matrices in DCM equations and should be clear to anyone familiar with matrix multiplication.

The one or more inputs that you specified previously can go to any area and to multiple areas. You can select
the receiving areas by selecting area indices in the \textit{C input} vector.

The \textit{B} matrix contains all gain modulations of connection
strengths as set in the $A$-matrices. These modulations model the
difference between the first and the other modelled evoked
responses. For example, for two evoked responses, DCM explains the
first response by using the $A$-matrix only. The 2nd response is
modelled by modulating these connections by the weights in
the $B$-matrix.

\section{Estimation}
When you are finished with model specification, you can hit the
\textit{estimate} button in the lower left corner. If this is the
first estimation and you have not tried any other source
reconstructions with this file, DCM will build a spatial forward
model. You can use the template head model for quick results. DCM will now
estimate model parameters. You can follow the estimation process by
observing the model fit in the output window. In the matlab command
window, you will see each iteration printed out with
expected-maximization iteration number, free energy $F$, and the
predicted and actual change of $F$ following each iteration
step. At convergence, DCM saves the results in a DCM file, by default
named 'DCM\_ERP.mat'. You can save to a different name, eg. if you are estimating multiple models, by pressing 'save' at
the top of the GUI and writing to a different name.

\section{Results}
After estimation is finished, you can assess the results
by choosing from the pull-down menu at the bottom (middle).

With \textit{ERPs (mode)} you can plot, for each mode, the data for
both evoked responses, and the model fit.

When you select \textit{ERPs (sources)}, the dynamics of each area are
plotted. The activity of the pyramidal cells (which is the
reconstructed source activity) are plotted in solid
lines, and the activity of the two interneuron populations are plotted
as dotted lines.

The option \textit{coupling (A)} will take you to a summary about the
posterior distributions of the connections in the $A$-matrix. In the
upper row, you see the posterior means for all intrinsic
connectivities. As above, element $(i,\; j)$ corresponds to a
connection from area $j$ to $i$. In the lower row, you'll find, for
each connection, the probability that its posterior mean is different
from the prior mean, taking into account the posterior variance.

With the option \textit{coupling(B)} you can access the posterior
means for the gain modulations of the intrinsic connectivities and
the probability that they are unequal to the prior means. If you specified several experimental effects, you will be asked which of them you want to look at.

With \textit{coupling(C)} you see a summary of the posterior
distribution for the strength of the input into the input receiving
area. On the left hand side, DCM plots the posterior means for each
area. On the right hand side, you can see the corresponding
probabilities.


The option \textit{Input} shows you the estimated input function. As
described by \cite{od_dcm_erp}, this is a gamma function with the
addition of low-frequency terms.

With \textit{Response}, you can plot the selected data, i.e. the data,
selected by the spatial modes, but back-projected into sensor space.

With \textit{Response (image)}, you see the same as under Results but
plotted as an image in grey-scale.

And finally, with the option \textit{Dipoles}, DCM displays an
overlay of each dipole on an MRI template using the posterior means of
its 3 orientation and 3 location parameters. This makes sense only if
you have selected an ECD model under \textit{electromagnetic model}.

Before estimation, when you press the button 'Initialise' you can
assign parameter values as initial starting points for the free-energy
gradient ascent scheme. These values are taken from another already
estimated DCM, which you have to select.

The button \textit{BMS} allows you do Bayesian model comparison of
multiple models. It will open the SPM batch tool for model selection. Specify a directory to write the output file to.  For the ``Inference method'' you can choose between ``Fixed effects'' and ``Random effects'' (see \cite{klaas_bms} for additional explanations). Choose ``Fixed effects'' if you are not sure. Then click on ``Data'' and in the box below click on ``New: Subject''. Click on ``Subject'' and in the box below on ``New: Session''. Click on models and in the selection window that comes up select the DCM mat files for all the models (remember the order in which you select the files as this is necessary for interpretating the results). Then run the model comparison by pressing the green ``Run'' button. You will see, at the top, a bar plot of the log-model evidences for all models. At the bottom, you will see the probability, for each model, that it produced the data. By convention, a model can be said to be the best among a selection of other models, with strong evidence, if its log-model evidence exceeds all other log-model evidences by at least 3.

\section{Cross-spectral densities}

\subsection{Model specification}
DCM for cross-spectral densities can be applied to M/EEG or intracranial data.

The top panel of the DCM for ERP window allows you to toggle through available analysis methods. On the top left drop-down menu, select 'CSD'. The second drop-down menu in the right of the top-panel allows you to specify whether the analysis should be performed using a model which is linear in the states, for this you can choose ERP or CMC. Alternatively you may use a conductance based model, which is non-linear in the states by choosing, 'NMM', 'MFM' or 'NMDA'. (see \cite{andre_sigmoid} for a description of the differences).

The steady state (frequency) response is generated automatically from the time domain recordings. The time duration of the frequency response is entered in the second panel in the time-window.  The options for detrending allow you to remove either 1st, 2nd, 3rd or 4th order polynomial drifts from channel data. In the subsampling option you may choose to downsample the data before constructing the frequency response. The number of modes specifies how many components from the leadfield are present in channel data. The specification of between trial effects and design matrix entry is the same as for the case of ERPs, described above.

\subsection{The Lead-Field}
The cross-spectral density is a description of the dependencies among the observed outputs of these neuronal sources. To achieve this frequency domain description we must first specify the likely sources and their location. If LFP data are used then only source names are required. This information is added in the third panel by selecting 'LFP'. Alternatively, x,y,z coordinates are specified for ECD or IMG solutions.


\subsection{Connections}
The bottom panel then allows you to specify the connections between sources and whether these sources can change from trial type to trial type.

On the first row, three connection types may be specified between the areas. For NMM and MFM options these are Excitatory, Inhibitory or Mixed excitatory and inhibitory connections. When using the ERP option the user will specify if connections are 'Forward', 'Backward' or 'Lateral'. To specify a connection, switch on the particular connection matrix entry. For example to specify an Inhibitory connection from source 3 to source 1, turn on the 'Inhib' entry at position (3,1).

On this row the inputs are also specified. These are where external experimental inputs enter the network.

The matrix on the next row allows the user to select which of the connections specified above can change across trial types. For example in a network of two sources with two mixed connections (1,2) and (2,1), you may wish to allow only one of these to change depending on experimental context. In this case, if you wanted the mixed connection from source 2 to source 1 to change depending on trial type, then select entry (2,1) in this final connection matrix.


\subsection{Cross Spectral Densities}
The final selection concerns what frequencies you wish to model. These could be part of a broad frequency range e.g. like the default 4 - 48 Hz, or you could enter a narrow band e.g. 8 to 12 Hz, will model the alpha band in 1Hz increments.

Once you hit the 'invert DCM' option the cross spectral densities are computed automatically (using the spectral-toolbox). The data for inversion includes the auto-spectra and cross-spectra between channels or between channel modes. This is computed using a multivariate autoregressive model, which can accurately measure periodicities in the time-domain data. Overall the spectra are then presented as an upper-triangular, s x s matrix, with auto-spectra on the main diagonal and cross-spectra in the off-diagonal terms.

\subsection{Output and Results}
The results menu provides several data estimates. By examining the 'spectral data', you will be able to see observed spectra in the matrix format described above. Selecting 'Cross-spectral density' gives both observed and predicted responses. To examine the connectivity estimates you can select the 'coupling (A)' results option, or for the modulatory parameters, the 'coupling (B)' option. Also you can examine the input strength at each source by selecting the 'coupling (C)' option, as in DCM for ERPs. The option 'trial-specific effects' shows the change in connectivity parameter estimates (from B) from trial to trial relative to the baseline connection (from A). To examine the spectral input to these sources choose the 'Input' option; this should look like a mixture of white and pink noise. Finally the 'dipoles' option allows visualisation of the a posteriori position and orientation of all dipoles in your model.

\section{Induced responses}
DCM for induced responses aims to model coupling within and between frequencies that are associated with linear and non-linear mechanisms respectively. The procedure to do this is similar to that for DCM for ERP/ERF. In the following, we will just point out the differences in how to specify models in the GUI. Before using the technique, we recommend reading about the principles behind DCM for induced responses \cite{cc_induced}.

\subsection{Data}
The data to be modelled must be single trial, epoched data. We will model the entire spectra, including both the evoked (phase-locked to the stimulus) and induced (non-phase-locked to the stimulus) components.

\subsection{Electromagnetic model}
Currently, DCM for induced responses uses only the ECD method to capture the data features. Note that a difference to DCM for evoked responses is that the parameters of the spatial model are not optimized. This means that DCM for induced responses will project the data into source space using the spatial locations provided by you.


\subsection{Neuronal model}
This is where you specify the connection architecture. Note that in DCM for induced responses, the $A$-matrix encodes the linear and nonlinear coupling strength between sources.


\subsection{Wavelet transform}
This function can be called below the connectivity buttons and allows one to transfer data into the time-frequency domain using a Morlet Wavelet transform as part of the feature extraction.  There are two parameters: The frequency window defines the desired frequency band and the wavelet number specifies the temporal-frequency resolution. We recommend values greater than 5 to obtain a stable estimation.

\subsection{Results}

\subsubsection{Frequency modes}
This will display the frequency modes, identified using singular value decomposition of spectral dynamics in source space (over time and sources).

\subsubsection{Time-Frequency}
This will display the observed time-frequency power data for all pre-specified sources (upper panel) and the fitted data (lower panel).

\subsubsection{Coupling (A-Hz)}
This will display the coupling matrices representing the coupling strength from source to target frequencies.

\section{Phase-coupled responses}

DCM for phase-coupled responses is based on a weakly coupled oscillator model of neuronal interactions. 

\subsection{Data}
The data to be modeled must be multiple trial, epoched data. Multiple trials are required so that the full state-space of phase differences can be explored. This is achieved with multiple trials as each trial is likely to contain different initial relative phase offsets. Information about different trial types is entered as it is with DCM for ERP ie. using a design matrix. DCM for phase coupling is intended to model dynamic transitions toward synchronization states. As these transitions are short it is advisable to use short time windows of data to model and the higher the frequency of the oscillations you are interested in, the shorter this time window should be. DCM for phase coupling will probably run into memory problems if using long time windows or large numbers of trials.

\subsection{Electromagnetic model}
Currently, DCM for phase-coupled responses will work with either ECD or LFP data.  Note that a difference to DCM for evoked responses is that the parameters of the spatial model are not optimized. This means that DCM for phase-coupled responses will project the data into source space using the spatial locations you provide.

\subsection{Neuronal model}
This is where you specify the connection architecture for the weakly coupled oscillator model. If using the GUI, the Phase Interaction Functions are given by $a_{ij} sin (\phi_i-\phi_j)$ where $a_{ij}$ are the connection weights that appear in the $A$-matrix and $\phi_i$ and $\phi_j$ are the phases in regions $i$ and $j$. DCM for phase coupling can also be run from a MATLAB script. This provides greater flexibility in that the Phase Interaction Functions can be approximated using arbitrary order Fourier series. Have a look in the \texttt{example\_scripts} to see how.

\subsection{Hilbert transform}
Pressing this button does two things. First, source data are bandpass filtered into the specified range. Second, a Hilbert transform is applied from which time series of phase variables are obtained. 

\subsection{Results}

\subsubsection{Sin(Data) - Region i}
This plots the sin of the data (ie. sin of phase variable) and the corresponding model fit for the $i$th region.

\subsubsection{Coupling (A),(B)}
This will display the intrinsic and modulatory coupling matrices. The $i,j$th entry in $A$ specifies how quickly region $i$ changes its phase to align with region $j$. The corresponding entry in $B$ shows how these values are changed by experimental manipulation. 
