Retinotopic Mapping Tutorial
- Stimulus and Scan
- Preprocessing the EPIs
- pRF Models
- Surfaces and Parameter Visualization
- Atlas Generation
- Bayesian Inference of Retinotopic Maps
Retinotopic mapping experiments are very common in visual fMRI research and are performed frequently enough in the Winawer lab that we have developed a rough pipeline for their preparation and analysis. This document will step through this pipeline using a dataset collected at NYU, demonstrating the various processing done at each stage. This document is written primarily for use with NYU’s Center for Brain Imaging (CBI). Details about the preprocessing steps performed by CBI for NYU researchers (e.g., B0 fieldmap correction) can be found on the CBI intranet. This tutorial assumes that you have already downloaded your data from Tesla (non-NYU users will need to perform some inference in figuring out how to adapt the first few steps to their own local setups).
The stimulus and scan protocol used in the experiment that is analyzed in this document was nearly identical to that used in the Human Connectome Project (HPC); see their page on protocols for more information. In the particular experiment analyzed here, we performed 8 runs of retinotopic mapping, each of which lasted 192 TRs (1 second each). We additionally collected 1 T1-weighted image at a resolution of 0.8 mm\(^3\).
This tutorial begins in my
~/Downloads/rscan_20180205 directory (using Mac OS 10.13.2) to which I
have already downloaded all scan files from the CBI’s Tesla file-server. For
information on how to do this, wee the related page on the Winawer-lab
wiki. At the start of the
analyses detailed below, I observed the following directory contents:
> pwd /Users/nben/Downloads/rscan_20180205 > ls 01+AAHead_Scout_64ch-head-coil 16+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 02+AAHead_Scout_64ch-head-coil_MPR_sag 17+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 03+AAHead_Scout_64ch-head-coil_MPR_cor 18+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 04+AAHead_Scout_64ch-head-coil_MPR_tra 19+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 05+FMRI_DISTORTION_AP_2mm_66sl 20+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 06+FMRI_DISTORTION_PA_2mm_66sl 21+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 07+AAHead_Scout_64ch-head-coil 22+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 08+AAHead_Scout_64ch-head-coil_MPR_sag 23+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 09+AAHead_Scout_64ch-head-coil_MPR_cor 24+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 10+AAHead_Scout_64ch-head-coil_MPR_tra 25+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 11+FMRI_DISTORTION_AP_2mm_66sl 26+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 12+FMRI_DISTORTION_PA_2mm_66sl 27+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 13+T1_MPRAGE 28+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 14+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 29+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 15+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 99+PhoenixDocument
For those not familiar with the CBI’s naming convention, each of these items listed above is a directory, and each contains two files: a text file whose contents are a reproduction of the original DICOM header for the scan, and a NifTI file whose contents are the actual 3D or 4D scan measurements. The directories containing the scans are numbered by the order in which they were collected.
Note. Looking closely at the above directory contents, one might notice that we performed multiple scout and distortion scans prior to performing a T1 and the 8 retinotopy scans. This was, in fact, the result of a problem getting our stimulus to display propertly. In an ideal world this would never happen, but in reality, it is quite common to have a scan or two that must be discarded or ignored in a study. In this case, we will discard scans 1-6 (and note that we will not use the scout scans 7-10 in these analyses). The commands we use below will reflect the fact that we are ignoring these scans.
FreeSurfer is a software suite for processing anatomical MRI data; it includes a large number of useful algorithms that, among other things, identify the brain, strip the skull, identify white- and gray-matter voxels, perform various normalizations, tesselate the white and pial surfaces, and perform cortical surface alignment to an average anatomical atlas.
The first step in processing a retinotopy experiment is generally to process the T1-weighted
anatomical image. This step is highly complicated, but that complication is entirely performed by
FreeSurfer, using a single command,
> ls 13+T1_MPRAGE/ wl_subj046_ColorRetinotopy+13+T1_MPRAGE-header.txt wl_subj046_ColorRetinotopy+13+T1_MPRAGE.nii > recon-all -subjid wl_subj046 -i 13+T1_MPRAGE/wl_subj046_ColorRetinotopy+13+T1_MPRAGE.nii -all &> ./wl_subj046_recon-all.log
Estimated Runtime: 3-12 hours, depending on your computer, FreeSurfer version, the subject, and the scan. I have never made a habit of timing this command and usually just assume that it will run overnight. I have heard that subjects with more curvature in their cortices require more time to process (which makes some intuitive sense given that FreeSurfer performs cortical surface tesselation), but I have never confirmed this myself.
In the above command, we call FreeSurfer’s
recon-all command, which handles the importing of
subjects into FreeSurfer’s subject directory by performing an enormous suite of processing on the
T1-weighted anatomical image. Documentation on FreeSurfer is not always up-to-date and
comprehensive, but official information about the recon-all process is available
here. In my experience, more can usually be
learned from the recon-all dev-table,
though understanding the steps documented there may require requesting
--help documentation from
the various FreeSurfer commands. For example, to understand the normalization procedure used in the
first stage of recon-all, note that the command
mri_normalize is used; the best information about
this processing step will likely come from running
Although FreeSurfer’s results are generally fine, it’s not a bad idea to check for obvious problems
in the resulting data files. An easy way to do this is to use the program
freeview, which comes
with FreeSurfer. Freeview itself has variety of visualization functions that it performs, the extent
of which is beyond the scope of this document. This
gives a brief tutorial and links to other resources, however.
To check our subject, we will first look at their
brain.mgz file to make sure it looks reasonable,
then overlay on it the
ribbon.mgz file to make sure that there aren’t major errors in the
segmentation. To start this, we can run
freeview from the command-line:
# switch to the subject's new FreeSurfer directory > cd /Volumes/server/Freesurfer_subjects/wl_subj046 > ls mri/ T1.mgz orig aparc+aseg.mgz orig.mgz aparc.a2009s+aseg.mgz orig_nu.log aseg.auto.mgz orig_nu.mgz aseg.auto_noCCseg.label_intensities.txt rawavg.mgz aseg.auto_noCCseg.mgz rh.ribbon.mgz aseg.mgz ribbon.mgz brain.finalsurfs.mgz segment.dat nu.mgztalairach.label_intensities.txt talairach.log brainmask.auto.mgz talairach_with_skull.log brainmask.mgz talairach_with_skull_2.log ctrl_pts.mgz transforms filled.mgz wm.asegedit.mgz lh.ribbon.mgz wm.mgz mri_nu_correct.mni.log wm.seg.mgz mri_nu_correct.mni.log.bak wmparc.mgz norm.mgz brain.mgz nu_noneck.mgz # Start Freeview; load the brain.mgz file > freeview -v mri/brain.mgz
The above command should open a window that looks like this:
By scrolling through this subject’s data, it should be easy to determine if FreeSurfer correctly
identified the location of the subject’s brain; in this case, it did a pretty good job. To check,
the segmentation, however, we must load the ribbon file. The ribbon contains the same voxels labeled
as either LH white-matter, LH gray-matter, RH white-matter, RH gray-matter, or none. If we click on
the “Load Volumes” button in the upper left corner of the FreeView window (the head with the green
plus icon), we can add a volume to the display. Click this button then select the same subject’s
ribbon.mgz file, also in the
mri/ directory. At the bottom of the import-file window is a list
of options for the display/color for the volume; select “Lookup Table” for this. Once you’ve done
this, you should have a FreeView window that looks something like this:
To make the ribbon more transparent, you can adjust the “Opacity” slider in the menu-bar on the
left. Additionally, you can check or uncheck the “ribbon” volume near the upper-left corner to
toggle display of the ribbon file altogether. It’s generally a good idea to look through the
subject’s occipital cortex to make sure that the segmentation assigned by the ribbon looks
reasonable with respect to the voxels in the
brain.mgz; segmentation errors in the ribbon can
result in significant errors during the cortical surface generation. Although repairing errors is
beyond the scope of this tutorial, it involves hand-editing the wm.mgz file and restarting
recon-all at an intermediate stage. The wm.mgz file can be hand-edited using
ITK-Snap, and this
page gives an
introduction to processing such repairs.
Finally, we can check that the cortical surfaces were constructed reasonably. To do this, we click
Surfaces tab in the upper-left corner of FreeView, then on the “Load Surface” icon (the
brain with the green plus in the upper left). This will bring up a surface-loading dialog-box in
which you can navigate to your subjects
surf/ directory (in their FreeSurfer subject directory)
and select the
lh.white surface file. This will yield a display that looks something like this:
Again, repair of these surfaces is beyond the scope of this document; though, in the author’s
experience, most cortical surface errors are due to segmentation errors and should be corrected
prior to surface generation. Other surface files that should be checked for errors are
The next step, after processing the anatomical scan, is to preprocess the EPI scans. Most of this step is performed by the prisma_preprocess.py script; though the final step is performed by the to_freesurfer.py script.
prisma_preproc.py script does the bulk of the preprocessing work for EPIs. The script, by
Serra Favila, encapsulates a number of complicated steps including the unwarping/undistorting of the
EPI images, motion correction, and calculation of the alignment to the subject’s
FreeSurfer-processed anatomy data. This script can be found at the Winawer lab github
page in the MRI_tools
repository. Detailed information about running
prisma_preproc.py can be found on this Winawer lab
wiki page that documents it.
Before running the
prisma_preproc.py script, it is worthwhile to setup a VistaSoft session
directory; for retinotopy, we typically keep our retinotopy sessions in the
directory on our lab webserver. This bash session demonstrates this setup:
# Make the subject directory... > mkdir /Volumes/server/Projects/Retinotopy/wl_subj046 > cd /Volumes/server/Projects/Retinotopy/wl_subj046 # ...and the session directory... > mkdir 20180205_ColorRetinotopy > cd 20180205_ColorRetinotopy # Make the various directories > mkdir Raw > mkdir Preproc > mkdir Stimuli > mkdir Outputs > mkdir Code # Copy over the prisma_preproc.py and to_freesurfer.py scripts; # the ~/Code/MRI_tools directory contains the MRI_tools repo. > cp ~/Code/MRI_tools/preprocessing/*.py Code/ # We will also want the retinotopy code later. > cp ~/Code/MRI_tools/retinotopy/*.py Code/ # Also copy over stimulus movie and parameter files; these # specific files are for the experiment we ran: > cp /Volumes/server/Projects/Retinotopy/ColorRetinotopyShared/stimuli/scan_*.mat Stimuli # Transfer over the raw data > mv ~/Downloads/rscan_20180205/* Raw/
After executing the above, we are ready to run the preprocessing script:
> cd Raw > ls 01+AAHead_Scout_64ch-head-coil 16+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 02+AAHead_Scout_64ch-head-coil_MPR_sag 17+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 03+AAHead_Scout_64ch-head-coil_MPR_cor 18+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 04+AAHead_Scout_64ch-head-coil_MPR_tra 19+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 05+FMRI_DISTORTION_AP_2mm_66sl 20+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 06+FMRI_DISTORTION_PA_2mm_66sl 21+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 07+AAHead_Scout_64ch-head-coil 22+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 08+AAHead_Scout_64ch-head-coil_MPR_sag 23+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 09+AAHead_Scout_64ch-head-coil_MPR_cor 24+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 10+AAHead_Scout_64ch-head-coil_MPR_tra 25+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 11+FMRI_DISTORTION_AP_2mm_66sl 26+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 12+FMRI_DISTORTION_PA_2mm_66sl 27+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 13+T1_MPRAGE 28+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 14+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0_SBRef 29+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 15+cmrr_mbepi_s6_2mm_66sl_PA_TR1.0 99+PhoenixDocument # Actually run the script... > python ../Code/prisma_preproc.py -subject wl_subj046 \ -datadir "$PWD" \ -outdir ../Preproc \ -epis 15 17 19 21 23 25 27 29 \ -sbref 14 \ -distortPE 12 \ -distortrevPE 11 \ &> ../Preproc/prisma_preproc.log
Estimated Runtime: 1-8 hours; this step will require several hours to run, depending on the number of EPIs, their length, their resolution, as well as your computing power and memory limits.
This last command will produce a significant amount of output, all redirected to the
prisma_preproc.log file in the session’s
Preproc directory. To understand the choice of
paramters, see the wiki page
documenting the script or the script’s help-text (
python prisma_preproc.ph -h).
When this script has finished running, we can look at the results in the
Preproc directory. Errors
that arise during the preprocessing script are documented on the wiki
page; to check for errors, use
less prisma_preproc.log to review the output from the script.
> cd ../Preproc > ls distort2anat_tkreg.dat timeseries_corrected_run03.nii.gz distortion_merged_corrected.nii.gz timeseries_corrected_run04.nii.gz distortion_merged_corrected_mean.nii.gz timeseries_corrected_run05.nii.gz prisma_preproc.log timeseries_corrected_run06.nii.gz sbref_reg_corrected.nii.gz timeseries_corrected_run07.nii.gz session.json timeseries_corrected_run08.nii.gz timeseries_corrected_run01.nii.gz workflow timeseries_corrected_run02.nii.gz > tail prisma_preproc.log sub: /Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Preproc/_merge_epis4/timeseries_corrected.nii.gz -> /Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Preproc/timeseries_corrected_run05.nii.gz 180207-10:52:35,118 interface INFO: sub: /Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Preproc/_merge_epis5/timeseries_corrected.nii.gz -> /Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Preproc/timeseries_corrected_run06.nii.gz 180207-10:52:42,585 interface INFO: sub: /Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Preproc/_merge_epis6/timeseries_corrected.nii.gz -> /Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Preproc/timeseries_corrected_run07.nii.gz 180207-10:52:49,510 interface INFO: sub: /Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Preproc/_merge_epis7/timeseries_corrected.nii.gz -> /Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Preproc/timeseries_corrected_run08.nii.gz 180207-10:52:58,230 workflow INFO: [Job finished] jobname: outfiles jobid: 11 180207-10:52:58,237 workflow INFO: Currently running 0 tasks, and 0 jobs ready. Free memory (GB): 43.20/43.20, Free processors: 8/8
As you can see, the end of the
prisma_preproc.log does not contain any error messages; this is a
good sign and indicates that the script probably ran fine. In addition to the log file, we can see
that there are several
timeseries_corrected_run*.nii.gz files containing the unwarped
motion-corrected time-series data for each EPI. Most of the other files can be ignored with the
exception of the
distort2anat_tkreg.dat file, which contains the affine transform that
FreeSurfer can use to align the corrected timeseries data to the subject’s FreeSurfer anatomical
data. Warning: this file is not what it seems–specifically, it is not a simple
transformation from the affine-matrix stored in the time-series files to the subject’s FreeSurfer
anatomy (see next section). This alignment can be performed by the
to_freesurfer.py can be found in the same directory and location as the
script. The purpose of the script is to apply the
distort2anat_tkreg.dat alignment file to the
timeseries files in the
Preproc directory after running the
prisma_preproc.py script and to
create resampled cortical-surface timeseries files from the timeseries volume files; these can be
used to analyze the time-sereis data on the cortical surface directly.
Note: If you are using VistaSoft to analyze your pRF data (as done below in this document), then you do not need to perform this step. The Matlab scripts used to process the timeseries data will perform the realignment automatically.
The first step in the
to_freesurfer.py script is to apply the
transformation. Note that this file contains a
tkreg transform, which is strikingly unintuitive if
you are not used to FreeSurfer’s way of thinking about MRI volumes. Typically, we assume that an
alignment transformation that aligns image
A to image
B would store that alignment as a
transformation from image
A’s orientation to image
B’s orientation: i.e., from affine matrix
would tell us how to change
A’s affine transorm to align it to
B. In a
tkreg alignment matrix,
however, the transformation that is stored aligns
A’s “tkreg” matrix to B.
In FreeSurfer, every MR image file has, in addition to any affine transformation that it stores, a
vox2ras-tkr (or “tkreg”) matrix. This matrix can be deduced from other header information, and can
be printed using the
> mri_info --vox2ras-tkr ./timeseries_corrected_run01.nii.gz -2.00000 0.00000 0.00000 104.00000 0.00000 0.00000 2.00000 -66.00005 0.00000 -2.00000 0.00000 104.00000 0.00000 0.00000 0.00000 1.00000
This matrix is not the same as the image’s
vox2ras or affine matrices; rather it is deduced
from the voxel size, image dimensions, and slice ordering. Notably, it does not change when the
affine transformation encoded in the file changes. Because this matrix is somewhat unintuitive, I
recommend applying the
distort2anat_tkreg.dat transformation using the
# to_freesurfer.py includes various options: > python ../Code/to_freesurfer.py --help usage: to_freesurfer.py [-h] [-t TAG] [-s] [-o OUTDIR] [-m METHOD] [-l LAYER] [-d SDIR] [-v] registration_file EPI [EPI ...] positional arguments: registration_file The distort2anat_tkreg.dat or similar file: the registration file, in FreeSurfer\'s tkreg format, to apply to the EPIs. EPI The EPI files to be converted to anatomical orientation optional arguments: -h, --help show this help message and exit -t TAG, --tag TAG A tag to append to the output filenames; if given as - then overwrites original files. By default, this is "_anat". -s, --surf If provided, instructs the script to also produce files of the time-series resampled on the cortical surface. -o OUTDIR, --out OUTDIR The output directory to which the files should be written; by default this is the current directory (.); note that if this directory also contains the EPI files and there is no tag given, then the EPIs will be overwritten. -m METHOD, --method METHOD The method to use for volume-to-surface interpolation; this may be nearest or linear; the default is linear. -l LAYER, --layer LAYER Specifies the cortical layer to user in interpolation from volume to surface. By default, uses midgray. May be set to a value between 0 (white) and 1 (pial) to specify an intermediate surface or may be simply white, pial, or midgray. -d SDIR, --subjects-dir SDIR Specifies the subjects directory to use; by default uses the environment variable SUBJECTS_DIR. -v, --verbose Print verbose output # (Optional) To preserve the original transformations, create a new directory > mkdir ../Preproc_FreeSurfer > cp ./distort2anat_tkreg.dat ./timeseries*.nii.gz ../Preproc_FreeSurfer > cd ../Preproc_FreeSurfer # It is typically run like this; this will only correct the timeseries headers > python ../Code/to_freesurfer.py -v ./distort2anat_tkreg.dat timeseries*.nii.gz Processing EPI ./timeseries_corrected_run01.nii.gz... - Correcting volume orientation... Processing EPI ./timeseries_corrected_run02.nii.gz... - Correcting volume orientation... Processing EPI ./timeseries_corrected_run03.nii.gz... - Correcting volume orientation... Processing EPI ./timeseries_corrected_run04.nii.gz... - Correcting volume orientation... Processing EPI ./timeseries_corrected_run05.nii.gz... - Correcting volume orientation... Processing EPI ./timeseries_corrected_run06.nii.gz... - Correcting volume orientation... Processing EPI ./timeseries_corrected_run07.nii.gz... - Correcting volume orientation... Processing EPI ./timeseries_corrected_run08.nii.gz... - Correcting volume orientation... # To also produce surface data-files > python ../Code/to_freesurfer.py -v -s ./distort2anat_tkreg.dat timeseries*.nii.gz Processing EPI ./timeseries_corrected_run01.nii.gz... - Correcting volume orientation... - Projecting to surface... Processing EPI ./timeseries_corrected_run02.nii.gz... - Correcting volume orientation... - Projecting to surface... Processing EPI ./timeseries_corrected_run03.nii.gz... - Correcting volume orientation... - Projecting to surface... Processing EPI ./timeseries_corrected_run04.nii.gz... - Correcting volume orientation... - Projecting to surface... Processing EPI ./timeseries_corrected_run05.nii.gz... - Correcting volume orientation... - Projecting to surface... Processing EPI ./timeseries_corrected_run06.nii.gz... - Correcting volume orientation... - Projecting to surface... Processing EPI ./timeseries_corrected_run07.nii.gz... - Correcting volume orientation... - Projecting to surface... Processing EPI ./timeseries_corrected_run08.nii.gz... - Correcting volume orientation... - Projecting to surface... > ls distort2anat_tkreg.dat rh.timeseries_corrected_run06.mgz lh.timeseries_corrected_run01.mgz rh.timeseries_corrected_run07.mgz lh.timeseries_corrected_run02.mgz rh.timeseries_corrected_run08.mgz lh.timeseries_corrected_run03.mgz timeseries_corrected_run01.nii.gz lh.timeseries_corrected_run04.mgz timeseries_corrected_run02.nii.gz lh.timeseries_corrected_run05.mgz timeseries_corrected_run03.nii.gz lh.timeseries_corrected_run06.mgz timeseries_corrected_run04.nii.gz lh.timeseries_corrected_run07.mgz timeseries_corrected_run05.nii.gz lh.timeseries_corrected_run08.mgz timeseries_corrected_run06.nii.gz rh.timeseries_corrected_run01.mgz timeseries_corrected_run07.nii.gz rh.timeseries_corrected_run02.mgz timeseries_corrected_run08.nii.gz rh.timeseries_corrected_run03.mgz rh.timeseries_corrected_run04.mgz rh.timeseries_corrected_run05.mgz
Estimated Runtime: a couple minutes per EPI.
Solving the pRF models based on the timeseries data is performed using VistaSoft in Matlab. The three scripts documented below are designed to be run on the Winawer-lab server using the directory structures shown in this document. For deviations from this narrow use-case, you will need to edit the scripts themselves. Mostly, these edits can probably be small (e.g., if you are using a different set of directory names), but for different pipelines, these scripts should be seen as a rough guide and not a solution. In these cases, I would recomment the VistaSoft Ernie tutorials and related VistaSoft tutorials on MRI data analysis.
All of these scripts can be found in the MRI_tools
repository in the
retinotopy directory. They require that the
ToolboxToolbox for Matlab be installed; this will
manage other dependencies, including VistaSoft. Recall that we already copied all of these scripts
into our session’s
Code directory when we setup the session directory. All three of these scripts
will deduce the subject, session, and directory structure when run as long as you are using a
similar structure as this set of demos. You can also edit the top parts of the scripts
init_vista script initializes the subject’s VistaSoft session; this primarily involves
initializing certain data structures, finding the preprocessed timeseries files, and applying the
appropriate alignments. If the subject does not yet have a VistaSoft anatomy directory in the
/Volumes/server/Projects/Anatomy directory (where we store VistaSoft anatomical sessions in the
Winawer lab), it will create one of these and initialize it from the subject’s FreeSurfer
init_vista, simply invoke it from inside your session’s
>> cd('/Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Code'); >> init_vista %% (This command generates significant output that is not reproduced here.)
Estimated Runtime: fast, usually less than a minute; could stretch up to 10 minutes if it has to initialize the Anatomy session and your computer is fairly slow.
solve_pRFs script calculates the actual solutions to pRF models for all of the subject’s
gray-matter voxels using the timeseries data. If you don’t edit the script, this creates a single
dataset named ‘Full’, which contains the average of all the timeseries data, and solves only this
dataset. If you wish to analyze separate subsets of the data, you must edit this script. That said,
adding new datasets is trivial: simply search your local copy of the script for
read the comment containing it.
solve_pRFs, simply invoke it from inside your session’s
>> cd('/Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Code'); >> solve_pRFs %% (This command generates significant output that is not reproduced here.)
export_niftis script writes the pRF model solutions out to a set of nifti files in your
Outputs directory. As with
solve_pRFs, the default behavior is to export a single
‘Full’ dataset; to change this, you will need to edit the script. As with
solve_pRFs, this is made
easy by instructions in a comment tagged with
export_niftis, simply invoke it from inside your session’s
>> cd('/Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Code'); >> export_niftis %% (This command generates significant output that is not reproduced here.)
The files output by
export_niftis contain the following data (here, as if the output prefix were
‘full’, which is also the default behavior):
full-ycrds.nii.gzcontain the coordinates of the pRF centers in degrees of the visual field.
full-vexpl.nii.gzcontains the proportion of variance explained by the model (a fraction between 0 and 1).
full-sigma.nii.gzcontains the effective pRF radius in degrees of the visual field.
The final script in the
retinotopy directory of the
MRI_tools repository is the
script. This script reads in the files exported by the
export_niftis.m script and creates volume
files for polar angle and eccentricity as well as surface files for all data types (polar angle,
eccentricity, x, y, variance explained, sigma). The script requires merely the FreeSurfer subject id
and the path to the
Outputs directory. It automatically operates on all datasets found there.
> pwd /Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy > ls Outputs all-sigma.nii.gz all-vexpl.nii.gz all-xcrds.nii.gz all-ycrds.nii.gz > python ~/Code/MRI_tools/retinotopy/postproc_pRFs.py wl_subj046 Outputs -v Processing Dataset all... - Importing parameters... - Creating polar angle/eccentricity images... - Projecting to surface... > ls Outputs all-angle.nii.gz all-xcrds.nii.gz lh.all-sigma.mgz rh.all-angle.mgz rh.all-xcrds.mgz all-eccen.nii.gz all-ycrds.nii.gz lh.all-vexpl.mgz rh.all-eccen.mgz rh.all-ycrds.mgz all-sigma.nii.gz lh.all-angle.mgz lh.all-xcrds.mgz rh.all-sigma.mgz all-vexpl.nii.gz lh.all-eccen.mgz lh.all-ycrds.mgz rh.all-vexpl.mgz
A simple way to check that the above steps ran correctly is to visualize the solved pRF parameters projected to the cortical surface; this essentially represents the culmination of all of the steps documented so far in these tutorials. FreeView is capable of performing this kind of visualization, as well as various other tools and libraries. For simple checks, FreeView is sufficient; for more complex visualization, I would suggest a combination of neuropythy and pycortex, or Neurotica if you use Mathematica.
To start FreeView, we want to load the subject’s inflated hemisphere, which makes visualization easiest; alternately, loading the pial or white hemispheres can be useful for orienting oneself relative to anatomical landmarks. Here we will look at only one hemisphere, but examining the other hemisphere should be straightforward given the commands for one.
# Open FreeView with the subject's LH inflated surface > freeview -f "$SUBJECTS_DIR"/wl_subj046/surf/lh.inflated
This should open FreeView with the inflated left hemisphere shown in the lower right display
window. To add an overlay, navigate to the
Surfaces tab in the upper left of the FreeView
window. When you click on
Surfaces, the left menu will change; among the surface controls revealed
is a dropdown menu labeled ‘Overlay’. One of the options in this list is ‘Load generic overlay…’.
Selecting this item will allow you to navigate to your subject’s
Outputs directory to load the
lh.all-angle.mgz file, containing the polar-angle measurements for the subject’s left
hemisphere. Once the overlay is loaded, it will likely need to be configured via the
Overlay button. These controls are not perfectly intuitive, but they are not difficult to figure
out by trial and error. When properly configured, you should see something like this:
Anatomical atlases are a common and useful way to examine how your subject compares with average
retinotopic mapping data. The atlases described here use cortical surface normalization (i.e., via
fsaverage_sym subjects) in order to describe the average expected
location of various ROIs and/or their average expected organizational properties. Here we will look
at two atlases:
- The “Benson-2014” atlas: Benson NC, Butt OH, Brainard DH, Aguirre GK (2014) Correction of Distortion in Flattened Representations of the Cortical Surface Allows Prediction of V1-V3 Functional Organization from Anatomy. PLOS Comput Biol 10(3):e1003538. This atlas describes the location and retinotopic organization of V1, V2, and V3, as well as several higher areas (though these areas are not considered authoritative by the authors as they have not been tested). Note that the Benson-2014 atlas predicts retinotopic maps from 0-90 degrees of eccentricity, though only the inner 20 degrees of eccentricity in V1, V2, and V3 has been validated. This atlas provides visual area label, polar angle, eccentricity, and pRF radius (sigma) predictions.
- The “Wang-2015” atlas: Wang L, Mruczek RE, Arcaro MJ, Kastner S (2015) Probabilistic Maps of Visual Topography in Human Cortex. Cereb Cortex 25(10):3911-31. This atlas describes the locations/boundaries of 25 visual areas with emphasis on dorsal visual areas. This atlas does not describe retinotopy but only average expected ROI labels.
The easiest way to apply the Wang-2015 atlas is to use the
occipital_atlas docker, information
about which can be found here. See this link for
documentation on using the Docker.
The easiest way to apply the Benson-2014 atlas is to, again, use
neuropythy. The library includes a builtin command,
benson14_retinotopy that places the atlas files in your subject’s
directories–surface files in the
surf/ directory and volumetric images in the
directory. It can be executed like so:
> python -m neuropythy benson14_retinotopy wl_subj046 -v Processing subject wl_subj046: - Interpolating template... - Exporting surfaces: - Exporting LH prediction file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/surf/lh.benson14_varea - Exporting LH prediction file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/surf/lh.benson14_eccen - Exporting LH prediction file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/surf/lh.benson14_angle - Exporting LH prediction file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/surf/lh.benson14_sigma - Exporting RH prediction file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/surf/rh.benson14_varea - Exporting RH prediction file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/surf/rh.benson14_eccen - Exporting RH prediction file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/surf/rh.benson14_angle - Exporting RH prediction file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/surf/rh.benson14_sigma - Exporting Volumes: - Preparing volume file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/mri/benson14_varea.mgz - Exporting volume file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/mri/benson14_varea.mgz - Preparing volume file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/mri/benson14_eccen.mgz - Exporting volume file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/mri/benson14_eccen.mgz - Preparing volume file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/mri/benson14_angle.mgz - Exporting volume file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/mri/benson14_angle.mgz - Preparing volume file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/mri/benson14_sigma.mgz - Exporting volume file: /Users/nben/Documents/WinawerLab/Freesurfer_subjects/wl_subj046/mri/benson14_sigma.mgz Subject wl_subj046 finished! > ls "$SUBJECTS_DIR"/wl_subj046/*/*benson14* /Volumes/server/Freesurfer_subjects/wl_subj046/mri/benson14_angle.mgz /Volumes/server/Freesurfer_subjects/wl_subj046/mri/benson14_eccen.mgz /Volumes/server/Freesurfer_subjects/wl_subj046/mri/benson14_sigma.mgz /Volumes/server/Freesurfer_subjects/wl_subj046/mri/benson14_varea.mgz /Volumes/server/Freesurfer_subjects/wl_subj046/surf/lh.benson14_angle /Volumes/server/Freesurfer_subjects/wl_subj046/surf/lh.benson14_eccen /Volumes/server/Freesurfer_subjects/wl_subj046/surf/lh.benson14_sigma /Volumes/server/Freesurfer_subjects/wl_subj046/surf/lh.benson14_varea /Volumes/server/Freesurfer_subjects/wl_subj046/surf/rh.benson14_angle /Volumes/server/Freesurfer_subjects/wl_subj046/surf/rh.benson14_eccen /Volumes/server/Freesurfer_subjects/wl_subj046/surf/rh.benson14_sigma /Volumes/server/Freesurfer_subjects/wl_subj046/surf/rh.benson14_varea
Estimated Runtime: 20 minutes or less; this step is not terribly intensive, but template interpolation can take awhile on a slower computer.
Note that the surface-file outputs of this command are in FreeSurfer’s ‘curv’ format (
if you are using
nibabel). If you prefer
mgz format, you can use the
python -m neuropythy benson14_retinotopy --help for more information. Similarly, the volume
files may be written as
.nii.gz nifti files instead of
.mgz files via the argument
Performing Bayesian inference on the retinotopic mapping data in order to generate an inferred map
prediction is as simple as a single command (albeit one with several arguments) executed by the
neuropythy library. The following code-block
demonstrates this. After running the command, files named similar to
describe the predicted polar angle, eccentricity, visual area (
varea) and pRF radius (
> pwd /Volumes/server/Projects/Retinotopy/wl_subj046/20180205_ColorRetinotopy/Outputs > ls all-angle.nii.gz all-xcrds.nii.gz lh.all-sigma.mgz rh.all-angle.mgz rh.all-xcrds.mgz all-eccen.nii.gz all-ycrds.nii.gz lh.all-vexpl.mgz rh.all-eccen.mgz rh.all-ycrds.mgz all-sigma.nii.gz lh.all-angle.mgz lh.all-xcrds.mgz rh.all-sigma.mgz all-vexpl.nii.gz lh.all-eccen.mgz lh.all-ycrds.mgz rh.all-vexpl.mgz > python -m neuropythy register_retinotopy wl_subj046 \ --verbose \ --surf-outdir=. --surf-format="mgz" \ --no-volume-export \ --lh-angle=lh.all-angle.mgz \ --lh-eccen=lh.all-eccen.mgz \ --lh-weight=lh.all-vexpl.mgz \ --lh-radius=lh.all-sigma.mgz \ --rh-angle=rh.all-angle.mgz \ --rh-eccen=rh.all-eccen.mgz \ --rh-weight=rh.all-vexpl.mgz \ --rh-radius=rh.all-sigma.mgz Processing subject: wl_subj046 Preparing RH Registration... Preparing LH Registration... Exporting files... Extracting RH predicted mesh... Extracting LH predicted mesh... # This command produces some new files: > ls all-angle.nii.gz lh.all-xcrds.mgz rh.all-vexpl.mgz all-eccen.nii.gz lh.all-ycrds.mgz rh.all-xcrds.mgz all-sigma.nii.gz lh.inferred_angle.mgz rh.all-ycrds.mgz all-vexpl.nii.gz lh.inferred_eccen.mgz rh.inferred_angle.mgz all-xcrds.nii.gz lh.inferred_sigma.mgz rh.inferred_eccen.mgz all-ycrds.nii.gz lh.inferred_varea.mgz rh.inferred_sigma.mgz lh.all-angle.mgz lh.retinotopy.sphere.reg rh.inferred_varea.mgz lh.all-eccen.mgz rh.all-angle.mgz rh.retinotopy.sphere.reg lh.all-sigma.mgz rh.all-eccen.mgz lh.all-vexpl.mgz rh.all-sigma.mgz
Estimated Runtime: Variable, but if you stick to the default number of steps and have a reasonably powerful computer, this should take less than an hour per hemisphere. If you have only one core, for example, this may take a few hours. With a powerful machine, this command can take as little as 15 minutes per hemisphere.