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
Attached FilesTo refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
You are not allowed to attach a file to this page.