1   This demos SPM5 (with cbu_updates) for distributed source 
   2   =========================================================
   3          localisation of CBU Neuromag MEG data
   4 	 =====================================
   5 
   6 Rik Henson, Sep 07
   7 ==================
   8 
   9 There are two demo localisations, on the same raw data.
  10 
  11 Demo 1 describes how to examine Left-hand vs Right-hand evoked responses (~0ms 
  12 relative to keypress), using the GUI
  13 
  14 Demo 2 describes how to examine Faces vs Scrambled faces evoked responses (or 
  15 evoked and induced across all trials), using command-line (ie batch scripts)
  16 
  17 
  18 FIRST
  19 =====
  20 
  21 You will need to copy the structural to a subdirectory, eg "sMRI", including a 
  22 mat file containing at least three reference points (fiducials) in MRI space (in 
  23 this example, nasion and left and right periaricular points), or you can create 
  24 this yourself, by simply displaying the MRI in SPM, creating a matrix of 3 rows 
  25 (nasion, left ear, right ear) and 3 columns (their x,y,z, in mm), and saving as 
  26 a mat file.
  27 
  28 You will also need the text files that recode the events of interest (there were 
  29 actually at least 8 different codes in the original FIF file, with rather odd 
  30 codings, corresponding to four different types of stimuli and four different 
  31 types of keypress). The website should provide these as "hands.txt" (for Demo 1) 
  32 and "faces.txt" (for Demo 2).
  33 
  34 Both demos require that the RAW FIF data is first run through MaxFilter 
  35 (including movement estimation every 200ms and compensation, HPI signal removal, 
  36 conversion to short, and downsampling by a factor of 4):
  37 
  38 /neuro/bin/util/maxfilter -gui -f 
  39 /megdata/cbu/henson/henson_rik/070412/faces_raw.fif -o raw.fif -ctc 
  40 /neuro/databases/ctc/ct_spars'bcdnste.fif -cal /neuro/databases/sss/sss_cal.dat 
  41 -format short -movecomp -hpistep 200 -hpisubt amp -hp raw_movement.pos -ds 4
  42 
  43 
  44 1 LEFT vs RIGHT HAND MEAN EVOKED RESPONSES (GUI)
  45 ================================================
  46 
  47 
  48 1.1 Type "spm 5" on a linux box, and toggle SPM to EEG mode (or type "spm 5 
  49 eeg").
  50 
  51 
  52 1.2 Press "Convert", select "FIF", select "raw.fif", select "FIF306_setup.mat" 
  53 (a layout file for display)...
  54 
  55 ...the Graphics window should now show the position and three (orthogonal) 
  56 orientations of each sensor (magnetometer and planar gradiometers) in 3D MEG 
  57 space, including sensor name (use cursor to rotate). The red dots are the 
  58 Polhemus digitised points...
  59 
  60 ...press return to select the default channel for the trigger codes ("STI101"), 
  61 press return to convert the whole raw data (ie start and end of timeseries)...
  62 
  63 ... then select the channel "EEG062" (second one) as the vertical EOG channel 
  64 (the other one will be made horizontal EOG; this depends how the electrodes were 
  65 attached)...
  66 
  67 ...this will now produce "raw.mat" and "raw.dat" which are the raw data in SPM 
  68 format. The mat file is a header file (you can load and peruse the structure "D" 
  69 inside); the dat file contains all the data (which will be memory-mapped later; 
  70 see spm_eeg_ldata). 
  71 
  72 There will also be some output in the Matlab window describing what was read, 
  73 and including the rescaling of the data. The gradiometers are downscaled 
  74 relative to the magnetometers (this can be overriden if you call 
  75 spm_eeg_rdata_FIF from the command line). This scaling is because the 
  76 gradiometers measure a gradient (T/m) and the magnetometers measure total flux 
  77 (T). They are therefore not commensurate. One can scale by the baseline SNR, 
  78 which NeuroMag report is typically 1:20 for grads:mags. I decided to scale the 
  79 gradiometers by 0.0595 (rather than 0.05), which is a function of the baseline 
  80 between the centres of two gradiometer coils, which is 16.8cm (if you want to 
  81 change this scaling, it is a hidden function call in spm_eeg_rdata_FIF; see 
  82 batch script that accompanies this README for an example). This scaling seems to 
  83 make the magnetometer and gradiometer data have similar amplitude range at 
  84 least. Alternative scalings are described in the NeuroMag manuals.
  85 
  86 
  87 1.3 Press "Epoching", select "raw.mat" file, enter "-300" for start of epoch and 
  88 "300" for end...
  89 
  90 ...then type "-2304 -1280 -768 -512" for the event-types to epoch. These strange 
  91 numbers simply correspond to the trigger codes for the four types of finger 
  92 press (left middle, left index, right index and right middle) - as distinct from 
  93 the codes for the visual stimuli (see Demo 2). (They are strange simply because 
  94 of a negative offset that was not set properly in the acquisition software for 
  95 this particular dataset)...
  96 
  97 ...then type "yes" to read a new event-list. If you press "no", it will epoch 
  98 the four trigger codes above, and keep those codes. However, for this Demo, we 
  99 want to collapse finger and just compare left vs right hand. This has been done 
 100 already in the text file "hands.txt". This file simply contains 1s and 2s for 
 101 each of the key press triggers (172 in total), and will recode each event into 1 
 102 (left hand) or 2 (right hand). This illustrates the type of procedure you will 
 103 need if you want to further conditionalise your original trigger codes on, for 
 104 example, subsequent behavioural analysis (eg correct vs incorrect)...
 105 
 106 ...It will prompt you to enter 172 new codes, but because the SPM input is 
 107 evaluated, you can type "spm_load" and then select the file "hands.txt" 
 108 described above. 
 109 
 110 This will produce two new files (e_raw.mat and e_raw.dat), which are the epoched 
 111 data.
 112 
 113 
 114 1.4 Press "Display", choose the "M/EEG" option, and select the "e_raw.mat" 
 115 file...
 116 
 117 ...the graphics window will now show all the channels. The MEG channels are 
 118 grouped into triplets with the magnetometer on top and the gradiometers below. 
 119 The two EOG channels are in the top left and right. If you left-click on the top 
 120 left channel, a new blow-up window will appear with the VEOG data (for the first 
 121 epoch). (Make sure the Matlab "Rotate3D" cursor is not still set, from the 
 122 Matlab "Tools" menu). You should see a bump around -200ms which is probably a 
 123 blink. This trial will be removed by artifact rejection next... 
 124 
 125 (You can explore the data further, eg other trials, scaling, topographies, etc, 
 126 but this is slow, and probably best done after the next step.)
 127 
 128 (Note also that if you select the "M/EEG fast" Display option instead, the 
 129 plotting will be slightly faster, and you can now use the left-mouse to "zoom 
 130 in" on channels, and the right-mouse to plot them).
 131 
 132 
 133 1.5 Press "Artefacts", select the e_raw.mat file, press "no" to read own list 
 134 (e.g, a list of rejected trials from another package; note you can also indicate 
 135 rejects manually through visual inspection of each trial via the Display 
 136 option)...
 137 
 138 ...press "no" for robust averaging (this is essentially an average weighted by 
 139 the typicality of each epoch, see Kilner et al, HBM06, but is very slow), then 
 140 press "yes" to Thresholding...
 141 
 142 You can now enter a single number to threshold all channels by that amount (in 
 143 absolute terms), or as many numbers as there are channels (306+2=308) if you 
 144 want a channel-specific threshold. There is no right or wrong answer here. Since 
 145 one rarely has a good idea what a max or min magnetic deflection is (in fT), you 
 146 can turn off thresholding of the MEG channels by giving a threshold of infinity. 
 147 Past experience for EOG suggests that blinks often exceed a vertical deflection 
 148 of 200uV. To implement this thresholding via EOG channels only, type
 149 
 150     "[ones(1,306)*Inf 200 200]"
 151 
 152 which simply says the above in Matlab.
 153 
 154 The Matlab window will now report that there were no "Bad" channels (channels 
 155 with more than 20% of trials above threshold), but 19 rejected trials (i.e, with 
 156 EOG deflections >200). Note we are currently working on more sophisticated 
 157 methods, eg for blink correction using ICA (called from EEGLAB).
 158 
 159 
 160 1.6 If you now select "Display: M/EEG" and select the "ae_raw.mat" file, you can 
 161 see which trials were rejected by the menu on the bottom right of the graphics 
 162 window. Inspection of these will show that they are all blinks picked up in the 
 163 VEOG channel around the time that the button was pressed (quite a common 
 164 occurrence: there were very few blinks in the epoch around the visual triggers 
 165 in Demo 2 below because the subject was instructed not to blink when the visual 
 166 stimulus was displayed).
 167 
 168 1.7 Now we want to create an average (excluding those trials marked as rejects 
 169 from 1.6). Press "Averaging", select the "ae_raw.mat" file. The Matlab window 
 170 will report that 96 trials were of type 1 (left hand0 and 57 of type 2 (right). 
 171 The new file created is called "mae_raw.mat" (m for mean).
 172 
 173 
 174 1.8 To aid visual inspection, you can now filter the data (though this is not 
 175 always advisable for analysis, depending on the frequency cutoff relative to 
 176 epoch length; filtering raw data is better to avoid end-effects). You could 
 177 press "Filter", select "low-pass", enter "40" as the cutoff, and you will get a 
 178 file "fmae_raw.mat" containing smoothed data. (This uses a high-order, optimised 
 179 Butterworth digital filter.)
 180 
 181 
 182 1.9 You now select "Display: M/EEG" and select the "fmae_raw.mat" file to see 
 183 the smoothed averages. You speed things up by selecting only a subset of 
 184 channels to display. Press "channels" in the Graphics window, and then "Deselect 
 185 all". Then select only the first 102 channels from the menu on the right using 
 186 Shift-Left-mouse. The first 102 channels are the magnetometers (also indicated 
 187 by the fact that their names end in "1" rather than "2" or "3", for the 
 188 gradiometers). You can also select the two EOG channels from the bottom of the 
 189 menu using Ctrl-Left-mouse.
 190 
 191 If you now press shift and left click trial-type 2, you will see a red line 
 192 added to the channels. If you click on one of the more central (magnetometer) 
 193 channels, you can see the left vs right hand differences around 0ms (key press).
 194 
 195 If you press "Topography", type 0ms (rather than 100), and select "2D", you will 
 196 see the evoked responses for the first trial-type. You should see red/blue 
 197 maxima on the right, corresponding to the left hand. If you first select trial-
 198 type 2, then press Topography, you will see two larger maxima, more on the left 
 199 (for right hand; the subject was right-handed).
 200 
 201 If you move the slider left and right, the topography will change over time 
 202 (precise timepoint shown in Matlab window).
 203 
 204 
 205 1.10 LOCALISATION. We are now going to localise the left vs right hand responses 
 206 in 3D. Press "3D reconstruction"...
 207 
 208 (The button "2D interpolation" would allow you to write out 2D images in sensors 
 209 space, if you want to do statistics in sensor space, using the usual machinery 
 210 of SPM, including RFT correction for multiple comparisons. However, topographies 
 211 are difficult to interpret for our MEG data with its three different sensor-
 212 types (unlike for EEG data for example), since they give a vector rather than a 
 213 scalar field. So any "topographies" displayed as above using all the sensors are 
 214 fairly meaningless.)
 215 
 216 ...in the new window, press "Load" and select the "mae_raw.mat" file (or fmae 
 217 one, if you prefer). Enter any comment you want...
 218 
 219 ...then start the first of four stages:
 220 
 221 Stage 1: Segment and normalise the MRI for this subject. Segmentation allows 
 222 creation of boundaries of scalp, inner skull and cortex. Normalisation allows 
 223 the canonical cortical mesh to be inverse-normalised to match this specific 
 224 subject. (See papers on WIKI)...
 225 
 226 (If you have no MRI, you can use a canonical MRI by pressing the "template" 
 227 button, with associated greater potential for errror)
 228 
 229 ...select the "smri.img" structural MRI (available from WIKI), then select 
 230 "7200" for the number of vertices in the mesh (ie maximum resolution)
 231 
 232 This call's SPM's normal segmentation, which will take some time. Put your feet 
 233 up; you deserve it. Having got this far... (you can watch some output in the 
 234 Matlab window).
 235 
 236 When it has finished, an image will appear in the Graphics window showing the 
 237 cortical, inner skull and scalp meshes.
 238 
 239 Stage 2: Coregister the MEG data to the MRI space. This is done in two stages: a 
 240 coarse co-registration by fitting the fiducials in MEG and MRI space (the latter 
 241 demarked by hand; see above), and a finer co-registration by fitting the 
 242 digitised headshape in MEG space to the MRI scalp mesh created during Stage 1.
 243 
 244 Press "Co-register", select the file "smri_fids.mat" (winzipped with the 
 245 structural on the WIKI), and watch the Graphics window. When it has finished, 
 246 you will see in the top panel the MRI meshes surrounded by the sensors (in 
 247 green; magnetometers and gradiometers collapsed per chip), plus the fiducials 
 248 according to MRI and MEG (purple and blue markers). The fiducial fit is not 
 249 excellent because the fiducials were only approximated by hand on the MRI, but 
 250 this does not matter because the precise fit is provided from the many more 
 251 points digitized on the headshape (and matched to MRI scalp). You can rotate the 
 252 image to see more. The bottom panel just shows the sensor names.
 253 
 254 Stage 3: Construct a forward model. This vital step creates the leadfield matrix 
 255 that expresses how each of the (7200) possible sources contributes to the (306) 
 256 sensors. Note that the sources are constrained to be normal to the cortical 
 257 surface; only their amplitude is estimated. This requires Maxwell's equations, 
 258 and for magnetic fields, can be approximated reasonably well by a single-sphere 
 259 model (Brainstorm code for a Boundary Element Model (BEM) is available, but not 
 260 yet integrated). This sphere is fit to the inner skull (show dynamically in the 
 261 Graphics window). All the matters from this fit is the centre of the sphere (and 
 262 the position and orientation of the sensors, which is already encoded in the 
 263 *.mat file as derived from the *.fif file).
 264 
 265 Stage 4: We now want to invert the forward model, given the MEG data. Press 
 266 "Invert", select "Classical" ("DCM" is for new methods for ECD), select "yes" to 
 267 localise both trial-types (left and right hand), and select "MSP"...
 268 
 269 ...MSP stands for "Multiple Sparse Priors", which is a novel approach unique to 
 270 SPM. It resembles weighted minimum L2-norm, but with multiple rather than one 
 271 prior. More details can be found in the paper on the WIKI. It has been shown to 
 272 be superior to the other two options: "MNM" (which is basic minimum norm, 
 273 minimising "energy" of solution) and "COH" (for "coherent", which maximises 
 274 spatial smoothness of solution, like LORETA but in the context geodesic distance 
 275 within a mesh)...
 276 
 277 ...select "Standard" to use standard parameter values (or "Custom to hand-
 278 craft).
 279 
 280 The Matlab window will show some details of the fitting (inversion via ReML). 
 281 The Graphics window will then show the results for condition 1 (left hand). It 
 282 shows by default the maximum in 3D space and time. The lines are the estimated 
 283 source amplitudes at the coordinate described (43, -36, 45; somewhere in right 
 284 motor or somatosensory cortex?). The red line is the condition currently 
 285 selected in the reconstruction window (ie condition 1; left hand); the grey 
 286 lines are any other conditions. The peak is at 28ms, and 87% of the total 
 287 variance has been explained (the log-evidence is actually a more useful, though 
 288 relative, measure, since one can explain all the variance in such 
 289 underdetermined problems, unless one explicitly models the noise. And yes, it 
 290 should be negative.). The bottom panel shows the MIP of posterior probabilities 
 291 of reliable sources.
 292 
 293 If you press "movie", and enter times from -300 to 300 (whole epoch), you can 
 294 see the estimated source amplitudes evolve over time. (Note right-lateralisation 
 295 around 0ms, but evidence of subsequent left and bilateral activity. Ask a motor 
 296 person.)
 297 
 298 If you enter 0 (or 200) into the time box and press MIP, you can see the MIP for 
 299 that time.
 300 
 301 If you press the "Condition 1" button it will toggle to the other conditions (ie 
 302 Condition 2; right hand), and then press MIP again (for 0ms), you will see more 
 303 left-lateralisation for condition 2 (right hand), though somewhat less focal, 
 304 and possibly several generators.
 305 
 306 If you press "Render" then you can view the data in lots more fancy ways. Mostly 
 307 self-explanatory, though as above, ignore the scalar sensor topographies which 
 308 are incorrect for our vector-like sensor data. Fancy, but perhaps not as 
 309 relevant as...
 310 
 311 ...To get an estimate of the sources over a particular timewindow and frequency 
 312 range, press "Window". This creates a time-freq contrast, eg by entering [0 50] 
 313 for the time and 0 for the frequency band, you will get a (Hanning windowed) 
 314 average source map for all frequencies....
 315 
 316 ...This map can now be written as a 3D Nifti image in MNI space by pressing the 
 317 "Image" button (and selecting a spatial smoothing of, eg, 12mm). The blobs in 
 318 this image will be displayed together with the structural MRI. This image can be 
 319 entered into (2nd-level) analyses just like for fMRI data, eg, by repeating this 
 320 process for multiple conditions/subjects to make population-level inferences. 
 321 Such 3D group-analyses are less easy in other packages.
 322 
 323 
 324 1.11 You can obviously explore other options, eg MNM or COH localisations (you 
 325 can avoid repeating Stages 1-3 above simply by pressing "Save" in the 
 326 reconstruction window, and then "New" (new localisation of these data). You will 
 327 notice for example that the MNM solution is more dispersed and superficial 
 328 (particular for the deeper, fusiform sources in Demo 2 below). Or you can 
 329 explore the relative importance of magnetometers vs gradiometers by artificially 
 330 switching off one or t'other by selecting them to be "Bad" (this has to be done 
 331 via the command line, as illustrated in the batch script that accompanies this 
 332 demo).
 333 
 334 
 335 1.12 Note that we have localised evoked responses (from trial-average). SPM can 
 336 also localise all responses (evoked and induced) simply by localising the 
 337 average covariance across trials (see paper on WIKI), rather than covariance of 
 338 the average. To do this, simply repeat the above localisation process, but 
 339 select the "ae_raw.mat" file (ie all epochs, pre-averaging) rather than the 
 340 "mae_raw.mat" file.
 341 
 342 
 343 2 FACES vs SCRAMBLED EVOKED RESPONSES (via batch script)
 344 ========================================================
 345 
 346 For this Demo, you can use the batch script provided "spm5_cbu_meg_batch.m" 
 347 (batching for Demo 1 is also included). Just cut and paste relevant sections. 
 348 
 349 The batch does 7 localisations (note: spm_eeg_weight_epochs used to create 
 350 contrasts of MEG data, eg differential ERFs): 
 351 1) differential ERF between faces and scrambled faces (all channels, whole 
 352 epoch)
 353 2) common (average) ERF (all channels, whole epoch)
 354 3) faces and scrambled faces separately (all channels, whole epoch)
 355 4) faces and scrambled faces separately (mags only, whole epoch)
 356 5) faces and scrambled faces separately (grads only, whole epoch)
 357 6) differential ERF between faces and scrambled faces (all channels, only 150-
 358 200 (M170)),
 359 7) faces and scrambled faces separately (all channels, only 150-200 (M170)),
 360 
 361 Things to note include: 
 362 a) in localisation 3 (all channels, faces and scrambled separately), there is 
 363 greater fusiform and temporal pole activity for faces than scrambled around 
 364 160ms
 365 b) in localisation 3 there is similar more posterior activity around 100ms 
 366 (localisations 3 and 2)
 367 c) localisation 1 of the differential ERF is not as clear (in this case), with 
 368 greatest early visual differences 
 369 d) gradiometers pick up posterior fusiform ok, but not temporal poles, and 
 370 magnetometers vice versa (cf localisations 4 and 5)
 371 e) reducing the timewindow localised from the whole epoch to 150:200 is 
 372 reasonable, but introduces possible ghost sources in superior regions, 
 373 particularly for differential ERF. This is likley to reflect the fewer data. If 
 374 the same reduced timewindow is used together with all trials (ie using the 
 375 e_faces.mat file), the localisations are much better - largely restricted to 
 376 fusiform - reflecting the better estimate of the data covariance. Now you are an 
 377 expert, you can attempt the latter (via a script)....
 378 
 379 ;-)
 380