%%%% Tutorial on the basics of convolving with an HRF, using Matlab %%%% Written by Rajeev Raizada, July 23, 2002. %%%% %%%% There is also a follow-up file about the structure of %%%% an fMRI-analysis design matrix: design_matrix_tutorial.m %%%% %%%% Neither file assumes any prior knowledge of linear algebra or Matlab. %%%% %%%% Probably the best way to look at this program is to read through it %%%% line by line, and paste each line into the Matlab command window %%%% in turn. That way, you can see what effect each individual command has. %%%% %%%% Alternatively, you can run the program directly by typing %%%% %%%% hrf_tutorial %%%% %%%% into your Matlab command window. %%%% Do not type ".m" at the end %%%% If you run the program all at once, all the Figure windows %%%% will get made at once and will be sitting on top of each other. %%%% You can move them around to see the ones that are hidden beneath. %%%% %%%% Anything in this program with a % sign in front of it is a comment. %%%% These will probably show up red in your Matlab Editor. %%%% Everything else is a Matlab command, that you can copy and %%%% paste into the Matlab command window. %%%% %%%% Please mail any comments or suggestions to: raj@nmr.mgh.harvard.edu % 1. The brain produces a fairly fixed, stereotyped blood flow response % every time a stimulus hits it. That's the HRF (haemodynamic response % function). % % 2. Once one of these resonses starts, it just does its thing until it % has finished. i.e. it peaks, it goes back down to zero, it undershoots % a bit, then it settles back to baseline. All this was kicked off as a % a result of the stimulus coming in, and it will run to completion % regardless of what stimuli, if any, come in later. % (This is just an approximation of what really happens in the brain, of % course, but it turns out that it isn't all that far from being true). % % 3. What we just described is exactly what the process called "convolution" % is. To convolve, you need two things. First, you need an "impulse % response function" (IRF), i.e. the chain of events that will be started % off every time an "impulse" happens, e.g. the brain sees a stimulus. % And in the brain, the response to each impulse is the haemodynamic % response function --- the HRF. Every time a stimulus comes in, one % of these HRFs gets started off, and then it runs it course (peak, % undershoot, baseline) over the next 16 seconds or so. % % Second, you need a vector to be convolved with that impulse response % function. A vector is just a bunch of numbers in a row. In this case, % our vector to be convolved is a list of whether a stimulus is shown % or not, at each moment in time. So, we have a bunch of numbers % lined up in a row, with time being represented by our position % along the row. % This type of vector is often referred to as a time-series. % % In this case, the time-series is describing whether a stimulus is % being shown at that moment in time or not. % If there's no stimulus, we put a zero. % If there is a stimulus, we put a 1. % The first number is what is happening at time t=0, % the second number is what is happening at t=1 etc. % % Here's a key point: every stimulus will kick off its own HRF, and so % it can often happen that the previous stimulus's HRF hasn't finished % by the time the next stimulus, with its new HRF comes in. % So, the two HRFs will overlap in time. % % In convolution, if you get different impulse response functions % overlapping in time, you can work out what the total response at % that moment in time will be simply by adding up the individual % responses. % % To say that you can get the total response just by adding up the % individual overlapping responses is exactly what it means to say % that the system is LINEAR. If something is linear, then its % total response to separate inputs is just the sum of what its % individual responses would have been to the individual inputs. % % Does the brain do this? Does it just "add up" blood flow responses % which are overlapping? % % Yes! (to a reasonable approximation). % And that is why convolution does a reasonable job of describing % what the brain does. % % To see more detail about the math of convolution, % look at the companion program math_of_convolution.m % % 4. In fMRI, the data that we get for each voxel is the blood % oxygenation at that voxel at each point in time. % We also know what the presentation times of our stimuli were. % So, what any fMRI-analysis program does is this: % It takes the time series of stimulus onsets, and convolves it % with the HRF. This gives a *prediction* of the blood flow response % that we should get a given voxel, if that voxel is responding to the % stimuli. % % Then, we take that prediction, and go and compare it against the measured % data. And we see how well they match. %%%%%% First let's create a made-up vector that looks like an HRF %%%%%% A vector is just a bunch of numbers in a row. %%%%%% In Matlab, we make a vector by putting the row of numbers %%%%%% inside square brackets [ ] %%% Give one value for each 1-second time point. %%% So, the 18 numbers in the vector below correspond %%% to the values of the HRF from time=0 to time=18 %%% %%% Note that the important thing about the HRF here is %%% just the overall *shape* that it has, not the exact values %%% of the numbers. So, it doesn't mean anything that the max %%% value in the numbers below is 9, or 9.2 or whatever. %%% It's the shape of the HRF curve over time that matters. %%% %%% [ Actually, this isn't completely true. %%% For some math-related purposes it turns out that it's convenient %%% to have the HRF sum to 1, and this is what SPM does. %%% But it would be a distraction to worry about that right now. %%% Let's just define an HRF with the right overall shape, and plot it. ] hrf = [ 0 0 1 5 8 9.2 9 7 4 2 0 -1 -1 -0.8 -0.7 -0.5 -0.3 -0.1 0 ] %%% When this command runs, Matlab will show the value of the %%% new variable "hrf" in the Command Window. %%% That's because there is no semi-colon at the end of the line. %%% If it did have a semi-colon at the end, Matlab wouldn't show it %%%%%%%%%%%%% Let's plot the HRF %%%% In Matlab, the command to do this is called "plot" %%%% You tell it three things: %%%% 1. The x-coordinates of the points to plot, lined up in a vector %%%% 2. The y-coordinates of the points to plot, lined up in a vector %%%% 3. (Optional) What styles of line and markers to use for the plot. figure(1); % Make a new figure window, Fig.1 clf; % Clear the figure %%%% Make a vector of time-values to use as the x-coordinates time_vector = [ 0:18 ]; %%% 0:18 means "A row of numbers, starting at zero, %%% and going up to 18". %%% The numbers increase in steps of 1, which is the default. %%% The square brackets make it a vector. %%% The semi-colon at the end means that this %%% new variable will NOT display in the command window. %%% If you want to see that value of it, type %%% time_vector %%% into the command window, without a semi-colon after it. %%%% Plot HRF against time, with time_vector as the x-coordinates, %%%% hrf as the y-coordinates, and with the line style being a %%%% solid line with circles on it. %%%% The way to specify "solid line with circles on it" is to %%%% make the third argument be 'o-' plot(time_vector,hrf,'o-'); grid on; % Overlay a dotted grid on the graph %%%% Label our axes and give the graph a title xlabel('Time (seconds)'); ylabel('fMRI signal'); title('The typical shape of an HRF'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Now let's make a vector of that is the time-series of our % stimulus onsets, so that we can convolve that vector with the HRF, % as described in the intro above. % % A vector is just a bunch of numbers in a row. In this case, % they are lined up in time. i.e. the first number is what is happening % at time t=0, the second number is what is happening at t=1 etc. % Our vector to be convolved is a list of whether a stimulus is shown % or not, at each moment in time. If there's no stimulus, we put a zero. % If there is a stimulus, we put a 1. % % Suppose that our scan lasts 60 seconds, % and that we use a visual stimulus: we flash a light at time t=10. % % Let's make an stimulus-time-series vector that has one element % for each second of time. % Then the vector for that light stimulus will % have a 1 in the 10th position, meaning t=10, % and will be zeros everywhere else. % % So, it's 9 zeros, then a 1, then 50 more zeros, % making the vector have 60 entries altogether, % for the 60 second scan. % % Note that by just putting a single 1 to represent the flash % of light, we are saying that the flash happened suddenly, % just at that moment of time, i.e. it was an event, as opposed % to a block (also called an "epoch"), in which the stimuli % are spread out over a longer period of time. E.g. a stimulus-block % might last twenty seconds, rather than an event, which happens % at a single moment. % % This sudden flash of light will spark off a sudden neural event, % in this case a sudden burst of neural firing in visual cortex. % % In Matlab, we can make nine zeros in a row, or 50 zeros in a row, % using the command "zeros". % The first argument is how many rows we want (just one row in this case), % and the second argument is how many columns we want ( nine or fifty ). % % To make the whole big vector with 9 zeros, a single 1, and then 50 zeros, % we can just put the separate vectors next to each other in a row, % and then put the whole thing inside square brackets, like this: first_light_stim_time_series = [ zeros(1,9) 1 zeros(1,50) ]; %%%%%%% Translation-guide % % Here we have made a time-series whether for every time-point (t=0,t=1,etc.) % we are explicitly saying whether or not a stimulus is being shown. % In this case, the light is shown at t=10, so put a 1 as the 10th element, % and for the rest of the time no lights are shown, % so we put 0 everywhere else. % % However, when using an fMRI-analysis package, people typically don't % explicity make the whole time-series. They simply tell the package % what the onset-times of the stimuli are, and let the package build % the time-series of 0's and 1's on its own. % % So, for the case where we have a single light stimulus at t=10, % the info entered into the analysis package would be something like: % % onset_times = [ 10 ]; % % Similarly, if we flashed one light at t=10 and another at t=30, % then we'd enter: % onset_times = [ 10 30 ]; % % In this tutorial, and also in design_matrix_tutorial, % we're going to explicitly make the stimulus-time-series vectors, % rather than only supplying the onset-times, because these time-series % vectors are the things that actually get convolved with the HRF. % In an fMRI-analysis package, these time-series vectors would % get made and convolved "behind the scenes", but we want to see % everything in full view! %%%% Now that we have the stimulus-time-series vector and the HRF vector, %%%% we can convolve them. %%%% In Matlab, the command to convolve two vectors is "conv", like this: hrf_convolved_with_stim_time_series = conv(first_light_stim_time_series,hrf); %%%%% Let's plot the result figure(2); % Make a new figure window, Fig.2 clf; % Clear the figure subplot(2,1,1); % This is just to make the plots line up prettily % The first number is how many rows of subplots we have: 2 % The second number is how many columns: 1 % The third number is which subplot to draw in: the first one % So, we end up with two plots stacked on top of each other, % and we draw in the first one (which is the upper subplot) stem(first_light_stim_time_series); % Show when the stimulus is presented. % We're using the command "stem" to plot here, instead % of the more standard command "plot". % "Stem" makes a nice looking plot with lines and circles. % This type of plot is good for showing discrete events, % such as stimulus onsets. grid on; % Overlay a dotted-line grid on top of the plot legend('Time-series of light stim'); % Make a box to say what the plot-lines are showing axis([0 60 0 1.2]); % This just sets the display graph axis size % The first two numbers are the x-axis range: 0 to 60 % The last two numbers are the y-axis range: 0 to 1.2 xlabel('Time (seconds)'); ylabel('Stimulus present / absent'); subplot(2,1,2); % Draw in the second subplot (the lower one) plot(hrf_convolved_with_stim_time_series,'rx-'); % The argument 'rx-' means: Draw in red (r), % use cross-shaped markers (x), and join them with a solid line (-) grid on; legend('Stimulus time-series convolved with HRF'); axis([0 60 -2 15]); xlabel('Time (seconds)'); ylabel('fMRI signal'); % You might be wondering why we didn't specify an x-coordinate % vector in the two plot commands just now. % % e.g. plot(hrf_convolved_with_stim_time_series,'rx-') % % We gave the y-coord values: the vector "hrf_convolved_with_stim_time_series" % But we didn't give any x-coord values. % When the plot command is given just one vector, it automatically % plots the first value at x=1, the second at x=2 etc. % Because the time-axis that we want in these plots % starts at 1 and goes up in steps of 1, this default is fine for us here. %%%%%% Display another light, now at time 30 %%%%%% Make a new stimulus-time-series vector for it, %%%%%% just like we did above. %%%%%% Except now there are 29 zeros, then a 1, then 30 zeros second_light_stim_time_series = [ zeros(1,29) 1 zeros(1,30) ]; %%%%%% Because the two vectors first_light_stim_time_series %%%%%% and second_light_stim_time_series are the same length %%%%%% (they are both 60 elements long, for the 60 second scan), %%%%%% we can add them together, element-by-element. %%%%%% %%%%%% This will make a new vector that has a 1 as the 10th element, %%%%%% and a 1 as the 30th element, and is zeros everywhere else. %%%%%% %%%%%% This represents us showing lights at t=10 and t=30 all_lights_time_series = ... % Three dots mean "continued on next line" first_light_stim_time_series + second_light_stim_time_series; %%%%% Do the convolution using the Matlab command "conv" hrf_convolved_with_all_lights_time_series = conv(all_lights_time_series,hrf); %%%%% Let's plot the result of this convolution, in a new figure figure(3); clf; % Clear the figure subplot(2,1,1); % This is just to make the plots line up prettily stem(all_lights_time_series); grid on; legend('Time-series of light stim'); axis([0 60 0 1.2]); xlabel('Time (seconds)'); ylabel('Stimulus present / absent'); subplot(2,1,2); plot(hrf_convolved_with_all_lights_time_series,'rx-'); grid on; legend('Stimulus time-series convolved with HRF'); axis([0 60 -2 15]); xlabel('Time (seconds)'); ylabel('fMRI signal'); %%%%%% Now let's try moving the two stimulus onsets closer together %%%%%% This will show how the two HRFs add together to make the measured signal second_light_stim_time_series = [ zeros(1,15) 1 zeros(1,44) ]; % A flash of light, with onset-time t=16 all_lights_time_series = ... % Three dots mean "continued on next line" first_light_stim_time_series + second_light_stim_time_series; % Onsets at t=10 and t=16 hrf_convolved_with_all_lights_time_series = conv(all_lights_time_series,hrf); hrf_from_first_light = conv(first_light_stim_time_series,hrf); hrf_from_second_light = conv(second_light_stim_time_series,hrf); figure(4); clf; % Clear the figure subplot(3,1,1); % This lines up three subplots on this figure window stem(all_lights_time_series,'bo-'); % Blue circles, solid line grid on; legend('Time-series of light stim'); axis([0 60 0 1.2]); ylabel('Stimulus present / absent'); subplot(3,1,2); plot(hrf_convolved_with_all_lights_time_series,'rx-'); % Red crosses, solid line grid on; legend('Stimulus time-series convolved with HRF'); axis([0 60 -2 15]); ylabel('fMRI signal'); subplot(3,1,3); hold on; % "Hold" is one way of putting more than one plot on a figure plot(hrf_from_first_light,'b-'); % Blue solid line plot(hrf_from_second_light,'m--'); % Magenta dashed line hold off; grid on; legend('HRF from first flash of light','HRF from second flash of light'); axis([0 60 -2 15]); xlabel('Time (seconds)'); ylabel('fMRI signal'); % Now let's put in a third onset, so we have onsets at t=10, t=13 and t=16. % Having all these trials following each other in a row is starting to % look like a blocked design. % When we plot the HRF that these closely-spaced trials evoke, in Figure 5, % note how the individual HRFs from each trial all start to lump together % into one big HRF. That's what the HRFs look like in a blocked design. third_light_stim_time_series = [ zeros(1,12) 1 zeros(1,47) ]; % A flash of light, with onset-time t=13 all_lights_time_series = ... % Three dots mean "continued on next line" first_light_stim_time_series + second_light_stim_time_series + ... third_light_stim_time_series; % Onsets at t=10, t=13 and t=16. hrf_convolved_with_all_lights_time_series = conv(all_lights_time_series,hrf); hrf_from_first_light = conv(first_light_stim_time_series,hrf); hrf_from_second_light = conv(second_light_stim_time_series,hrf); hrf_from_third_light = conv(third_light_stim_time_series,hrf); %%%%%%%%%%%%% Let's plot the results of convolving the time-series of %%%%%%%%%%%%% these three closely spaced flashes of light with the HRF figure(5); clf; subplot(3,1,1); % This lines up three subplots on this one stem(all_lights_time_series,'bo-'); grid on; legend('Time-series of light stim'); axis([0 60 0 1.2]); ylabel('Stimulus present / absent'); subplot(3,1,2); plot(hrf_convolved_with_all_lights_time_series,'rx-'); % Red crosses, solid line grid on; legend('Stimulus time-series convolved with HRF'); axis([0 60 -2 16]); ylabel('fMRI signal'); subplot(3,1,3); hold on; % This is one way of putting more than one plot on a figure plot(hrf_from_first_light,'b-'); % Blue solid line plot(hrf_from_second_light,'m--'); % Magenta dashed line plot(hrf_from_third_light,'g-'); % Green solid line hold off; grid on; legend('HRF from first flash of light','HRF from second flash of light', ... 'HRF from third flash of light'); axis([0 60 -2 15]); xlabel('Time (seconds)'); ylabel('fMRI signal');