SPM for programmers

When you have already preassumptions with regard to the coupling strength of the neural parameters or the hemodynamic parameters you can start the estimation of a DCM with a certain starting position.

All you have to do is, create a variable DCM.M.P and assign new values to it, e.g:

load dcm_file
DCM.M.P = DCM.M.pE;
% alter DCM.M.P
save dcm_file DCM

The EM-algorithm will begin with the parameters specified in DCM.M.P.

DCM is implemented in a framework that can be used to develop and test other models using the Expectation Maximization algorithm and Bayesian techniques. To understand how to use other models within this framework, I can recommend to look at spm_nlsi.m. After the return-command, at the end of the script, there is a nice short example, how to setup everything in order to test your own model.

And by the way: the code for the dynamic causal model can be found in spm_fx_dcm.m (the neural signal is in y(:,1), some hemodynamic parameters are stored in y(:,2:5)) and spm_gx_dcm.m (the hemodynamic forward model).

An interesting read about Psychophysiological Interactions and some issues related to deconvolution and orthogonality of the regressors can be found here at the FMRIB Centre.

Gerrit from the blog Sustainable Research posted a nice article about the screen-command for the linux console that allows you to run scripts in the background without an open console window. This can be also useful, when you want to run SPM-scripts without a user interface and remote on a server. You just start the script within a screen, exit the screen and close the remote connection. And with Putty for Android or iPhone you can supervise your scripts from everywhere :).

I suppose that this is not the most efficient way to do this, but it works. With the following commands you can easily write a script that displays any kind of nifti-image (raw files, mask, beta-maps) across all subjects and saves it as a JPG image. This can be handy for quality assurance:

spm_image('init','niftifile.nii') % display a nifti-file
spm_orthviews('Reposition', [22 0 -21]) % jump to a certain coordinate
fs = get(0,’Children’);
res = getframe(fs(1));
imwrite(res.cdata, ‘niftifile.jpg’); % save the window as JPG

This should work for SPM5 and SPM8. Any suggestions how to improve the “save as jpg” procedure?

I guess I am not the only one who uses matlab scripts to analyze large amounts of fMRI data. To get rid of the popups of SPM during the processing, you can start Matlab either with

matlab -nodesktop


to solely start matlab on the command line, or you start it with

matlab -noFigureWindows


which starts Matlab as normal but prevents the depiction of any additional windows.

This post was inspired by a question that I received by email. Thanks!

The high-pass filter can have a great effect on the subsequent parameter estimation. To explore, how a certain setting (default is 128) influences the time course of the regressor you can use the spm_filter function. First, we build a HRF-regressor:

reg = [zeros(10,1); ones(10,1); zeros(30,1); ones(10,1); zeros(30,1); ones(10,1); zeros(30,1)];


bf = spm_get_bf %TR is 2.2 in our case
U.u = reg;
U.name = {'reg'};
% create the regressor convolved with the HRF
convreg = spm_Volterra(U, bf.bf);


In your SPM.mat you can find this kind of regressors under SPM.xX.X.

spm_filter is a function that either creates the low frequencies to be removed or removes those frequencies from the time course. The input variable K defines all necessary parameters:

K.RT = 2.2;
K.row = 1:130; % 130 corresponds to the length of convreg
K.HParam = 128; % cut-off period in seconds


Here, I defined K.row as it also happens internally. But I guess, ones(1,130) should also work, as it seems that K.row just serves as a mask. Now, we compute the low frequencies that should be removed from convreg.

nK = spm_filter(K);


nK.X0 contains the frequencies that can be removed. A second call of spm_filter removes these frequencies.

Y = spm_filter(nK, convreg);


In SPM.mat these regressors are stored in SPM.xX.nKX and are depicted in the design matrix. As you can see, the high-pass filter cut-off period can greatly affect the regressors that are used for the parameter estimation:

K.HParam = 256;
figure; plot(spm_filter(spm_filter(K), convreg));


K.HParam = 64;
figure; plot(spm_filter(spm_filter(K), convreg));


If you experiment with DCM, it might be interesting to know that you can build synthetic input data for DCM with the spm_dcm_create-function. Just type spm_dcm_create on the console, press enter and follow the steps. Other interesting units for DCM are: spm_dcm_estimate (the estimation procedure) and spm_dcm_compare (SPM8) for comparing two or more models.

With the help of two SPM-functions you can manually convolve time series in Matlab. This might help to understand how the convolution with the HRF (or any other function) affects the time series. Here is an example:

% create a time series
reg = [zeros(10,1); ones(10,1); zeros(10,1); ones(10,1); zeros(10,1)];
plot(reg)


% create a basis function, in this case the hrf
bf = spm_get_bf
% with a tr of 2.2 and HRF as basis function you receive
dt: 2.2000
name: 'hrf'
length: 33
order: 1
bf: [15x1 double]
% these are necessary definitions for the spm_Volterra-function
U.u = reg;
U.name = {’reg’}; % U needs a name, but the string has no meaning
% convolve reg with the hrf
convreg = spm_Volterra(U, bf.bf);
plot(convreg)


To convert matrix coordinates for a given Nifti/Analyze-image into the image coordinates and vice versa, you need to load the transformation matrix of the image into the workspace:

vol = spm_vol('image.img')
vol.mat % displays the transformation matrix


% trasform from matrix/voxel coordinates to image coordinates (the fourth element has to be 1)
vol.mat * [x y z 1]‘
% transform image coordinates (e.g. mni space) to matrix coordinates
inv(vol.mat) * [x y z 1]‘