Today we will begin a several part lab exercise series devoted to the analysis of functional MRI data. In our next session, we will use the FSL (FMRIB Software Library) package created by the Oxford group to analyze our data in the manner typical of the field. However, FSL hides most of the analysis steps from direct observation, and describes those steps in complex statistical terms. So today I want you to analyze your data 'by hand' using Matlab commands that are (relatively) easy to understand and to follow. I want you to interact closely with your data so that you have a good understanding of the signal and noise characteristics of fMRI data. All analysis, no matter how sophisticated, starts with the raw signal. That is where we begin today!

• Preprocessing fMRI data,

• Explore different fMRI data sets to observe how a simple task alters voxel intensity in a fMRI time series.
  • Visually inspect a data set to identify activated voxels.
  • Use a simple correlational approach to identify activated voxels using one or two templates.

• n/a

We will use different datasets throughout the lab. They will be described in detail in the relevant parts of the lab. What they all have in common is that they are what we would call Block designs.

A “block design” is an an fMRI experiment in which a stimulus (or task) is presented for several seconds, often followed by several seconds of rest. Let's say we were running a block design experiment in which we wanted to investigate color sensitivity in the visual system. We would present blocks of color photographs and blocks of greyscale photographs. Importantly, within a given block the participant would probably be shown several individual stimuli, but they would all be of the same category. So within a 10 second greyscale block the participant might see 10 different greyscale pictures, each displayed for 1 sec.

We will start here by looking at a visual evoked response, followed by a motor task. Our goal for each task will be to find a voxel in the brain that is activated by a given task.

Data for Part 1 of the lab can be found in /Users/hnl/Desktop/input/fmri/high_res/. This is high-resolution functional data (2mm x 2mm x 2mm) that will allow us to roughly identify relevant anatomy. Data were collected for two tasks:

• visual response task
• motor task

Each of the tasks used a block design (see above) that had the following common construction:

• 115 volumes with TR = 2 sec
• Each run of each task had the following construction
  • Task A ⇒ Task B ⇒ Task A …
  • Each Task block lasted 12 seconds.
  • There were 7 blocks for each task for a total of 14 blocks

• Task A checkerboard images presented to the LEFT hemifield.
• Task B checkerboard images presented to the RIGHT hemifield.

• Task A participants were asked to repeatedly squeeze their LEFT hand.
• Task B participants were asked to repeatedly sqeeuze their RIGHT hand.

Note: These tasks are essential the same as those used in the three original 1992 fMRI papers.

1. Open Terminal and navigate to the date directory

 cd ~/Desktop/input/fmri/high_res

2. Open FSLeyes

fsleyes

3. Load the data file:

• File → Add from file
• select tb9611_checkers_run01.nii.gz

By default, FSLeyes displays the very first of the full-brain acquisitions (aka volumes). In the Location box at the bottom of the screen you will see that the Volume is 0.

4. To appreciate that we now have a time dimension, let's watch a movie of the data.

• Click on the settings button in the upper left corner (see image below)
• Slide the Movie update rate slider so that it's about three quarters of the way to the right
  • This controls the speed at which the movie will play
  • Close the View settings window
• Click on the Movie mode button (see image below) and the movie will being to play
  • To have a better view, you might want to turn off the crosshairs by clicking on the Crosshair button (see image below)

5. Click on a voxel somewhere in the medial occipital region.

Knowing what you know about the task parameters, can you see the task-related variation in the intensity changes in the movie? Look carefully in visual cortex. Do you see the activation signal?

You might notice some pulsation going on but probably find very it difficult to pinpoint regions with task-related activity. This should give us a hint that the signal-to-noise ratio is not too favorable or given the massive amounts of different voxels - we cannot easily detect a signal with the naked eye.

6. Stop movie mode by clicking on the Movie mode button.

A moment ago you viewed a movie of the signal in each voxel (brightness) changing over time. Here, we will instead look at the time-series from specific voxels. A time-series for this study would simply be a list of 115 values. At each voxel, we'd have a signal intensity for each of the 115 volumes. We actually have 624,000 separate time-series; one for each voxel (104 x 100 x 60). But fret not, we'll only be looking at one at a time.

If we found a responsive voxel-one that reflected a response to the task-what do you think its time-series would look like? Remember that each hemifield is stimulated in 12-second blocks. What part of the brain do you expect to respond to this stimulation? What do you expect that response to look like?

7. Let's look at the time-series data from individual voxels

• View → Time series
  • You'll be doing this a lot tonight, so you might want to remember the keyboard shortcut to open a timeseries: command + 3

You should now have a new window at the bottom of the screen that displays the voxel time-series. The y-axis represents the intensity of the BOLD signal and the x-axis represents the different time-points.

8. Now click around in the visual cortex until you find a really nice looking time-series. You might consider focusing your search in and around the calcarine sulcus, which is where V1 and V2 are located. Find a voxel with a highly responsive time-series in the right hemisphere. In other words, you should see the waveform go up and down with the timing of the task design. Check with me or a classmate to be sure you've found a good one.

In this task, participants were asked to squeeze their hands at particular times. Thus, we should find particular regions in the brain (e.g., motor regions) that will show increases in brain activity when the participant is squeezing their hand.

9. Load the data file:

• File → Add from file
• Select tb9611_handsqueeze_run01.nii.gz
• Make the checkerboard data invisible by clicking on the blue eye in the Overlay list. Alternatively, you can highlight it and click the minus button to remove it altogether.

You can focus your search for a good-looking time-series to the motor cortex.

I'm afraid that finding the hand area of motor cortex will be a bit harder than finding the early visual areas in the previous part of the lab. Below are some MRI images indicating the location of this region (click to embiggify). Note: most of these images are 1 mm³ T1-weighted images, whereas you're data are 8 mm³ T2*-weighted images.

Another hint to finding the hand area is to first find the superior frontal sulcus, which runs anterior-posterior, indicated by the arrow heads in the image below. The central sulcus with the omega shape (knob) should then be behind it. The star in this image shows the hand area in the left hemisphere. The omega shape is clearly seen in the mirror location of the right hemisphere.

Even after finding the correct brain region, you might have some trouble finding a good voxel. At this point, I would suggest you use the keyboard arrow keys to move the crosshair by one voxel in a systematic way.

Here's a slice that includes the hand-area…can you spot the motor cortex hand region?

Here you will try to find task-activated voxels in the functional MRI data assigned to you.

Over the course of the next several exercises, you will work with the data from the same localizer experiment.

• The data are located in this directory:

/Users/hnl/Desktop/input/fmri/loc/data/nifti

• There are 17 subdirectories within this nifti subdirectory, each containing the data for one subject

2545 2552 2553 2554 2585 2731 2766 2767 2814 3850 3851 3855 3866 5743 5744 5769 5770

• Within each of these 17 subdirectories are data for as many as four localizer tasks. We will only be using the first two: a face localizer (face) and a motor localizer (motor)
  • the face task: viewing pictures of scenes (houses) vs. pictures of faces
  • the motor task: left vs. right hand movement

• For each of these localizer tasks and subjects, there were 2 or 3 runs. For example, subject 2545 has three face runs, each saved to a separate file:
  • 2545_face_run1.nii.gz
  • 2545_face_run2.nii.gz
  • 2545_face_run3.nii.gz
    • note, All subjects only have two motor runs. Most subjects have 3 runs of the remaining conditions.

Each of the localizer tasks used a block design (see info box above) that had the following common construction:

• 150 volumes with TR = 2 sec
  • 153 volumes were initially collected, but the first 3 volumes were deleted to allow the spins to reach a steady-state magnetization.
• Each run of each task had the following construction
  • Task A ⇒ Rest ⇒ Task B ⇒ Rest ⇒ Task A …
  • Each Task and Rest block lasted 12 seconds.
  • In each run, there were
    • 6 Task A blocks
    • 6 Task B blocks
    • 12 intervening Rest blocks.
  • Task A started at the 6-second mark of each run with the first brain volume occurring at time zero.

• Subjects were shown a visual display consisting of one of the following three symbols
  • ««
  • ===
  • »»
• «« was Task A and indicated that the subject should make rapid alternating button press responses with the index and middle finger of their left hand
• === indicated that the subject should rest and make no button presses
• »» was Task B and indicated that the subject should make the same alternating press responses with the fingers of their right hand

• Task A consisted of a series of eight pictures of scenes.
• Task B consisted of a series of eight pictures of faces.
• Rest consisted of a fixation cross on a gray background.
• For both Face and Scene blocks, subjects covertly counted the number of immediate picture repetitions - a so-called “one-back” task.
  • Subjects reported their total count of repetitions at the conclusion of the run.
  • The purpose of this task was simply to ensure that participants continued to pay attention to the pictures.

1. Use the steps above (see here for a refresher) to load your subject into FSLeyes

• You might want to close FSLeyes and then reopen it so you have a clean slate.
• You can look at run1, run2, or run3 (if it exists) for these first analyses.

2. Knowing what you know about the task and its temporal structure (i.e., when blocks come on and off) - try to locate individual voxels whose time courses vary with the task timing. Think about where a motor or face task should activate the brain.

I can assure you that there are voxels that clearly show task variation, and you may be amazed when you find them. However, don't spend too much time looking because we are also going to use a simple statistical procedures in the next part to find these voxels

3. If you find a voxel or voxels that appear to be activated by the task, call me over and show it to me.

You should now have a pretty good idea that it is not easy to locate a functional MRI activation by simply eye-balling the data. The change in the raw MR intensity for active voxels (the 'signal') is not much greater than the change in raw MR intensity for unactivated voxels (the 'noise'). That is, there is poor signal-to-noise for fMRI activation.

We will use a MATLAB script, fmri_lab_script2_2025.m, to help us find voxels at which the signal correlates with our task. You should have already copied the scripts for today's lab to your scripts directory.

1. Edit this script using the following command in MATLAB.

  edit '/Users/hnl/Desktop/scripts/fmri_lab_script2_2025.m'

What do we expect as our expected signal? Well, in the absence of any better ideas - why don't we enter a waveform that looks like the task timing.

Inside your script is a variable called template. It is a 1-D matrix (a vector) that contains 150 zeros - one entry for each volume of data within the run. Remember, one volume was acquired every 2s and we have a total of 150 of these volumes in each run.

2. Create an expected waveform for your task using 1's and 0's. Put a 0 when you expect there to be no activation in a volume, and put a 1 when you expect there to be activation in a volume. Note that template takes up a few lines of the script. The , … at the end of each line means that the line is continued to the next line. Thus, template is a single vector with 150 elements.

You may also note that I put the lines together in such a way to facilitate your task. The first and last line of 0s represent rest while there are six main lines in between that represent the six repetitions of the Task A - Rest - Task B - Rest structure.

To complete the template below, think about

1. the duration each block of Task A, each block of Task B, and each block of Rest
2. during which of those block periods we might expect brain activity?
   1. For now, let's just discriminate between Task and Rest. That is, don't worry about treating Task A and Task B differently; we'll do that later in the lab.

The initial template should look like this. Fill in the 1s where appropriate.

template = [0 0 0,... % Rest
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,... % Task A - Rest - Task B - Rest
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,... % Task A - Rest - Task B - Rest
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,... % Task A - Rest - Task B - Rest
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,... % Task A - Rest - Task B - Rest
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,... % Task A - Rest - Task B - Rest
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,... % Task A - Rest - Task B
    0 0 0];   % Rest

3. Run the script by calling the function from the MATLAB command line

fmri_lab_script2_2025('SUBJ','TASK',1)

fmri_lab_script2_2025('9999','motor',1)

The script will produce useful output:

• In your MATLAB command window you will see the Number of voxels exceeding minimum correlation of 0.30 ⇒ .
  • This is a very rough indication of how many voxels were 'activated' (i.e., identified by correlation with your template).
• FSLeyes will plot the functional data with the correlation results overlaid.
  • You'll probably want to press option + r to recenter the brain.
  • The contrast limits are set to .3 to .8 with a red-yellow color map. This means that voxels in which the expected time-series and the actual time-series correlated with a coefficient (r) of at least .3 are colored in. The more yellow the voxel, the larger the correlation coefficient.
  • By default, the cursor will be placed on the voxel that had the strongest correlation with your template.

Your display should look something like this:

Let's take a closer look at the relationship between the expected (i.e., your template) and the actual time-series.

4. Plot the time-series

• Make sure to highlight the functional MRI data in the Overlay list (this is the file that does not have corr in the name).
• Press command + 3
• Change the Plotting Mode to Normalised
  • See the red arrow two figures down

Now you should see the time-series from the voxel plotted at the bottom of the screen. The y-axis represents the intensity of the BOLD signal and the x-axis represents different time-points:

You can probably see that the signal goes up and down in time with the task. But we can compare this more directly by overlaying your template.

5. Overlay your model.

• Press the button to Import data series from a text file
  • See the green arrow in the figure above
• Select the text file you just created with your MATLAB script.
  • It will be in your data directory and named SUBJ_TASK_model1.txt
• Click Ok in the Scaling factor pop-up window

You should now see your model (the template you created) and the actual time-series, like this:

Your results probably look pretty good, but not great. Let's think about how we can improve things …

We created a template that faithfully represented the timing of the task, so our results may not be optimal. One reason is that there is a ~4-8 sec lag between the onset of the task and the onset of this slow blood flow related response. You can account for the hemodynamic lag by time-shifting your template waveform. You can edit your expected activation template by adding zeros to the front end, and removing the same number of zeros from the back end. This has the effect of shifting the expected activation template to the right. To shift to the left, remove zeros from the front end and add them to the back end. Remember, however, you must have 150 elements in this vector.

6. Account for the hemodynamic lag by shifting the template to the right or left by adding and removing zeros from the beginning and end of the vector (again, this shifts it in time).

• Rerun the analysis with the new model (aka template)
• Observe how this changes the number of activated voxels.

As you probably (should have) observed, the statistics are VERY sensitive to getting the expected stimulus waveform right and for accounting for the hemodynamic delay.

Now examine the results with the results that generated the largest number of activated voxels.

Have you noticed how noisy some of the raw time waveforms appear - even when coming from “significant” voxels? We know that our raw MR time series are composed of signal and noise. Some of the noise is of a higher frequency than the signal, so perhaps we can get rid of some of it through the imposition of a simple low pass filter. If we suppress the noise, we should make our raw time series more like our expected waveform, and our correlations and probabilities should improve.

7. Included in your script is a bit of code that applies a very simple moving average filter. The actual operation is the line in which we calculate the mean of the time-series over time points j-len to j+len; where len = length of kernel.

timeseries = squeeze(func.data(x,y,z,1:tdim));
if(mean(timeseries) > threshold)
  if(movavg ~= 0)
    temp = timeseries;
    for j = len+1:tdim-len
      temp(j) = mean(timeseries(j-len:j+len));
    end
    timeseries(len+1:tdim-len) = temp(len+1:tdim-len);
end

Rather than create a new script to implement this moving average filter, we can modify a variable at the top of the script called movavg (see below).

• If movavg = 0 ('false' in a Matlab logical expression), then the moving average filter will NOT be applied.
• If movavg = 1 ('true' in a Matlab logical expression), then the moving average filter will be applied.

Try setting movavg = 1 and then save your script (see image below). Then run the script to visualize the effect of smoothing.

This may take a minute or two. Look at the bottom left-hand corner of the MATLAB window (near “Start”) to see if it says “Busy.” If so, it's still computing…

In fMRI studies, we are usually interested in comparing the activations evoked by different tasks - something that cannot be done with a single template like the one we've been using. In this section, you will use fmri_lab_script3_2020.m and create two templates - one for Task A and one for Task B.

1. Edit fmri_lab_script3_2025.m and modify the template1 and template2 vectors to correspond to the timing of Task A and Task B for your assigned demonstration experiment. Use what you learned in Part 3 to optimize your templates.

2. Then run the script on the same data as you did above.

3 Plot the time-series

• Make sure to highlight the functional MRI data in the Overlay list (this is the file that does not have corr in the name).
• Press command + 3
• Change the Plotting Mode to Normalised

4. Overlay your models (Load modelA, then repeat these steps to load modelB)

• Press the button to Import data series from a text file
• Select the text file you just created with your MATLAB script.
  • It will be in your data directory and named SUBJ_TASK_modelA.txt
• Click Ok in the Scaling factor pop-up window

• Note that the activations identified by template1 and template2 are color-coded in the final output.
  • Voxels that correlated with template1 (Task A) are red-yellow, whereas those that correlated with template2 (Task B) are blue-lightblue.

The motor task was run twice and the face task was usually run three times in each subject. This was done to increase the sample size and thus increase signal-to-noise. The figure below shows the time-series from a single voxel for run1 (red), run2 (blue), run3 (green), and the average all three (black). You can see that the addition of two runs helps smooth out the data and increase SNR.

1. As our final exercise, use fmri_lab_script5_2025.m to apply your two template model to the average of two runs.

• You can simply copy and paste the ideal templates (template1 and template2) from fmri_lab_script3_2025.m into fmri_lab_script5_2025.m.
• Don't forget to look at the time-series with both of your templates overlaid.

fmri_lab_scripts5_2025('SUBJ','TASK')