How to incorporate head movements in MEG analysis
Description
Changes in head position during MEG sessions may cause a significant error in the source localization. Besides, the mixture of different head positions over time adds variance to the data that is not accounted for by the experimental manipulation. Thus head movements may deteriorate statistical sensitivity when analyzing MEG on both sensor and source levels. It is therefore recommended to incorporate head movements in the offline MEG analysis, see Stolk et al., NeuroImage 2013.
Continuous head localization information is stored in HLC channels (Head Localization Channels) in CTF MEG system. An example script here shows how to read these channels in FieldTrip and estimate the amount of movement offline. Information from these channels can also be used to track the head position in real time. The data used in this example script can be obtained here
In general there are multiple ways that you can use the continuous head localization information.
 you can discard a subject or trial(s) from subsequent analysis if he/she moved too much
 you can regress out the movements from the processed data
 you can compensate the raw data for the movements
 you can correct the forward model (i.e. the leadfield) for the spatial blurring that is due to the movements
The first way of dealing with it requires that you visualize and decide on the movements. This is demonstrated in the first half of the example script. A more general description of the data, also containing information on the head movements, is made using ft_qualitycheck. At the donders, this function is run automatically on each new MEG dataset every night.
The second way of dealing with the movements means that you perform ft_timelockanalysis, ft_freqanalysis or ft_sourceanalysis with the option keeptrials=yes. This will give trial estimates of the ERF, the power or the source strength for each trial. The effect that the variable head position has on those singletrial estimates can be estimated and removed from the data using ft_regressconfound. This method has been found to significantly improve statistical sensivity following head movements, up to 30%, and is therefore demonstrated in the second half of the example script.
The third way of dealing with the movements requires that you make a spatial interpolation of the raw MEG data at each moment in time, in which you correct for the movements. In principle this could be done using the ft_megrealign function, but at this moment (May 2012) that function cannot yet deal with withinsession movements.
The fourth way of dealing with the movements is implemented in the ft_headmovement function. It is not explained in further detail on this example page.
Readingin and visualizing the head localization
Prepare configuration to define trials:
cfg = [];
cfg.dataset = 'TacStimRegressConfound.ds';
cfg.trialdef.eventtype = 'UPPT001';
cfg.trialdef.eventvalue = 4;
cfg.trialdef.prestim = 0.2;
cfg.trialdef.poststim = 0.3;
cfg.continuous = 'yes';
cfg = ft_definetrial(cfg);
Read the data with the following HLC channels:

HLC00n1 X coordinate relative to the dewar (in meters) of the nth head localization coil

HLC00n2 Y coordinate relative to the dewar (in meters) of the nth head localization coil

HLC00n3 Z coordinate relative to the dewar (in meters) of the nth head localization coil
cfg.channel = {'HLC0011','HLC0012','HLC0013', ... 'HLC0021','HLC0022','HLC0023', ... 'HLC0031','HLC0032','HLC0033'}; headpos = ft_preprocessing(cfg);
Determine the mean (per trial) circumcenter (the center of the circumscribed circle) of the three headcoils and its orientation (see subfunction at the bottom of this page)
% calculate the mean coil position per trial
ntrials = length(headpos.sampleinfo)
for t = 1:ntrials
coil1(:,t) = [mean(headpos.trial{1,t}(1,:)); mean(headpos.trial{1,t}(2,:)); mean(headpos.trial{1,t}(3,:))];
coil2(:,t) = [mean(headpos.trial{1,t}(4,:)); mean(headpos.trial{1,t}(5,:)); mean(headpos.trial{1,t}(6,:))];
coil3(:,t) = [mean(headpos.trial{1,t}(7,:)); mean(headpos.trial{1,t}(8,:)); mean(headpos.trial{1,t}(9,:))];
end
% calculate the headposition and orientation per trial (for function see bottom page)
cc = circumcenter(coil1, coil2, coil3)
Now you can plot the head position relative to the first value, and compute the maximal position change.
cc_rel = [cc  repmat(cc(:,1),1,size(cc,2))]';
% plot translations
figure();
plot(cc_rel(:,1:3)*1000) % in mm
% plot rotations
figure();
plot(cc_rel(:,4:6))
maxposchange = max(abs(cc_rel(:,1:3)*1000)) % in mm
The figure illustrates head position changes during 1hour MEG session (data used for this plot are different from those used in the example above). You may decide to exclude a subject from the subsequent analysis if the head movement exceeds a certain threshold.
Regressing out headposition confounds
MEG experiments typically involve repeated trials of an evoked or induced brain response. A mixture of different head positions over time adds variance to the data that is not accounted for by the experimental manipulation, thus potentially deteriorating statistical sensitivity. By using a general linear model, head movement related trialbytrial variance can be removed from the data, both at the sensor and source level. This procedure involves 3 step

Preprocess the MEG data, for instance pertaining to an ERF analysis at the sensor level. Note the keeptrials = ‘yes’ when calling ft_timelockanalysis.
% define trials cfg = []; cfg.dataset = 'TacStimRegressConfound.ds'; cfg.trialdef.eventtype = 'UPPT001'; cfg.trialdef.eventvalue = 4; cfg.trialdef.prestim = 0.2; cfg.trialdef.poststim = 0.3; cfg.continuous = 'yes'; cfg = ft_definetrial(cfg); % preprocess the MEG data cfg.channel = {'MEG'}; cfg.demean = 'yes'; cfg.baselinewindow = [0.2 0]; cfg.dftfilter = 'yes'; % notch filter to filter out 50Hz data = ft_preprocessing(cfg); % timelock analysis cfg = []; cfg.keeptrials = 'yes'; timelock = ft_timelockanalysis(cfg, data);

Create trialbytrial estimates of head movement. Here one may assume that the head is a rigid body that can be described by 6 parameters (3 translations and 3 rotations). The circumcenter function (see below) gives us these parameters. By demeaning, we obtain the deviations. In other words; translations and rotations relative to the average head position and orientation.
% define trials cfg = []; cfg.dataset = 'TacStimRegressConfound.ds'; cfg.trialdef.eventtype = 'UPPT001'; cfg.trialdef.eventvalue = 4; cfg.trialdef.prestim = 0.2; cfg.trialdef.poststim = 0.3; cfg.continuous = 'yes'; cfg = ft_definetrial(cfg); % preprocess the headposition data cfg.channel = {'HLC0011','HLC0012','HLC0013', ... 'HLC0021','HLC0022','HLC0023', ... 'HLC0031','HLC0032','HLC0033'}; headpos = ft_preprocessing(cfg); % calculate the mean coil position per trial ntrials = length(headpos.sampleinfo) for t = 1:ntrials coil1(:,t) = [mean(headpos.trial{1,t}(1,:)); mean(headpos.trial{1,t}(2,:)); mean(headpos.trial{1,t}(3,:))]; coil2(:,t) = [mean(headpos.trial{1,t}(4,:)); mean(headpos.trial{1,t}(5,:)); mean(headpos.trial{1,t}(6,:))]; coil3(:,t) = [mean(headpos.trial{1,t}(7,:)); mean(headpos.trial{1,t}(8,:)); mean(headpos.trial{1,t}(9,:))]; end % calculate the headposition and orientation per trial cc = circumcenter(coil1, coil2, coil3) % demean to obtain translations and rotations from the average position and orientation cc_dem = [cc  repmat(mean(cc,2),1,size(cc,2))]';

Fit the headmovement regressors to the data and remove variance that can be explained by these confounds.
% add head movements to the regressorlist. also add the constant (at the end; column 7) confound = [cc_dem ones(size(cc_dem,1),1)]; % regress out headposition confounds cfg = []; cfg.confound = confound; cfg.reject = [1:6]; % keeping the constant (nr 7) regr = ft_regressconfound(cfg, timelock);
Figure
Example statistical results in a singlesubject (baseline vs. task activity contrasts). With ft_regressconfound, sensorlevel statistical sensitivity was increased after tactile stimulation (4050 ms; note the more extreme tscores in the upper panel). In a similar vein, sourcelevel statistical sensitivity was increased after visual stimulation (0500 ms; 65Hz; lower panel).
Figure taken fromStolk et al., NeuroImage 2013.
Practical issues
Some features of this GLMbased compensation method need emphasizing. These points are described in more detail in the ‘Testing the offline GLMbased head movement compensation’ section of Stolk et al., NeuroImage 2013.
First, ft_regressconfound can be applied to timelock, freq, and source data. The estimation of regression coefficients (beta weights of the head position data) is performed separately for each channel and each latency, in the case of timelock data. Consequently, after compensation, the sensor level data cannot be used anymore for source modeling. To employ the GLM based compensation on the source level, single trial estimates for the cortical locations of interest have to be made from the original sensor level data, preferably using a common spatial filter based on all trials. The beta weights are subsequently estimated for each cortical location and the variance in source amplitude over trials that is explained by the head movement is removed. It is therefore recommended to use ft_regressconfound as a final step prior to calling ft_timelockstatistics/ft_freqstatistics/ft_sourcestatistics.
Second, the same trials in the headposition data have to be selected as those present in the MEG data since these two will be fitted. And more or less related; this general linear modeling (GLM) approach only affects the signal variance and not the signal mean over trials (because the constant remains in the data). So when performing a group study, taking the subject mean to the group level statistics will not change these statistics. To benefit from improved statistical sensitivity after using ft_regressconfound, it is advised to take a measure that incorporates the consistency (over trials) of a neural effect to the group level. For instance, the tdescriptive, as obtained using an independent samples ttest on trials of one condition versus that of another. These tvalues can then be tested at the group level for rejecting the nullhypothesis of no difference between conditions (T=0).
Finally, note that the circumcenter function is a helper function that calculates the position (geometrical center of the three localizer coils) and orientation of the head. This saves some degrees of freedom (df=6) as compared to taking into account the x,y,zcoordinates of each coil separately (n=3) as regressors (df=9). If you want to also use the squares, cubes, and derivatives as regressors (to account for nonlinear effects of head motion on the MEG signal), this can save quite a bit of degrees. However, too large a number of covariates can reduce statistical efficiency for procedures. In that case, MATLAB will produce the Warning ‘Rank deficient’. A rule of thumb is to roughly have 10% of the sample size (based on chapter 8 of Tabachnick & Fidell (1996)).
Please cite this paper when you have used the offline head movement compensation in your study:
Stolk A, Todorovic A, Schoffelen JM, Oostenveld R. Online and offline tools for head movement compensation in MEG. Neuroimage. 2013 Mar;68:3948. doi: 10.1016/j.neuroimage.2012.11.047.
Appendix: circumcenter
function [cc] = circumcenter(coil1,coil2,coil3)
% CIRCUMCENTER determines the position and orientation of the circumcenter
% of the three fiducial markers (MEG headposition coils).
%
% Input: X,y,zcoordinates of the 3 coils [3 X N],[3 X N],[3 X N] where N
% is timesamples/trials.
%
% Output: X,y,zcoordinates of the circumcenter [13 X N], and the
% orientations to the x,y,zaxes [46 X N].
%
% A. Stolk, 2012
% number of timesamples/trials
N = size(coil1,2);
%% x, y, and zcoordinates of the circumcenter
% use coordinates relative to point `a' of the triangle
xba = coil2(1,:)  coil1(1,:);
yba = coil2(2,:)  coil1(2,:);
zba = coil2(3,:)  coil1(3,:);
xca = coil3(1,:)  coil1(1,:);
yca = coil3(2,:)  coil1(2,:);
zca = coil3(3,:)  coil1(3,:);
% squares of lengths of the edges incident to `a'
balength = xba .* xba + yba .* yba + zba .* zba;
calength = xca .* xca + yca .* yca + zca .* zca;
% cross product of these edges
xcrossbc = yba .* zca  yca .* zba;
ycrossbc = zba .* xca  zca .* xba;
zcrossbc = xba .* yca  xca .* yba;
% calculate the denominator of the formulae
denominator = 0.5 ./ (xcrossbc .* xcrossbc + ycrossbc .* ycrossbc + zcrossbc .* zcrossbc);
% calculate offset (from `a') of circumcenter
xcirca = ((balength .* yca  calength .* yba) .* zcrossbc  (balength .* zca  calength .* zba) .* ycrossbc) .* denominator;
ycirca = ((balength .* zca  calength .* zba) .* xcrossbc  (balength .* xca  calength .* xba) .* zcrossbc) .* denominator;
zcirca = ((balength .* xca  calength .* xba) .* ycrossbc  (balength .* yca  calength .* yba) .* xcrossbc) .* denominator;
cc(1,:) = xcirca + coil1(1,:);
cc(2,:) = ycirca + coil1(2,:);
cc(3,:) = zcirca + coil1(3,:);
%% orientation of the circumcenter with respect to the x, y, and zaxis
% coordinates
v = [cc(1,:)', cc(2,:)', cc(3,:)'];
vx = [zeros(1,N)', cc(2,:)', cc(3,:)']; % on the xaxis
vy = [cc(1,:)', zeros(1,N)', cc(3,:)']; % on the yaxis
vz = [cc(1,:)', cc(2,:)', zeros(1,N)']; % on the zaxis
for j = 1:N
% find the angles of two vectors opposing the axes
thetax(j) = acos(dot(v(j,:),vx(j,:))/(norm(v(j,:))*norm(vx(j,:))));
thetay(j) = acos(dot(v(j,:),vy(j,:))/(norm(v(j,:))*norm(vy(j,:))));
thetaz(j) = acos(dot(v(j,:),vz(j,:))/(norm(v(j,:))*norm(vz(j,:))));
% convert to degrees
cc(4,j) = (thetax(j) * (180/pi));
cc(5,j) = (thetay(j) * (180/pi));
cc(6,j) = (thetaz(j) * (180/pi));
end