Data Simulation

Convolving a Simulated fMRI Time Series

Scenario: A multilayer feedforward network was created in MikeNet with manually-set weights to produce an artificial system of interconnected regions that would behave in a predictable manner according to inputs of different types (e.g., input of type A would activate some regions, deactivate other regions, and have no influence on a third set of regions; the simultaneous input of A and B was required to maximally activate one of the regions, etc.). A continuous time series of simulated activations was read from each of the network regions as inputs A and B were given, where each region's activation could range from 0.0 to 1.0, and was dependent on the activity of the region that connected into it at the previous time point. This introduced the concept of a lag into the network, such that regions directly receiving input were activated faster than those that were further down the connectivity train.

Problem: to make the synthetic activations as realistic as possible, I wished to convolve the simulated time series with SPM's canonical HRF, under the assumption of additivity of activity. It seemed that this might be possible using the spm_Volterra command, but alas, no. In the end, I accomplished this in two steps:

Step 1: Compute the HRF basis function

With SPM in your MATLAB path, we can call the spm_get_bf function:

%first, create a data structure with some defining parameters
xBF.name='hrf';
xBF.length=16;
xBF.dt=0.25;
xBF.order=2;
xBF.T=16;
xBF.dt=2;
xBF=spm_get_bf(xBF); %if you plot(xBF.bf) you will see the characteristic canonical HRF

Having the HRF helps, but that's the immediate response to one event. I wish to convolve this function with an entire simulated time series, similar to how the spm_Volterra function convolves this function with a time series generated from a vector of event onsets.

Step 2: Multiply the HRF basis function by each time point

Now at this point, there's an HRF, and there's one or more time series in a matrix, M, which could all be 0,1 or be values in between, representing the intensity of the event at time t. Matrix multiplication doesn't get the job done. Instead, it seems the way to do it is to make a big matrix, 1 row per time point, and then multiply the value at that time point by the HRF, and store it in the next open row. At the end of it all, you will have a time points × (time points + HRF window) matrix. Think of each row as being a moving window. Summing these windows, we collapse our time points together to get what seems to be the convolution of the HRF with the time series.

cols=[4:11]; %supposing column i of the time series matrix M contains the time series you wish to convolve
bf=xBF.bf;
winsize=length(bf); %the window seems to be xBF.length + 1
simTS=zeros(size(M,1)+winsize, size(M,1));
allsimTS=[]; %keep appending rows to the empty matrix
for i=cols
   for t=1:size(M,1)
       thisscale=M(t,i);
       tbf=(bf*thisscale)';
       simTS(t,t:t+winsize-1)=tbf;
   end
   simTS=sum(simTS); %summing all the temporally overlapping time series windows to collapse into a single time series
   allsimTS=[allsimTS;simTS];
end

Did it work? Kind of (?)