1. Anatomical Data Acquisition and Processing

1. Acquire one magnetization-prepared rapid gradient echo (MPRAGE) MRI scan of the subject. This may take 5-10 min depending on which specific scanning protocol is used.
2. Acquire two additional fast low-angle shot (FLASH) MRI scans (flip angles = 5° and 30°) if EEG data are used for inverse imaging analysis, as FLASH sequences provide different tissue contrast from the standard MPRAGE sequences ¹.
3. Use FreeSurfer software (see Table) ^{2, 3} to reconstruct the cortical surface and to set up individual M/EEG dipole source space.
   1. This source space is constrained to the grey/white matter boundary segmented from the MPRAGE scan. Each hemisphere contains roughly 100,000 potential vertices, spaced at ~1 mm. For dipole amplitude estimation (see below), use a grid spacing of 7 mm, which yields ~3,000 dipoles per hemisphere.
4. Reconstruct the skin, outer skull and inner skull surfaces from the MPRAGE and FLASH images using MNE (see Table) and FreeSurfer. Use these surfaces to generate a three-layer boundary element model (BEM).

2. M/EEG Data Acquisition

1. Prepare subject for M/EEG recording.
   1. Refer to Liu et al ⁴ for details of electrooculogram and reference electrode preparation as well as digitization of the subject's fiducial landmarks, head-position indicator coils (HPI) and EEG electrodes.
2. Once subject is seated in the MEG, measure the head position using the HPI coils.
3. Start recording. Begin presentation of auditory and visual stimuli.
   1. Many hardware and software solutions are available to perform stimulus presentation (e.g., Presentation, E-Prime). We use a Tucker-Davis Technologies RZ6 for auditory stimulus presentation and trigger stamping, with Psychtoolbox ⁵ for visual stimulus presentation, both controlled by MATLAB. Testing the auditory and visual latencies using a microphone and photodiode attached to the screen, and subsequently ensuring there is no observable jitter (which may necessitate setting the presentation projector to its native resolution) prior to the experiment helps ensure timing integrity.
4. Subject responds to auditory and visual stimuli via an optical button box while performing behavioral task.
5. Save all stimuli, experimental parameters and data files for off-line analysis.

3. M/EEG Co-registration with MRI Scan and Data Processing

1. Using the MNE software, load digitizer data and subject's reconstructed MRI head model. Pick fiducial landmarks to initiate co-registration process and proceed to use automatic alignment procedure to complete coordinate transformation (Figure 2).
2. To relate the location of each dipole in the source space with the location of each sensor, combine recorded HPI data (see 2.2) to compute a forward solution with the three-layer BEM (see 1.4)
3. Inspect all recorded M/EEG data and identify channels that have exceptionally high variance or are completely flat. Set these channels as bad channels.
4. Use signal-space projection ⁶ or other noise reduction techniques (such as signal space separation ⁷) to project or separate out spatial field patterns originated from ambient environmental field contamination or other undesirable physiological signals, such as those associated with eye-blinks and cardiac artifacts (Figure 3).
   1. Apply time-domain artifact removal (e.g., removing epochs containing abnormally high amplitude signals due to spiking of a channel) and frequency-domain artifact removal (e.g., band-notch filtering at 50 or 60 Hz line-frequency) to further increase signal-to-noise ratio.
5. Identify a baseline period in which the subject was not performing any task (e.g., 200 ms period prior to the onset of each trial). Generate an average of these baseline epochs in order to obtain a noise estimate (also known as the covariance matrix).
6. Identify epochs of interest (e.g., only collecting epochs with correct behavioral responses) and define conditions for experimental contrasts (e.g., epochs associated with subject having switched their auditory attention to the opposite hemifield as originally cued - "Switch" condition - versus subject maintaining attention on the original hemifield - "Hold" condition). Generate an average response for each of the condition defined.
   1. These averages can be baseline-corrected or not depending on experimental parameters (see ⁸); data shown here are baseline-corrected.
7. Combine the covariance matrix (3.5) and the computed forward solution (3.2) to obtain a distributed cortically constrained minimum-norm inverse operator that relates the sensor measurements to dipole current estimates in the source space.
   1. You can either approximately constrain or fix the dipole orientation to the cortical normal direction ⁹.
8. Generate a "brain movie" of the distributed dipole estimate (i.e., the current estimate at each dipole location in the source space in time) for each experimental condition (Figure 4).
   1. Depending on the temporal characteristics of your experimental design, you can bin your data in time by averaging current estimates using non-overlapping temporal windows.

4. Statistical Inference Based on a Common Surface-based Coordinate System

1. Morph the "brain movies" for each subject onto a common (average) cortical space based on a surface-based coordinate system that optimally aligns individual sulcal-gyral patterns ³. This allows us to compare or average cortical activities across all subjects. (Figure 5).
2. There are many different statistical inference approaches. We will highlight three possible approaches here. Approaches that are not implemented in software package can be written using custom software-in our examples we use MATLAB to perform the non-parametric spatio-temporal clustering permutation test. Despite the high dimensionality (Space x Time x Subjects) of these data, all of these approaches can be performed using standard modern desktop computer hardware in seconds (ROI; 4.3 approach) to hours (non-parametric clustering; 4.5).
3. Region-of-interest (ROI) approach
   1. You can define the ROI anatomically (e.g., defined by automatic parcellation algorithm ¹) and / or functionally (e.g., by recording a functional localizing task, such as a Go/No Go saccade task to identify the oculomotor regions).
   2. You can further constrain your analysis to a specific time-of-interest that is appropriate to your experimental paradigm (e.g., a time period immediately before and after the onset of the sound stimuli). You can also use other statistical inference associated with time-series analysis.
4. Whole-brain Bonferroni or False-Discovery-Rate (FDR) correction
   1. Employ Bonferroni or FDR correction if you require whole-brain, all-time analysis.
   2. Generate a statistical map at each dipole location and each time point using appropriate test statistics, such as a t-test or within-subjects ANOVA for approximately normally-distributed data. For example, z-scores from dynamic statistical parameter mapping of the MNE estimates for fixed dipole sources ¹⁰ can be used when paired with a correction for correlations in the estimates (such as the conservative Greenhouse-Geisser correction).
   3. For Bonferroni correction, obtain significant space-time points by thresholding at a significance level of 0.05 divided by the number of comparisons (number of dipoles multiplied by number of time points). For a less conservative approach, use FDR p-value correction ¹¹.
5. Non-parametric spatiotemporal clustering
   1. Use this method (based on a simple extension of ¹²) to find regions of large, consistent spatial and temporal activation while being less conservative than Bonferroni correction, and less prone to Type I statistical errors than FDR by controlling for the family-wise error rate.
      1. Because this approach uses permutation or Monte Carlo resampling techniques, it does not rely on assumptions of normality of the data, and only assumes that the condition labels are exchangeable under the null hypothesis. Although it is more computationally intensive than the previous two approaches, it can still be performed in hours on a single machine using modern desktop computer hardware.
   2. Generate a statistical map at each dipole location and each time point using the appropriate test statistics, such as a t-test.
   3. Threshold this map at a preliminary significance threshold, e.g. p < 0.05.
   4. Cluster these putative significant points based on spatio-temporal proximity, e.g. significant points within 5 ms and 5 mm geodesic distance of one another are put in the same cluster. Score each resulting cluster using hypervolume or total significance (e.g. sum of t-scores of points in the cluster).
   5. Perform a standard permutation resampling (or Monte Carlo resampling for larger datasets, e.g. number of subjects N > 10, to save on computation) test with a maximal statistic (see ¹² for permutation test examples). Briefly, for a random subset of subjects (choosing anywhere from 0 to N subjects), relabel the conditions being compared before obtaining the statistical map, perform clustering on the new statistical map, and obtain the maximal cluster score for that relabeling. Perform this procedure on new random relabelings for up to 2^N permutations to obtain a distribution of the maximal statistic; performing all 2^N possible relabelings yields the permutation test and utilizing a random subset of fewer than 2^N relabelings yields a Monte Carlo (or random) permutation test.
   6. Obtain the significance of a given original cluster (from the original labeling) by determining the proportion of time the maximal cluster sizes were greater than that of the original cluster, e.g. clusters that were larger than 95% of the maximal statistic clusters can be declared significant.
      1. For an in-depth discussion on statistical inference in MEG distributed source imaging, see ¹³.
6. The resulting data files can be visualized in many ways, including using the formats natively utilized by MNE software to store spatio-temporal cortical estimates, namely STC files. These, alongside labels that can be produced corresponding to the significant regions, can be generated using MNE toolboxes provided for MATLAB and Python.

5. Representative Results

Figure 6 shows a set of representative results using the behavioral paradigm outlined in Figure 4. Using the non-parametric spatiotemporal clustering procedure (4.5), the right FEF is found to be significant when a subject is performing a reorientation task compared to a standard task (Figure 6 left). Using the ROI approach (4.3), the time course of the right FEF is shown, along with the time period that these two conditions are significantly different.

Figure 1. Workflow for generating a "brain movie" using cortically-constrained minimum-norm dipole estimates (cf., Figure 1 in Liu et al., 2010).

Figure 2. MNE software used to facilitate EEG channels and HPI locations co-registration onto one subject's MRI co-ordinate space.

Figure 3. MEG data before and after using SSP to remove cardiac (highlighted in orange) and eye-blinks (highlighted in blue-green) artifacts and lowpass filtering to remove line-frequency. Click here to view larger figure.

Figure 4. A "brain movie" on subject's native cortical space and the timing of the audio-visual presentation (with auditory stimuli presented at 600 ms and a visual stimulus presented at -600 ms) in one experimental paradigm (Note: this will be presented as a movie in the final movie clip)

Figure 5. Comparison between a hypothetical ROI mapped on a subject's native cortical space and after morphed onto a common cortical space.

Figure 6. Representative spatio-temporal cluster and time-course associated with the two experimental conditions tested.

In order to estimate the dipole activation on the cortex from the acquired MEG/EEG data, we need to solve an inverse problem, which does not have a unique stable solution unless appropriate anatomically and physiologically sound constraints are applied. Using the anatomical constraint acquired for individual subjects using MRI and adopting the minimum-norm as our estimation criterion, we can arrive at an inverse cortical current source estimate that agrees with the sensor measurements. This approach has proved useful in studies of not only auditory processing ¹⁴ but also other domains such as visual ¹⁵ and language processing ¹⁶ .

There are many other inverse approaches. However, all these methods can be summarized into two categories: localization (e.g., equivalent current dipole modeling) or imaging (e.g., MNE, beamforming techniques). Furthermore, each inverse approach has its tradeoff (see ¹⁷ for an in-depth discussion). For example, the current estimate using the approach presented here must necessarily be distributed in space due to its minimum-norm constraint. This minimum-norm estimate approach is well suited for tasks that recruit a distributed cortical network. Mapping early responses to stimuli that evoke focal source activity, such as those in audition that are believed to be localized in and around bilateral auditory cortex (e.g., N1m and awareness related negativity ¹⁸), can also be improved by using fMRI co-constraints ¹⁴.

Spectral domain analysis, e.g., investigating the role of different cortical rhythms involved in attention, across the cortex can also be performed after using any of the aforementioned inverse techniques. Furthermore, this type of analysis can easily be extended to address questions related to functional connectivity between distinct regions in the brain.