First-level analyses (using FSL)#

This notebook provides the first part (first-level analysis) of an introduction to multilevel models and how to implement them in FSL. The next notebook, run_level_analysis.ipynb, will continue with run-level analysis, and next week we will conclude this topic with group-level analyses. If you haven’t done so already, please go throught linux_and_the_command_line notebook first!

You actually know most about “first-level analyses” already, as it describes the process of modelling single-subject (single-run) timeseries data using the GLM, which we discussed at length in week 2. In this notebook, you’ll learn how to use FSL for (preprocessing and) first-level analyses. Note that there are several other excellent neuroimaging software packages, such as SPM (Matlab-based), AFNI, BROCCOLI, Freesurfer, and Nilearn. In fact, in week 7, there is a notebook on how to use Nilearn to fit statistical models on fMRI data.

We like FSL in particular, as it’s a mature and well-maintained package, free and open-source, and provides both a graphical interface and a command line interface for most of its tools.

What you’ll learn: after this lab, you’ll be able to …

  • visualize and inspect (f)MRI in FSLeyes;

  • set up a first-level model in FSL FEAT

Estimated time needed to complete: 2-4 hours

# Import functionality
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Using FSLeyes#

In the previous notebook about Linux and the command line, we discussed how to use built-in Linux commands and FSL-specific commands using the command line. Now, let’s focus on using a specific FSL tool called “FSLeyes”, which is a useful program to view (f)MRI files (basically any nifti-file). We can start this graphical tool from the command line with the command: fsleyes &. We append the & to the FSLeyes command because this makes sure FSLeyes will run “in the background” and thus allows us to keep using the command line (if we would not use the &, FSLeyes would run be active “in your terminal”, and we wouldn’t be able to use that anymore).

ToDo (not graded): Start FSLeyes from the command line.

This should open a new window with FSLeyes. Now, let’s open a file!

ToDo (not graded): In the menu-bar, click File, Add from file, navigate to the week_5 folder, and select the file sub-02_desc-preproc_T1w.nii.gz.

Now, FSLeyes should look similar to the figure below.


Image from the FSLeyes documentation

As you can see, FSLeyes consists of multiple views, panels, and toolbars. You should see an anatomical MRI image shown in three perspectives (sagittal, coronal, and axial) with green crosshairs in the “Ortho view”.

ToDo (not graded): Try clicking at different locations in either of the three views, which should update the crosshairs and the two other views. Note that you can also use the arrow keys on your keyboard to "scroll through" the image within a particular view, although you might experience some lag doing so. Also try to zoom in using the slider in the "Ortho toolbar".

Notice the labels on the sides of the image: P (posterior), A (anterior), I (inferior), S (superior), L (left), and R (right).

ToThink: How do you think FSLeyes knows what is left and what is right (or anterior vs. posterior or inferior vs. superior)? (Think about the lab from week 1!)

In the “Overlay toolbar”, you can adjust both the brightness and contrast of the image. One could argue that the default settings are a little bit too dark (for this scan at least).

ToDo (not graded): Try increasing the brightness and the contrast by dragging the sliders to the right. Note that the "Min." and "Max." values also change if you change the brightness/contrast. You can also set these values directly.

Importantly, while moving the location of the crosshairs, you see in the lower right corner (under the “Location panel”) that some numbers change. The left side of the panel shows you both the location of the crosshairs in “world coordinates” (relative to the isocenter of the scanner, in millimeters) and in “voxel coordinates” (starting from 0).

The right side of the panel shows you the voxel coordinates of the crosshairs together with the signal intensity, e.g., \([156\ 123\ 87]: 600774.75\). Try changing the crosshairs from a voxel in the CSF, to a voxel in gray matter, and a voxel in white matter and observe how the signal intensity changes.

ToDo (1 point): What is the size of the three dimensions of this anatomical scan? In other words, how many slices are there in the X (saggital), Y (coronal), and Z (axial) direction? Hint: check whether FSL starts counting from 1 or from 0. Fill in your answer below (by replacing the Nones).
# Replace the Nones with the correct numbers
X = None
Y = None
Z = None

raise NotImplementedError()
''' Tests the above ToDo.'''
from niedu.tests.nii.week_5 import test_how_many_slices_fsleyes    
test_how_many_slices_fsleyes(X, Y, Z)

Now, let’s check out some functional data (fMRI) in FSLeyes!

ToDo (not graded): Add the sub-02_task-flocBLOCKED_space-T1w_desc-preproc+trimmed_bold.nii.gz files from the week_5 directory. Note that this file has already been preprocessed by Fmriprep (including motion correction, signal distortion correction, and a rigid body registration to the anatomical file). By adding this file, you are creating an overlay, which you can see in the overlay list and choose to be visible or invisible by ticking the box next to it.

FSLeyes neatly plots the functional data on top of the anatomical image. Note that the resolution differs between the anatomical (\(1 \times 1 \times 1\) mm) and functional (\(2.7 \times 2.7 \times 2.97\) mm)! FSLeyes still allows you to view these images on top of each other by interpolating the images.

Now, let’s make our view a little bit nicer.

ToDo (not graded): Set the minimum value ("Min." in the "Overlay toolbar") of the fMRI file to 8000. Make sure the sub-02_task-flocBLOCKED_space-T1w_desc-preproc+trimmed_bold.nii.gz file is actually the "active" overlay (i.e., it should be highlighted in blue in the "Overlay list").

Setting a minimum value will threshold the image at this value, which will mask most of the background. Now, let’s also change the color scheme for the func file to be able to distinguish the func and anat file a little better:

ToDo (not graded): In the "Overlay toolbar", change the colormap from "Greyscale" to "Red-Yellow".

Now you can clearly see that setting a threshold did not mask all non-brain tissue, as some skull seems to still show “activity”!

ToDo/ToThink (not graded): Set your crosshair at $[X = 27, Y = 28, Z = 17]$ (in voxel coordinates). Note that you can also fill in these values in the "Location panel" directly, which will update your crosshairs automatically. In this particular "Ortho view", which two regions seem to be affected by "signal loss"? To get a better view at the underlying anatomical image, you can either adjust the opacity (transparency) of the fMRI file in the "Overlay toolbar" or toggle the fMRI file on and off in the "Overlay list" by clicking on the Eye-icon left to the highlighted sub-02_task-flocBLOCKED...bold file.

Again, just like with the anatomical scan, the “Location panel” shows you both the (voxel and world) coordinates as well as the signal intensity.

ToDo/ToThink (not graded): Scroll through your functional file. Where does the signal intensity seem to be the highest? Does this surprise you given the contrast that is usually most dominant in fMRI scans? (Tip: it helps if you decrease the "Max." value to around 90,000, set the opacity slider all the way to the right, and it's best visible in the axial view!)
ToDo (not graded): Try viewing a couple of other volumes (increase/decrease the number next to "Volume"). Visually, you probably won't see much of a difference between different volumes.
ToThink: Would you expect such small intensity differences between volumes? Why (not)?

Instead of manually adjusting the volume number, you can press the “Video” icon (next to the big “+” icon, which you can use to toggle the crosshairs on or off) to start “Movie mode” in which FSLeyes will loop through all volumes. Due to the lag on the server (or a software bug in FSLeyes, I’m not sure), however, FSLeyes is not fast enough to render the volumes, so you’ll likely see a black screen…

That said, looping through your volumes is usually a great way to visually check for artifacts (e.g. spikes) and (extreme) motion in your data! We highly recommend doing this for all your data!

Another useful way to view your fMRI data is using the “Time series” view.

ToDo (not graded): Open the Timeseries viewer (ViewTime series). This should split the window in two, with the Timeseries viewer on the bottom. Check out the time series of different voxels.
ToDo/ToThink (not graded): Go to voxel $[X=23, Y=4, Z=19]$ and look at the time series. Can you deduce what type of experimental design was used in this task? (Hint: the filename kind of gives it away ...)
ToDo/ToThink (not graded): Go to voxel $[X=26, Y=27, Z=4]$ and look at the time series. Given the location, what do you think the time series of this voxel represents?

FSLeyes has many more features than we discussed in this section (e.g., working with standard atlases/parcellations). We’ll demontrate these features throughout the rest of this lab where appropriate!


When viewing the anatomical scan in the previous section, you might have noticed that the scan still contained the subject’s skull. Usually, the skull is removed because it improves registration to the subject’s anatomical T1-scan (especially if you use BBR-based registration). This process is often called “skullstripping” and FSL contains the skullstripping-tool bet (which stands for “Brain Extraction Tool”). We’ll demonstrate how to use bet using its command-line interface.

Like all FSL tools with command-line interfaces, you can get more info on how to use it by calling the command with the flag --help.

ToDo (not graded): Call bet --help and check out the usage info.

You see that the “bet” command should be called as follows:

bet <input> <output> [options/flags]
ToDo (not graded): Skullstrip the sub-02_desc-preproc_T1w.nii.gz file using bet and name the output-file sub-02_desc-preproc+bet_T1w.nii.gz (don't use any additional flags for now).

We recommend visually checking the result of skullstripping for imperfections (e.g. part of the brain that are excluded or parts of non-brain that are included). Let’s do this in FSLeyes!

ToDo (not graded): If there are any active images in FSLeyes at this point, remove them (Overlay → Remove all). Also, close the Timeseries view (click on the cross in the upper right corner of the timeseries panel).

Then, load the non-skullstripped anatomical file (sub-02_desc-preproc_T1w.nii.gz) first. Then, overlay the skullstripped image (the sub-02_desc-preproc+bet_T1w.nii.gz file you just created). The difference between the non-skullstripped and the skullstriped image is, now, however not clear; to make this clearer, we can pick a different colormap for the skullstripped brain! For example, change the “Grayscale” colormap to “Red-Yellow”. Also, change the opacity to about 50% to better see the underlying non-skullstripped anatomical scan.

Scroll through the image to see whether bet did a good job.

ToDo (not graded): In fact, bet seemed to have removed some parts of the brain that are actually gray matter tissue (e.g., $[X=97, Y=84, Z=142]$ and $[X=154, Y=187, Z=109]$).

Also, it failed to remove some non-brain tissue. For example, what non-brain structure do you think is located at \([X=143, Y=178, Z=72]\)? And at \([X=140, Y=124, Z=44]\)?

You could, of course, just skullstrip all your participants without looking at the resulting skullstripped brain images. However, as we just saw, this may lead to suboptimal results, which may bias subsequent preprocessing steps, such as functional → anatomical registration and spatial normalization to a standard template. As such, we recommend visually inspecting the skullstrip result for each participant in your dataset!

If you find that bet did not do a proper job, you can adjust several parameters to improve the skullstripping operation. The two most useful parameters to adjust are “robust brain centre estimation” (flag -R), “fractional intensity threshold” (flag -f <value>), and the vertical gradient (flag -g <values>). We recommend always enabling “robust brain centre estimation” (which iterates the operation multiple times, which usually gives more robust results).

The proper value for fractional intensity threshold (-f) and the vertical gradient (-g), however, you need to “tune” yourself (although the default, 0.5, is a good starting value). Smaller values for fractional intensity threshold give you larger brain outline estimates (i.e., “threshold” to be used for image is lower). With the vertical gradient you can tell bet to be relatively more “aggressive” (with values < 0) or “lenient” (with values > 0) with removing tissue at the bottom (inferior side) of the brain. So, e.g., using -g -0.1 will remove relatively more tissue at the bottom of the brain.

To include these parameters, you can call bet with, for example, a “fractional intensity threshold” of 0.4, a vertical gradient of -0.1, and robust brain centre estimation enabled as follows:

bet <input> <output> -f 0.4 -g -0.1 -R
ToDo (1 point)

For this ToDo, try out some different parameters for fractional intensity threshold and vertical gradient, and enable robust brain centre estimation. Note that the output file will be overwritten if you use the same output-name (which is fine if you don’t want to save the previous result). Also, there may not be a combination of settings that creates a “perfect” skullstrip result (at least for this brain).

Once you’re happy with the segmentation, overlay it on the original anatomical file in FSLeyes (and use a different colormap). Then, make a snapshot using the camera icon and save it in your week_5 directory (name it “my_skullstripped_image.png”). After you created and saved the snapshot, add the following line to the answer-cell below:


After adding this to the text-cell below and running the cell, the snapshot should be visible in the notebook. (This allows us to check whether you actually did the ToDo; there is not necessarily a “correct” answer to this ToDo — we just want you to mess around with bet and the command line).


Note that adjusting the parameters of bet may in fact not completely solve all the imperfections in the skullstripping procedure. If you really care about accurate segmentation (e.g., when doing anatomical analyses such as voxel-based morphometry or surface-based analyses), you may consider editing your skullstrip result manually. FSLeyes has an “edit mode” (ToolsEdit mode) in which can manually fill or erase voxels based on your own judgment what is actual brain tissue and what is not.

Skullstripping is actually the only preprocessing operation you have to do manually using the command line (which might be a good idea, anyway, as you should inspect the result manually). The rest of the preprocessing steps within FSL are actually neatly integrated in the FSL tool “FEAT” (FMRI Expert Analysis Tool), which is FSL’s main analysis workhorse (implementing first-level and grouplevel analyses). We’ll discuss doing first-level and fixed-effect analyses with FEAT later in this tutorial, but first let’s talk about the dataset that we’ll use for the upcoming lab!

The “NI-edu” dataset#

So far in this course, we have used primarily simulated data. From this week onwards, we’ll also use some real data from a dataset (which we’ll call “NI-edu”) acquired at our 3T scanner at the University of Amsterdam. The NI-edu study has been specifically “designed” for this (and the follow-up “pattern analysis”) course. The experimental paradigms used for this dataset are meant to investigate how the brain processes information from faces and the representation of their subjective impressions (e.g., their perceived attractiveness, dominance, and trustworthiness). In addition to tasks that investigate face processing specifically, the dataset also includes several runs with a “functional localizer” task, which can be used to identify regions that preferentially respond to a particular object category (like faces, bodies, characters, etc.). In this course, we’ll mostly use the localizer data, because the associated task (design) is relativelly simple and generates a robust signal.

Below, we download the data needed for this tutorial, which is the anatomical data and functional data from a single run of the flocBLOCKED task (see below) from a single subject (sub-03). The data will be saved in your home directory in a folder named “NI-edu-data”.

import os
data_dir = os.path.join(os.path.expanduser("~"), 'NI-edu-data')
print("Data (+- 110MB) will be saved in %s ..." % data_dir)
print("Go get some coffee. This may take a while.\n")

# Using the CLI interface of the `awscli` Python package, we'll download the data
# Note that we only download the preprocessed data from a single run from a single subject
# This may take about 2-10 minutes or so (depending on your download speed)

# Download anatomical files
!aws s3 sync --no-sign-request s3:// {data_dir} --exclude "*" --include "sub-03/ses-1/anat/*"

# And functional files
!aws s3 sync --no-sign-request s3:// {data_dir} --exclude "*" --include "sub-03/ses-1/func/*task-flocBLOCKED*"

# And a brain mask (needed for FEAT, later)
!aws s3 sync --no-sign-request s3:// {data_dir} --exclude "*" --include "derivatives/fmriprep/sub-03/ses-1/anat/sub-03_ses-1_acq-AxialNsig2Mm1_desc-brain_mask.nii.gz"

Data (+- 110MB) will be saved in /Users/lukas/NI-edu-data ...
Go get some coffee. This may take a while.

Anyway, let’s check out the structure of the directory with the data. The niedu package has a utility function (in the submodule utils) that returns the path to the root of the dataset:

if os.path.isdir(data_dir):
    print("The dataset is located here: %s" % data_dir)
    print("WARNING: the dataset could not be found here: %s !" % data_dir)
The dataset is located here: /Users/lukas/NI-edu-data

Contents of the dataset#

Throughout the course, we’ll explain the contents of the dataset and use them for exercises and visualizations. For now, it suffices to know that the full dataset contains data from 13 subjects, each scanned twice (i.e., two sessions acquired at different days). Each subject has it’s own sub-directory, which contains an “anat” directory (with anatomical MRI files) and/or a “func” directory (with functional MRI files).

Info: The exact number of subjects, tasks, and sessions included in your dataset depends on the version of the dataset you downloaded. The "minimal" version only includes the data absolutely necessary for this course (a single subject), while the "complete" version contains all data (13 subjects; but is substantially larger in size).

We can “peek” at the structure/contents of the dataset using the listdir function from the built-in os package (which stands for “Operating System”):

import os
['derivatives', 'sub-03']

The sub-03 directory contain the “raw” (i.e., not-preprocessed) data of subject 3. Let’s check out its contents. We’ll do this using a custom function show_tree from the niedu package:

from niedu.global_utils import show_directory_tree

# To create the full path, we use the function `os.path.join` instead of just writing out
# {data_dir}/sub-03. Effectively, `os.path.join` does the same, but takes into account your
# operating system, which is important because Windows paths are separated by backslashes
# (e.g., `sub-03\ses-1\anat`) but Mac and Linux paths are separated by forward slashes 
# (e.g., `sub-03/ses-1/anat`).
sub_03_dir = os.path.join(data_dir, 'sub-03')
show_directory_tree(sub_03_dir, ignore='physio')

As you can see, the sub-03 directory contains a single session directory, “ses-1”. “Sessions” ususally refer to different days that the subject was scanned. Here, we only downloaded data from a single session (while, in fact, there is data from two sessions) in order to limit the amount of data that needed to be downloaded.

Within the ses-1 directory, there is a func directory (with functional MRI files) and an anat directory (with anatomical MRI files). These directories contain a bunch of different files, with different file types:

  • (Compressed) nifti files: these are the files with the extension nii.gz, which contain the actual MRI data. These are discussed extensively in the next section;

  • Json files: these are the files with the extension json, which are plain-text files containing important metadata about the MRI-files (such as acquisition parameters). We won’t use these much in this course.

  • Event files: these are the files ending in _events.tsv, which contain information about the onset/duration of events during the fMRI paradigm.

This particular organization of the data is in accordance with the “Brain Imaging Data Structure” (BIDS; read more about BIDS here). If you are going to analyze a new or existing MRI dataset, we strongly suggest to convert it to BIDS first, not only because makes your life a lot easier, but also because there exist many tools geared towards BIDS-formatted datasets (such as fmriprep and mriqc)!

MRI acquisition parameters#

For each participant, we acquired one T1-weighted, high-resolution (\(1 \times 1 \times 1\) mm.) anatomical scan. Our functional scans were acquired using a (T2*-weighted) gradient-echo echo-planar imaging (GE-EPI) scan sequence with a TR of 700 milleseconds and a TE of 30 milliseconds and a spatial resolution of \(2.7 \times 2.7 \times 2.7\) mm. These functional scans used a “multiband acceleration factor” of 4 (which is a “simultaneous multislice” technique to acquire data from multiple slices within a single RF pulse, greatly accellerating acquisition time).

Experimental paradigms#

In this dataset, we used two types of experimental paradigms: a “functional localizer” (which we abbreviate as “floc”) and a simple face perception task (“face”). (We also acquired a resting-state scan, i.e. a scan without an explicit task, but we won’t use that for this course.) The two different floc paradigms (flocBLOCKED and flocER) were each run twice: once in each session. Each session contained 6 face perception runs (so 12 in total).

Functional localizer task#

The functional localizer paradigm we used was developed by the Stanford Vision & Perception Neuroscience Lab. Their stimulus set contains images from different categories:

  • characters (random strings of letters and numbers);

  • bodies (full bodies without head and close-ups of limbs);

  • faces (adults and children);

  • places (hallways and houses);

  • objects (musical instruments and cars).

These categories were chosen to target different known, (arguably) category-selective regions in the brain, including the occipital face area (OFA) and fusiform face area (FFA; face-selective), parahippocampal place area (PPA; place/scene selective), visual word form area (VWFA; character-selective) and extrastriate body area (EBA; body-selective).

In the functional localizer task, participants had to do a “one-back” task, in which they had to press a button with their right hand when they saw exactly the same image twice in a row (which was about 10% of the trials).

Importantly, we acquired this localizer in two variants: one with a blocked design (“flocBLOCKED”) and one with an event-related design (“flocER”; with ISIs from an exponential distribution with a 2 second minimum). The blocked version presented stimuli for 400 ms. with a fixed 100 ms. inter-stimulus interval in blocks of 12 stimuli (i.e., 6 seconds). Each run started and ended with a 6 second “null-blocks” (a “block” without stimuli) and additionally contained six randomly inserted null-blocks. The event-related version showed each stimulus for 400 ms. as well, but contained inter-stimulus intervals ranging from 2 to 8 seconds and no null-blocks apart from a null-block at the start and end of each run (see figure below for a visual representation of the paradigms).


These two different versions (blocked vs. event-related) were acquired to demonstrate the difference in detection and estimation efficiency of the two designs. In this course, we only use the blocked version.

Face perception task#

In the face perception task, participants were presented with images of faces (from the publicly available Face Research Lab London Set). In total, frontal face images from 40 different people (“identities”) were used, which were either without expression (“neutral”) or were smiling. Each face image (from in total 80 faces, i.e., 40 identities \(\times\) 2, neutral/smiling) was shown, per participant, 6 times across the 12 runs (3 times per session). We won’t discuss the data from this task, as it won’t be used in this course.

Using FSL FEAT for preprocessing and first-level analysis#

Most univariate fMRI studies contain data from a sample of multiple individuals. Often, the goal of these studies is to (1) estimate some effect (e.g., \(\beta_{face} >0\)) at the group-level, i.e., whether the effect is present within our sample, and subsequently, taking into account the variance (uncertainty) of this group-level effect, (2) to infer whether this effect is likely to be true in the broader population (but beware of the tricky interpretation of null-hypothesis significance tests).

A primer on the “summary statistics” approach#

When you have fMRI data from multiple participants, you are dealing with inherently hierarchical (sometimes called “multilevel”, “mixed”, or simply “nested”) data: you have data from multiple participants, which may be split across different runs or sessions, which each contain multiple observations (i.e., volumes across time).

Often, people use so called “mixed models” (or “multilevel models”) to deal with this kind of hierarchical data. These type of models use the GLM framework for estimating effects, but have a special way of estimating the variance of effects (i.e., \(\mathrm{var}[c\hat{\beta}]\)), which includes variance estimates from all levels within the model (subject-level, run-level, and group-level). In these mixed models, you can simply put all your data (\(y\) and \(\mathbf{X}\)) from different runs, sessions, and subjects, in a single GLM and estimate the effects of interest (and their variance).

Unfortunately, this “traditional” multilevel approach is not feasible for most univariate (whole-brain) analyses, because these models are quite computationally expensive (i.e., they take a long time to run), so doing this for all > 200,000 voxels is the brain is not very practical. Fortunately, the fMRI community came up with a solution: the “summary statistics” approach, which approximates a proper full multilevel model.

This summary statistics approach entails splitting up the analysis of (hierarchical) fMRI data into different steps:

  • Analysis at the single-subject (time series) level

  • Analysis at the run level (optional; only if you have >1 run!)

  • Analysis at the group level

The GLM analyses you did in week 2-4 were examples of “first-level” analyses, because they analyzed the time series of voxels from a single subject. This is the first step in the summary statistics approach. In this section, we’re going to demonstrate how to implement first-level analyses (and preprocessing) using FSL FEAT.

Later in this lab, we’re going to demonstrate how to implement analyses at the run-level. Next week, we’ll discuss group-level models. For now, let’s focus on how to run first-level analyses in FSL using FEAT.

First-level analyses (and preprocessing) in FEAT#

Alright, let’s get started with FEAT. You can start the FEAT graphical interface by typing Feat & in the command line (the & will open FEAT in the background, so you can keep using your command line).*

* Note that the upcoming ToDos do not have test-cells; at the very end of this section, you’ll save your FEAT setup, which we’ll test (with only hidden tests) afterwards.

ToDo (not graded): Open FEAT (Feat &). Make sure you are in your home folder (which you can navigate to by typing cd) when you do.

The “Data” tab#

After calling Feat & in the command line, the FEAT GUI should pop up. This GUI has several “tabs”, of which the “Data” tab is displayed when opening FEAT. Below, we summarized the most important features of the “Data” tab.


ToDo (not graded): In this section, we'll analyze a single run of the flocBLOCKED paradigm of a single participant (sub-03). The data from this participant is located in the following directory (not in your public directory jupyter-nim_XX/Public/):

We are going to analyze the flocBLOCKED data from the first session (ses-1): sub-03/ses-1/sub-03_ses-1_task-flocBLOCKED_bold.nii.gz.

sub_03_file = os.path.join(sub_03_dir, 'ses-1', 'func', 'sub-03_ses-1_task-flocBLOCKED_bold.nii.gz')
print("We are going to analyse: %s" % sub_03_file)
We are going to analyse: /Users/lukas/NI-edu-data/sub-03/ses-1/func/sub-03_ses-1_task-flocBLOCKED_bold.nii.gz
ToDo (0.5 points)

Let’s fill in the Data tab in FEAT. Click on “Select 4D data”, click on the “folder” icon, and navigate to the file (sub-03_ses-1_task-flocBLOCKED_bold.nii.gz). Once you’ve selected the file, click on “OK”. Note: within FEAT, you can click on .. to go “up” one folder. This way, you can “navigate” to the location of the file that you want to add.

Then, set the output directory to your home folder + sub-03_ses-1_task-flockBLOCKED. Note: this directory doesn’t exist yet! By setting the output directory to this directory, you’re telling FSL to create a directory here in which the results will be stored.

After selecting the 4D data, the “Total volumes” and “TR” fields should update to 325 and 0.700000, respectively. You can leave the “Delete volumes” field (0) and high-pass filter (100) as it is (no need to remove volumes and 100 second high-pass is a reasonable default).

The “Pre-stats” tab#

Now, let’s check out the “pre-stats” tab, which contains most of the preprocessing settings. Most settings are self-explanatory, but if you want to know more about the different options, hover you cursor above the option (the text, not the yellow box) for a couple of seconds and a green window will pop up with more information about that setting!

Note that the checkbox being yellow means that the option has been selected (and a gray box means it is not selected).

Ignore the “Alternative reference image”, “B0 unwarping”, “Intensity normalization”, “Perfusion subtraction”, and “ICA exploration” settings; these are beyond the scope of this course!

ToDo (0.5 points)

Let’s fill in the Pre-stats tab! Select the following:

  1. Select “MCFLIRT” for motion-correction (this is the name of FSL’s motion-correction function);

  2. Set “Slice timing correction” to “None” (we believe that slice-time correction does more harm than good for fMRI data with short TRs, like ours, but “your mileage may vary”);

  3. Turn on “BET brain extraction” (this refers to skullstripping of the functional file! This can be done automatically, because it’s less error-prone than skullstripping of anatomical images);

  4. Set the FWHM of the spatial smoothing to 4 mm;

  5. Turn off “intensity normalization” (you should never use this for task-based fMRI analysis);

  6. Turn on “highpass filtering”; leave “perfusion subtraction” unchecked;

  7. Leave “ICA exploration” unchecked

The “Registration” tab#

Alright, let’s check out the “Registration” tab, in which you can setup the co-registration plan (from the functional data to the T1-weighted anatomical scan, and from the T1-weighted anatomical scan to standard space; we usually don’t have an “expanded functional image”, so you can ignore that).

ToDo (0 points): Enable the "Main structural image" option. You should see an empty field (where you need to fill in the path to the skullstripped structural file) and below it a "search" option and the particular registration method. We'll leave it to "Normal search" and BBR ("Boundary-Based Registration", a particular functional-to-T1 registration technique).

Now, for the “Main structural image”, you actually need a high-resolution T1-weighted anatomical file that has already been skullstripped. You could do that with, for example, bet for this particular partipant. Below, we do this slightly differently, by applying an existing mask, computed by Fmriprep, to our structural image. We use the Python package Nilearn for this (you will learn more about this package in week 7)! No need to fully understand the code below.

from nilearn import masking, image

# Define anatomical file ...
anat = os.path.join(sub_03_dir, 'ses-1', 'anat', 'sub-03_ses-1_acq-AxialNsig2Mm1_T1w.nii.gz')

# ... and mask
mask = os.path.join(data_dir, 'derivatives', 'fmriprep', 'sub-03', 'ses-1', 'anat', 'sub-03_ses-1_acq-AxialNsig2Mm1_desc-brain_mask.nii.gz')

# Affine of mask is different than original anatomical file (?),
# so let's resample the mask to the image
mask = image.resample_to_img(mask, anat, interpolation='nearest')

# Set everything outside the mask to 0 using the apply_mask/unmask trick
anat_sks = masking.unmask(masking.apply_mask(anat, mask), mask)

# And finally save the skullstripped image in the same directory as the
# original anatomical file
out = anat.replace('.nii.gz', '_brain.nii.gz')
if not os.path.isfile(out):
ToDo (0.5 points): Select the skullstripped file of sub-03 (the one ending in _brain.nii.gz in the anatomy subfolder) for the "Main structural image". In fact, FSL also needs the original (non-skullstripped) T1-file, but it will assume that the original file is in the same directory and has the same name, but without the _brain suffix (which is the case for our file).

You can leave the file under “Standard space” as it is (this is the MNI152 2mm template, one of the most used templates in the fMRI research community). Turn on the “Nonlinear” option in the “Standard space” box and do not change the default options of non-linear registration.

The “Stats” tab#

Almost there! Up next is the “Stats” tab, in which you specify the design (\(\mathbf{X}\)) you want to use to model the data.

ToDo (0.5 points): First, make sure "FILM prewhitening" is turned on ("FILM" is just the name of FSL's prewhitening implementation). Then, select "Standard Motion Parameters" in the box below, which instructs FEAT to use the six motion parameters as nuisance predictors in the design. You don't have to do anything with the "Voxelwise Confound List" and "Add additional confound EVs" options.

Then, to specify our predictors-of-interest, we need to click the “Full model setup” button. This creates a new window in which we can create our predictors; see the image below for an explanation of the options.

The most often-used option to define the event onsets and durations is the “Custom (3 column format)” option. This allows you to select a standard text-file with three columns (separated by tabs), in which the first column indicates the event onsets (in seconds), the second column indicates the event duration (in seconds), and the third column indicates the event “weight” (which may vary in case of parametric designs; for non-parametric designs, you should set the weight-value to 1). Different rows define different events of a particular condition.

Our data only has a BIDS-compatible “events” file, which contain the onset, duration, and trial type of trials from all conditions, while FSL requires a separate events file per condition (and without a header). Below, we run a function that creates these separate condition files for us:

from niedu.utils.nii import create_fsl_onset_files
ToDo (1 point):

Let’s set up the design! For our flocBLOCKED run, we should model each of the image conditions (“face”, “body”, “place”, “character”, and “object”) as separate predictors, assuming that they have a different effect on voxels (in different brain regions). Let’s model the responses (button presses) as well.

So, create six predictors with the names “face”, “body”, “place”, “character”, “object”, and “response” (in this exact order). To save you some work, we have already created the onset-files associated with these six different conditions (which end in: *condition-{name of condition}_events.txt).

For each predictor:

  • set the convolution to “Double-Gamma HRF”

  • turn off the “Add temporal derivative” and “Orthogonalise” options

  • leave the “Apply temporal filtering” option checked (this applies the high-pass filter to the design, \(X\), as well).

ToThink (1 point): Here, we chose not to include temporal derivatives, because we believe that, given the experimental design, including them is unlikely to lead to a better model fit. Why do you think this might be the case?


If you did it correctly, you can press the “View design” button and a new window with the design should pop up (with the different predictors as squiggly lines from top top bottom). The red line of the left side of the design matrix represents the (length of the) high-pass filter.

ToThink (1 point): If you look closely at the predictors (columns) of the visualized design matrix, you'll notice that the first and last couple of time points are not "straight" (i.e., at zero) but are slightly increasing or decreasing. What causes this effect? (Hint: the answer can be found in the 'temporal preprocessing' tutorial!)


Now let’s create some contrasts! In the “Contrasts and F-tests” tab, you can specify the contrasts that you want to test. Importantly, FSL makes a difference between “Original EVs” and “Real EVs”. Original EVs are the predictors-of-interest that you defined in the “EVs” tab. Real EVs are the original EVs + the extra predictors that are generated due to the specific HRF-model. For example, if you’d choose a gamma basis set instead of a double-gamma HRF, the original canonical HRF of each predictor would be “Original EVs” while all predictors making up the HRF-components (canonical and temporal and dispersion derivatives) would be the “Real EVs”. In our case, the original and real EVs are the same.

To create contrasts (and F-tests, i.e., collections of contrasts), click on the “Contrasts & F-tests” tab. Here, you define the actual contrast vectors. For each contrast, you can give it a name (under “Title”) and set the contrast weight for each predictor (or “EV”). To create F-tests, set the number of F-tests that you want to do (i.e., increase the number next to the “F-tests” field), and check, per F-test, which contrasts to want to include in them. For example, if you want to include contrasts 1 and 2 in F-test, make sure that these checkboxes are selected under “F1” (or “F2”, “F3”, etc.).

ToDo (4 points, 0.5 per contrast/F-test)

Given the six original EVs (“face”, “body”, “place”, “character”, “object”, and “response”), define the contrasts corresponding to the following null (\(H_{0}\)) and alternative (\(H_{a}\)) hypotheses (give them whatever name you want, but adhere to this order):

  1. \(H_{0}: \beta_{face} = 0\)
    \(H_{a}: \beta_{face} > 0\)

  2. \(H_{0}: \beta_{place} = 0\)
    \(H_{a}: \beta_{place} > 0\)

  3. \(H_{0}: \beta_{response} = 0\)
    \(H_{a}: \beta_{response} > 0\)

  4. \(H_{0}: 4\times\beta_{face} - \beta_{body} - \beta_{place} - \beta_{character} - \beta_{object} = 0\)
    \(H_{a}: 4\times\beta_{face} - \beta_{body} - \beta_{place} - \beta_{character} - \beta_{object} > 0\)

  5. \(H_{0}: 4\times\beta_{place} - \beta_{body} - \beta_{face} - \beta_{character} - \beta_{object} = 0\)
    \(H_{a}: 4\times\beta_{place} - \beta_{body} - \beta_{face} - \beta_{character} - \beta_{object} > 0\)

  6. \(H_{0}: \beta_{face} = \beta_{body} = \beta_{place} = \beta_{character} = \beta_{object} = 0\)
    \(H_{a}: (\beta_{face} > 0)\ \&\ (\beta_{body} > 0)\ \&\ (\beta_{place} > 0)\ \&\ (\beta_{character} > 0)\ \&\ (\beta_{object} > 0)\)

  7. \(H_{0}: \beta_{face} = \beta_{body} = \beta_{place} = \beta_{character} = \beta_{object} = 0\)
    \(H_{a}: (\beta_{face} < 0)\ \&\ (\beta_{body} < 0)\ \&\ (\beta_{place} < 0)\ \&\ (\beta_{character} < 0)\ \&\ (\beta_{object} < 0)\)

Also, define the following F-test (where the vertical line stands for “or”):

  1. \(H_{0}: \beta_{face} = \beta_{place} = 0\)
    \(H_{a}: (\beta_{face} \neq 0)\ |\ (\beta_{place} \neq 0)\)

After you’re done defining the contrasts and F-test, you can press “View design” again, which will now also include the contrasts/F-tests that you defined. Moreover, if you click on “Efficiency”, a new window opens with information about the effiency of the design. On the left, a correlation matrix of the predictors is shown, where darker off-diagonal cells means lower correlation between predictors (which is a good thing!). On the right, the diagonal of the matrix shows you the “eigenvalues” of the singular value decomposition of the design matrix, which should decrease from top-left to bottom-right. In general, the brighter the values, the more efficient the design matrix. In the lower-right corner, it lists the estimated effect size (as % signal change) required for a significant effect. These estimates depend on your design efficiency (\(\mathbf{c}(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{c}^{T}\)) and the assumed level of noise (\(\hat{\sigma}^2\)): the lower these values, the better!

ToDo (ungraded): click "Done". If you implemented the contrasts/F-test correctly, you shouldn't see error messages, and you'll see the design matrix pop up again.

The “Post-stats” tab#

Alright, the last tab! The “Post-stats” tab contains options to threshold the resulting statistics-maps. The only interesting options for our purposes are those in the “Thresholding” box. Technically, we don’t need to threshold any of our statistics, because we are only interested in inference at the group-level (not the first-level analysis; we’re not aiming to infer anything from the results of a single-subject!). But using a p-value threshold of 0.05 “uncorrected” (which means: not corrected for multiple comparisons — a topic we’ll discuss next week!) is usually used to detect any potential problems during the experimental, preprocessing, and analysis process. For example, if we don’t see significant visual cortex activation in response to a task with visual stimuli at \(p < 0.05\) (uncorrected), then it is very likely that something went wrong somewhere in the preprocessing/analysis pipeline (e.g., stimulus onsets are incorrect).

ToDo (0.5 points): Set the thresholding type to "Uncorrected" and the p-value to 0.005 (not the default 0.05).

Alright, we’re done setting up FEAT for our first-level analysis! One last thing: let’s save the setup (which we’ll check and grade).

ToDo (not graded, but necessary to grade the other things): Click the "Save" button in FEAT's main menu and save the setup in your week_5 directory with the name setup_feat.fsf.

We actually ran the first-level analysis for this participant already, so you don’t have to do this yourself! We’ll check it out in the next section. (The following cells all test your setup_feat.fsf file; as they all have hidden test, you may ignore them.)

''' Tests your setup_feat.fsf file: initial check + Data tab'''
import os
import os.path as op

par_dir = op.basename(op.dirname(op.abspath('.')))
if par_dir != 'solutions':  # must be a student-account
    fsf = 'setup_feat.fsf'
    if not op.isfile(fsf):
        raise ValueError("Couldn't find a 'setup_feat.fsf' file in your week_5 directory!")
print("Only hidden tests!")
ValueError                                Traceback (most recent call last)
/var/folders/0c/2qv1yrms7f56lxzzw28nt_440000gn/T/ipykernel_15549/ in <cell line: 6>()
      7     fsf = 'setup_feat.fsf'
      8     if not op.isfile(fsf):
----> 9         raise ValueError("Couldn't find a 'setup_feat.fsf' file in your week_5 directory!")
     11 print("Only hidden tests!")

ValueError: Couldn't find a 'setup_feat.fsf' file in your week_5 directory!
''' Tests your setup_feat.fsf file: Pre-stats tab'''

print("Only hidden tests!")
Only hidden tests!
''' Tests your setup_feat.fsf file: Registration tab'''

print("Only hidden tests!")
Only hidden tests!
''' Tests your setup_feat.fsf file: Stats tab (filt) '''

print("Only hidden tests!")
Only hidden tests!
''' Tests your setup_feat.fsf file: Stats tab (EVs) '''

print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 1)'''

print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 2)'''

print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 3)'''

print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 4)'''

print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 5)'''

print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 6)'''

print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab, (con 7)'''

print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Stats tab (F-test)'''

print("Only hidden tests!")
Only hidden tests!
''' Tests your feat_setup.fsf file: Post-stats tab'''

print("Only hidden tests!")
Only hidden tests!
Note: No need to actually run this first-level analysis — we did it already for you!

Making sense of FEAT output#

When FEAT is running the analyis, it will create a results-directory (at the location you specified) ending in .feat, which contains all the information and files necessary for the analysis as well as the results. Fortunately, it also generates an informative report.html file with a nice summary of the preprocessing and analysis results. We already ran a first-level analysis like with the setup that you specified in the previous section, which we download below:

print("Starting download of FEAT directory (+- 285MB) ...\n")
!aws s3 sync --no-sign-request s3:// {data_dir} --exclude "*" --include "derivatives/fsl/sub-03/flocBLOCKED/ses-1.feat/*"
Starting download of FEAT directory (+- 285MB) ...

The downloaded data is located here:

feat_dir = os.path.join(data_dir, 'derivatives', 'fsl', 'sub-03', 'flocBLOCKED', 'ses-1.feat')
ToDo (ungraded): Navigate to the FEAT directory using the command line (with cd!) and, once you're in the right directory, open the report.html file with:

firefox report.html

Checking the report.html file#

The “Registration” tab shows you several figures related to the (two-step) co-registration process.

ToDo (ungraded): Check out the Registration tab.

In each figure, the red outline shows the target (i.e., where the data is registered to). Each figure contains two subplots, which show the co-registration in different “directions”. The reason they plot both “directions” is that a single “direction” may give you a little information about whether the registration completed without problems.

The first image shows you the result from the complete registration process: from our functional file (“FMRI”) to the MNI152 template (“standard space”). Here, the target (in red) is the outline of the MNI152 template. The other figures show the intermediate results in the two-step co-registration process. The second figure (“Registration of example_func to highres”) shows you first registration step: the BBR-based registration of (a reference volume of) our functional data to the high-resolution T1-weighted scan. Note that, in the upper row, the red outline represents the white matter boundaries of our T1-weighted scan, because this is wat BBR uses as the registration target (rather than the entire volume).

The third image show you the second registration step: the affine + non-linear registration of our high-resolution T1-weighted scan to the MNI152 template. The fourth image shows you the result of the entire two-step registration procedure (functional → template; essentially the same as the first image).

ToThink (ungraded): We didn't use distortion correction in this particular preprocessin/analysis pipeline. From the images within the Registration page, where in the brain does geometric distortion cause the most severe registration inaccuracies? (Hint: it's best visible in the first image.)
ToDo/ToThink (ungraded): Now check out the Pre-stats tab. This primarily shows you the estimated motion parameters (translations, rotations) from this run. Can you deduce from the translation/rotation plots which volume was chosen as the reference volume for motion correction?
ToDo (ungraded): Check out the Stats tab. Notice that the 6 motion parameters (indicated by the name conf, for 'confounding' variables) were added to the design by FEAT. (Also, we cut off the part of the design figure with the contrast definitions — otherwise the previous ToDos would be a bit too easy.)
ToDo (not graded): Check out the Post-stats tab. This shows you the results of the different contrasts (and F-tests) by plotting back the parameter estimates (here: "zstats", which are just t-values converted to z-values), which are thresholded at $p < 0.005$ (uncorrected for multiple comparisons), on the brain.

This overview is an excellent way to check whether the results are roughly what you’d expect. For example, you’d expect to see clear activation in visual cortex for our the contrasts against baseline (zstat1: face > 0 and zstat2: place > 0 and zstat6: stim).

Apart from the z-statistic maps, FSL shows you for each contrast the time series of the voxel that has the highest z-value. The plots show the data (i.e., \(y\); in red), the model-fit (i.e., \(\hat{y}\); in blue/purple), and the “partial model fit” from that contrast (i.e., \(X_{c}\hat{\beta}_{c}\) for that particular contrast \(c\); in green).

Checking out the results from FEAT (in FSLeyes)#

While the report.html file provides a nice summary of the results, it’s also good to visualize and look at output of FEAT is more detail. First, we’ll walk you through the contents of the FEAT directory and then we’ll take a look at the results in FSLeyes.

ToDo (ungraded): Open (a new tab in) a terminal and navigate to the FEAT directory. Then, inspect the contents using ls.

You should see a large set of different files: text-files, nifti-files, and subdirectories. Many of these files represent intermediary output from FEAT or data/images necessary for the report*.html files (which we inspected earlier). The most important subdirectories are the reg and stats directories. These contain the registration parameters and GLM results, respectively, that we later need for our fixed-effects and group-level analyses.

ToDo (ungraded): Navigate to the reg directory and list its contents (using ls).

The .mat files contain affine matrices for the different registration steps: from functional to the high-resolution anatomical scan (example_func2highres.mat) and from the high-resolution anatomical scan to the standard template (highres2standard.mat) (and the inverse matrices: highres2example_func.mat and standard2example_func.mat). In addition, it contains the “warpfield” (files ending in _warp.nii.gz) computed by the non-linear registration step, i.e., how much each voxel needs to be moved in each direction (\(X, Y, Z\)).

ToDo (ungraded): Print the contents of the example_func2highres.mat file to the terminal using cat. Does it look like what you expected?
ToDo (ungraded): Navigate into the stats folder (do you remember how to go "up" one directory?) and check the contents (with ls).

Basically, all the nifti files you see here are 3D files (“brain maps”), representing different parts of the t-value formula, as is visualized in the figure below:

In fact, to show you that the these files are simply parts of the \(t\)-value formula, let’s calculate a t-stat map ourselves using the cope and varcope files.

Assignment (2 points):

Below, we load the cope1.nii.gz (containing the \(\mathbf{c}\hat{\beta}\) values per voxel for contrast 1) and varcope1.nii.gz (containing the \(\hat{\sigma}^{2}\mathbf{c}(\mathbf{X}^{T}\mathbf{V}^{-1}\mathbf{X})^{-1}\mathbf{c}^{T}\) values per voxel for contrast 1). Calculate the corresponding 3D tstat map (containing all t-values per voxel for contrast 1) using numpy array operations, so no for loop! Store this in a variable named tstat1.

Hint: all voxels outside the brain (in both the cope1.nii.gz and the varcope.nii.gz files) will have the value 0. When calculating the t-values, this will lead to nan values (“not a number) because you’re dividing 0 by 0 (which will give a RuntimeWarning). Make sure set these nan values back to 0 (the function np.isnan might be useful ….).

# Implement your ToDo here!
import os
import nibabel as nib
import numpy as np

cope1_path = os.path.join(feat_dir, 'stats', 'cope1.nii.gz')
cope1 = nib.load(cope1_path).get_fdata()
varcope1_path = os.path.join(feat_dir, 'stats', 'varcope1.nii.gz')
varcope1 = nib.load(varcope1_path).get_fdata()

raise NotImplementedError()
''' Tests the above ToDo. '''
# Some (visible) tests to check minor mistakes
if np.nan in tstat1:
    raise ValueError('There are still nans in your data!')

from niedu.tests.nii.week_5 import test_tvalue_from_fsl
test_tvalue_from_fsl(tstat1, feat_dir)

Alright, now let’s continue with the more interesting part: visualizing the results.

ToDo (ungraded): If there are any active overlays, close them (Overlay → Remove all). Then, open the entire FEAT directory by clicking: File → Add from directory → navigate to the FEAT directory and click "Open". Note that you can use the "backspace" key to go up one directory.

FSLeyes will automatically open the filtered_func_data.nii.gz (the preprocessed functional data) file in the Ortho view.

ToDo: FSLeyes has a specialized viewing "mode" for inspecting FEAT directories. To activate this, go to View → Layouts → FEAT mode. Activating this layout will additionally open the time series viewer and "Cluster browser" on the right. You may close the Cluster browser for now.

The time series viewer now shows you both the data (\(y\)) and the “full model fit” (\(\hat{y}\),) for the voxel highlighted by the crosshairs.

ToDo/ToThink: Go to voxel $[31, 16, 22]$. You'll see that this voxel nicely follows our blocked design, but doesn't seem to differentiate between our conditions (face vs. house vs. place, etc.). Does this surprise you, given the particular brain area that this voxel is part of?

Let’s load some statistic files to see the (whole-brain) effects of a particular contrast.

ToDo (ungraded): Let's view the $z$-values from the first contrast ($\beta_{face} > 0$): File → Add from file → go to the stats directory → select zstat1.nii.gz.
ToThink (ungraded): Usually, people show brain maps of $z$-values instead of $t$-values in their figures. Why do you think this is the case?
ToDo (ungraded): To make the data easier to interpret, set the colormap of zstat1 to "Red-Yellow", set the minimum value to 2.6 (this approximately corresponds to $p<0.005$), and set the Opacity slider to about 50%. (Note that you could also have loaded the thresh_zstat1.nii.gz file from the main FEAT directory, which is already thresholded using the threshold that you set in the post-stats tab of the FEAT GUI.)

In the time series viewer, we can also add additional time series, such as the fitted “COPE” corresponding to our first contrast. The fitted “COPE” basically represents the predicted time series for a particular contrast, i.e., \(\mathbf{X}\mathbf{c}^{T}\mathbf{c}\hat{\beta}\). (Usually, however, the intercept added to the predicted time series regardless of whether it’s included in the contrast.)

For example, suppose that we want to compute the fitted COPE corresponding to the \(\beta_{face} > 0\) contrast, we would do the following:

# Let's load in the actual design matrix from the FEAT dir
dm_path = os.path.join(feat_dir, 'design.mat')
X = np.loadtxt(dm_path, skiprows=5)
X = X[:, :6]  # remove motion preds for clarity of example

# Define face > 0 contrast vector
# Add new axis to make sure matrix product works
# Order: face, body, place, character, object, response
c = np.array([1, 0, 0, 0, 0, 0])[np.newaxis, :]

# Create some hypothetical beta vector
est_betas = np.array([3.3, 2.5, 1.4, -2.1, 0.8, 0.1])
fitted_cope_face = X @ c.T @ c @ est_betas

plt.figure(figsize=(15, 3.5))
plt.xlim(0, X.shape[0])
plt.ylabel("Activity (A.U.)", fontsize=20)
plt.xlabel("Time (volumes)", fontsize=20)
plt.title(r"Fitted COPE for $\beta_{face} > 0$", fontsize=25)
plt.ylim(-2.5, 3.5)
ToDo (ungraded): Change the est_betas vector in the cell above to a couple of other contrasts (e.g., $\beta_{character} > 0$) and observe how the fitted COPE time series changes.
ToDo (ungraded): Back to FSLeyes. To add the fitted COPE time series of the first contrast, click the "wrench" icon in the time series panel, go to Time series settings for ses-1: filtered_func_data, and check the box next to Plot COPE1 fit (face).

In the time series viewer, you should see a new line corresponding to “COPE1 fit: face”.

ToDo (ungraded): Set your crosshairs at $[45, 17, 19]$. Notice that this voxel's $z$-value is very high: 15.39! However, if you look at the COPE1 fit, you'll notice that — while this voxel indeed seems to activate for faces — it does not selectively activate for faces (it seems also to activate in response to the other conditions).

This is exactly the reason why we also defined the \(4\times\beta_{face} - \beta_{body} - \beta_{place} - \beta_{character} - \beta_{object} > 0\) contrast (contrast nr. 4 in our analysis)! This contrast allows us to identify regions in the brain that are (relatively more) specific for faces than for other types of images. In fact, in the context of our functional localizer task, you could argue that the \(\beta_{face} > 0\) contrast is largely useless.

ToDo (ungraded): Remove the zstat1 file from the overlays (Overlay → Remove). Then, open the already thresholded thresh_zstat4.nii.gz from the ses-1.feat subfolder (File → Add from file → select thresh_zstat4.nii.gz → Open). Change to colormap to "Red-Yellow".

Many voxels throughout the brain seem to specifically activate in response to images of faces (but beware: these brain images are not yet “properly” thresholded, so many of these voxels might actually be false positives, but that’s the topic of next week!). Two cluster of voxels, however, seem especially strongly activated: one centered at \([24, 19, 21]\) and another one centered at \([23, 30, 13]\). Note that these clusters are present in both hemispheres, but the cluster in the right hemisphere is clearly more strongly activated, which is consistent with the apparent right-hemisphere bias in face processing (see e.g., Kanwisher et al., 1997).

ToThink (ungraded): Which functional regions do you think that these two clusters represent? (Hint: check out the Face perception wiki page.)
ToDo (ungraded): Go to the voxel at location $[23, 30, 13]$ and notice that the voxel signal seems to respond (relatively) more selectively to faces than to other conditions. But, at the risk of interpreting noise rather than signal, you can also see that this voxel is seemingly not *completely* specific to faces (e.g., there is some apparent face-unrelated activity around volume 20-40). Indeed, whether there are truly face-specific regions in the brain at all is still a matter of debate (see e.g., Tarr & Gauthier, 2000).
ToDo (ungraded): Remove the thresh_zstat4.nii.gz file (or hide this overlay by clicking on the eye-icon in the Overlay list). Then, add the thresh_zstat3.nii.gz file, corresponding to the $\beta_{response} > 0$ contrast, and set the color map to "Red-Yellow" (or any other colormap that you prefer).
ToThink (1 point): Based on this brain map, you could with reasonable confidence conclude that the participant did not completely do the task correctly (or was instructed incorrectly). Explain why.


ToDo (ungraded): Remove the thresh_zstat3.nii.gz file (or hide this overlay). Then, add the thresh_zstat7.nii.gz file, corresponding to the $\beta_{face} = \beta_{body} = \beta_{place} = \beta_{character} = \beta_{object} < 0$ contrast, and set the color map to "Red-Yellow" (or any other colormap that you prefer).

This contrast shows you which voxels are more active during periods in which stimuli were absent than when stimuli were presented (regardless of condition). The pattern of voxels that are often found to be more active during “rest” than during (external) “stimulation” has been named the default mode network (DMN), which may reflect daydreaming, episodic or autobiographical information processing, or self reflection (amongst others; the exact “function” of the DMN is still unclear).

ToDo/ToThink (1 point): Go to voxel $[35, 62, 11]$ and add the time series of COPE6 (reflecting the contrast $\beta_{face} = \beta_{body} = \beta_{place} = \beta_{character} = \beta_{object} > 0$) and COPE7 (reflecting the "rest/DMN" contrast). You'll see that these time series ($c\hat{\beta}$) are exactly the same! Explain why this is the case.


ToDo (not graded): Finally, let's check out the data that we couldn't explain: the residuals ($y - \hat{y}$) and the noise term ($\hat{\sigma}^{2}$). First, close all overlays (Overlay → Remove all). Then, open the res4d.nii.gz file (File → Open from file → select res4d.nii.gz). The time series viewer should now show the residual time series from the voxel highlighted by your crosshairs. Also add the sigmasquareds.nii.gz file (i.e., the $\hat{\sigma}^{2}$ term for every voxel; File → Add from file → select sigmasquareds.nii.gz). Set the colormap for the "sigmasquareds" overlay to "Red-Yellow" and the Max. value to 40000.

This sigmasquareds map (and the underlying “raw” residual time series) can given you a good overview of the noise in your data and potential violations of GLM (e.g., unequal or correlated error variance). This is a great way to check for unexpected issues (e.g., artifacts) with your data, that will, if unmodelled, end up in your noise term.

ToDo (not graded): If you scroll through the sigmasquareds map, you'll notice that some regions contain more "unmodelled noise" than others. For example, if you go to voxel $[37, 23, 2]$, located in the brain stem, you can see that the corresponding residual time course contains a strong periodic ("oscillating") component. Because our scan has a relatively fast TR (0.7 seconds), this time course represents most likely signal related to pulsatile cerebral bloodflow caused the participant's heartbeat.

This particular voxel lies most likely within the Circle of Willis: a “hub” of multiple arteries supplying blood to the brain. As BOLD MRI measures the (relative) amount of oxygenated blood in the brain, it is no surprise that arteries (and veins!) produce a strong heartbeat-related signal, which will, if left unmodelled, make up much of the unexplained variance. Indeed, in this sigmasquareds map, most of voxels with high values (i.e., those that are relatively yellow) are located within or near arteries or veins. For example, you can clearly “see” the sagittal sinus (see image below).

From Wikimedia commons (link)

Apart from the veins/arteries, the cerebrospinal fluid (CSF) also seem to contain a lot of unmodelled noise. In fact, most of this unmodelled variance also likely represents the pulsative cerebral bloodflow, as shown (exaggerated) in the image below:


From Wikipedia,
ToDo (not graded): Go to voxel $[37, 23, 2]$ and press the big "+" icon in the header of the Time series panel. This will "pin" this particular time series in the viewer. Click the wrench icon and under Plot settings for res4d, change the Colour to red (which will change the color of the currently highlighted time series line).

Now, move your crosshair across different voxels with high values (i.e., relatively yellow ones), for example voxel \([37, 48, 17]\) (located in the CSF). You can see that this voxel’s residual time series is highly (but not perfectly) correlated to the original voxel in the Circle of Willis, the common cause being the pulsative cerebral blood flow. (Note, though, that different voxels may show time series that look slightly shifted in time relative to the original Circle of Willis voxel, especially in more superior regions of the brain, which is caused by different slice onsets!)

Also note that the absolute signal amplitude of voxels within veins/arteries is much higher than voxels in gray matter!

Alright, time for something else: let’s delve into how you can “combine” run or session data in the next notebook (run_level_analyses.ipynb).

Tip! Before handing in your notebooks, we recommend restarting your kernel (KernelRestart & Clear Ouput) and running all your cells again (manually, or by CellRun all). By running all your cells one by one (from "top" to "bottom" of the notebook), you may spot potential errors that are caused by accidentally overwriting your variables or running your cells out of order (e.g., defining the variable 'x' in cell 28 which you then use in cell 15).