pRF methods

Dumoulin and Wandell (2008) developed a model-based method that uses the BOLD time course to estimate models of population receptive fields (pRF). Here we describe the mrVista implementation of this method.

• There is a pRF Tutorial page which walks through a particular sample pRF data set. 'this is a zip file from white... needs to be changed to RDT and probably new data: 'TODO for Jon and Nathan'

• There is also a Tutorial on older "traveling wave" analyses of retinotopy data.

The pRF method works well in V1/V2. We have discovered that the assumption of spatial summation fails in much of extrastriate cortex. We are developing an approach to managing based on a | Compressive spatial summation model (Kay et al 2013). Visit the CSS model page for more information.

To see how to run the pRF stimuli at the CNI scanner, go to the pRF_stimuli page.

Table of Contents Introduction Preprocessing Analysis Overview Set stimulus parameters HRF Detrending Other parameters Model estimation View results Retinotopy model: Main functions Using the command line to fit pRF models Fitting procedures Creating a new stimulus Creating a new 'make stimulus' file Using a stored image matrix Jittering images How to build a new pRF model Algorithms PRF model simulation Programming notes Script Notes

Introduction

The pRF analysis estimates the parameters of an explicit model of the population of receptive fields within an MR voxel. The simplest pRF model is a two-dimensional circular Gaussian. In addition, we have implemented a difference of circular Gaussian an oval Gaussian model, sum of positive Gaussians, sum of unsigned Gaussians and others. Some of these models are accessible for the graphical interface, but others require evaluation from the command line. We anticipate the list of models will expand over the years.

In the first part of this tutorial, we describe the procedure for fitting the circular, symmetric pRF. For this case the field map is defined by the estimated center positions of the Gaussians, (x,y), and the population receptive field size is determined by the Gaussian standard deviation (sigma). Other pRF summaries, such as the proportion of the receptive field in the contra- and ipsilateral visual field (laterality) can be derived and visualized from the estimated parameters.

Later in the tutorial we explain how to run the code from the command line and write code for more general models. We describe the procedure using the data set which can be downloaded from Not Yet Implemented. You can follow through this page using the data set.

Preprocessing

The general preprocessing steps (motion correction, slice-timing correction, registration to new coordinate frames) are the same as any analysis. The code runs only in the Gray view, so that you must have a segmentation of the subject's brain (see Gray-White Segmentation Pipeline).

this code can now also be run in the inplane, though I am not sure this functionality has been pushed to the master

The current implementation is built and tested in the "Gray" view. Transform your data into this "Gray" view.

To transfer your data to the "Gray" view open both "Inplane" and "Gray" window. You can transform the raw data by selecting in the "Gray" window: Xform | Inplane -> Volume | tSeries (all scans) | trilinear interpolation.

An advantage of using the "Gray" view is that you can easily add other fMRI runs from different sessions or days. The fMRI data from different session can be ported into the current fMRI session by selecting in the "Gray" window: Xform | Other Session -> This Session | Import tSeries only.

The pRF analysis analyses all scans within the selected data-type (e.g., Original, Averages). If you want analyze only on a subset of the scans, create a new data-type with only those scans. One way to do this is using the averageTSeries command multiple times:

VOLUME{1} = averageTSeries(VOLUME{1},[scansToaverage],'myNewDataTypeName');

And choosing just one scan to average each time, while keeping myNewDatatTypeName the same.

There should be a simple routine for selecting scans and creating a new type. Any volunteers?

Analysis Overview

In the test data set that we have provided, open up the Gray view and select the Data Type 'Averages'.

You can perform the One Gaussian retinotopy analysis from the user interface you would

• Set the parameters (Analysis | Retinotopic Model | Set Parameters)
• Load the parameters (Analysis | Retinotopic Model | Load Parameters)
• View stimulus (Plot | Retinotopic Model | View Stimulus)
• Run analysis (Analysis | Retinotopic Model | Run ) or in Unix you can call (Analysis | Retinotopic Model | Run in background )

We explain these steps below

Set stimulus parameters

The pRF method uses the experimental stimulus as part of the estimation process. Certain classes of stimuli (e.g., rings, wedges and bars) are used in many experiments. These stimuli can be recreated from the analysis program if you specify the parameters.

You can specify these stimuli from: Analysis | Retinotopic Model | Set Stimulus Parameters. You will get a window like the one shown below. In this example data set, the parameters are already set. Just view them and click 'OK'.

Set stimulus parameters

For each Scan number specify a stimulus type. Three predefined stimuli (wedge, ring or bars) are available. For each of these types you are asked to specify the stimulus size, which is the radius of the most peripheral point of the contrast signal (in deg from the fovea). You will be asked for other parameters that may be specific to that stimulus.

If you did not use one of the predefined stimuli (wedge, ring or bars), you must define the stimulus type from a stored image matrix and parameters file. To do this, choose 'StimFromScan' as the "stimulus type". You then need to choose an image file and a parameters file from the drop-down fields at the bottom of the GUI. See [#Creating].

If you used bars, the width should be set using the following formula: 360/(params.radius/params.ringDeg). Often this will be 90. The params are saved out (optionally) when running ret on the stimulus laptop. The formula above is a hack since the 8 bars code is based on originally using rings and wedges.

The stimulus window also asks you about some important non-stimulus parameters. These are essential because they depend on the stimulus and subject -

Hemodynamic response function (HRF) model (default = one gamma Boynton-style)

Maximum DCT frequency to be removed for detrending (default = 1)

To verify that you set the same stimulus as the one used in the fMRI section use: Plots | Retinotopic Model | View stimulus aperture.

To make a standard type of stimulus, you can create a function that makes the stimulus. Named the function make[stimName].m. For an example see makeWedge.m or make8Bars.m.

HRF

It is possible to define all of the stimulus-related (and other) parameters from the command line. Two particular important terms are the HRF and detrending parameters.

If you have specific knowledge about the subject's HRF, you can adjust either the function or function parameters.

Detrending

The DCT frequency determines how the detrending is performed. The basic analysis removes low frequency DCT terms. The number you would like to remove depends on the stimulus parameters.

Other parameters

There are many additional parameters. An extensive set of comments and detailed description is provided in the rmDefineParameters.m file. Important parameters and examples are also summarized in rmMain.m. All the parameters have defaults in those functions. These can be changed using command line arguments.

Model estimation

Only the One Gaussian model can be run directly from the window interface. To run this model select

 '''Analysis | Retinotopic Model | Coarse to fine minimization | Run Now '''. 

This executes rmMain. This can take a long time on the whole brain.

You can run the analysis on only a region of interest. If ROIs are loaded the analysis will be performed for the current ROI only; when no ROIs exist the analysis will run on the entire dataset.

(On unix systems, you can run the analysis in the background by using Analysis | Retinotopic Model | Run in background (unix) . This function will run the code in a separate shell script - the script will not die if you quit your current matlab session or even when you logout of the machine (use kill -9 pid and typically you can find the pid using ps or top ).

You can also run the model from the command line. This is useful for scripting purposes and allows you to changing various default parameters. For example:

VOLUME{1} = rmMain(VOLUME{1},[],'coarse to fine','model',{'one gaussian','difference of gaussians'});

runs the coarse-to-fine analysis on the default model ('one gaussian') as well as a a more complex model ('difference of gaussians'). For more options see rmMain.m help rmMain and rmDefineParameters.m.

View results

To view the results load a retinotopic model file using: File | Retinotopic Model | Select and Load Model.

You can clear this model by selecting File | Retinotopic Model | Clear.

To get basic information about the model and on when and where it was estimated try: File | Retinotopic Model | Model information.

You can visualize the parameters using File | Retinotopic Model | Load Model Parameter. You will see an interface as shown at the right.

250 px

To load the parameters from the model you specify three fields

a) Which model (if you estimated more than one), b) The parameter (e.g., eccentricity, polar angle, pRF size, variance explained) c) Which data field to store the parameter (co , amp, ph or map).

Selection (c) is obscure - sorry. The mrVista data fields have historical names that were developed for the traveling wave. We are simply storing the pRF fields into those structures. At some point we will make a more sensible arrangement.

If you don't know what the latter statement means simply leave it at the default. The default automatically stores the parameter in the most appropriate field. The default assignments are:

• co <- Variance Explained or coherence
• ph <- Polar angle
• map <- Eccentricity, pRF size, anything else

To visualize the relevant data using the conventional pulldown information from View | [Coherence]

To plot the fits and the pRF for a given voxel in the current ROI use: Plot | Retinotopic Model | Receptive field and model fit | fitted timeseries.

To get the average values of a parameter in the current ROI, such as pRF size, use: Plot | Retinotopic Model | ROI stats. The average is weighted by the percent variance-explained at each voxel.

Retinotopy model: Main functions

The retinotopy modeling functions generally begin with rm, as in rmMain, or rmDefineParameters.

Other important routines are the two different search strategies:

• rmGridFit
• rmSearchFit

These are typically called in sequence. The grid fit is a coarse analysis. The search fit refines the grid fit estimates.

Using the command line to fit pRF models

Here is some example code showing how to fit the pRF models from the MATLAB command line. this assumes that you are in the session directory and have already opened a gray view (VOLUME{1} below).

     view = rmMain([view],[roiFileName],[wSearch],[varargin]);
     vw = VOLUME{1};
     roiFileName = 'leftMT';  % If empty, the entire volume is fit - takes a while
     searchType  = 'coarse to fine';
     outFileName = 'myFile'
     prfModels = {'one gaussian','difference of gaussians'};
     VOLUME{1} = rmMain(vw,roiFileName,searchType,'matFileName', outFileName,'model',prfModels);

You can adjust many of the parameters of the fit using the calling arguments. The main calling parameters you can adjust are shown by typing doc mrMain. A complete list is in the function rmDefineParameters.

There are many other parameters that are set by default. The defaults can be overwritten by a command line call (see Model estimation section).

Fitting procedures

Creating a new stimulus

There are two ways to create a new stimulus, either by creating a new 'make' file or by using an image matrix from your scan.

It is important to note that the stimulus we make for use with retinotopic model code is not identical to a stimulus used in a scan. We use simplified stimuli for solving models. For example, there is only one image per volume acquisition (TR), and the images usually show only the stimulus aperture (e.g., a wedge), and not the carrier, or internal contrast patterns (e.g., a checkerboard). This is why we need to use special make functions.

Creating a new 'make stimulus' file

A make* file should be named make[StimulusDescription].m and placed in the directory StimulusDefinitions. The new stimulus definition will automatically appear in the gui. Current examples include make8Bars, makeWedge, and makeRing. The make* file generates stimulus images. It takes params, a struct with stimulus parameters, as input and output, adding the field "images" to the struct. The stimulus parameters can be generated from the set parameters ui. There is one image per temporal frame, and the images are vectorized for ease of computation. The images are usually binary, with a 1 in pixel locations with non-zero contrast and 0s elsewhere. See existing make* files for further details, and rmMakeStimulus to see how the make* files are called.

Using a stored image matrix

Alternatively, you can make your stimulus starting with the raw images used in the scan. To do this, you need two files, an image file and a params file. The file names should contain the strings 'image' and 'params', respectively, and the files should be placed in a directory called 'Stimuli', inside the main project directory.

ExpTools Stimulus GUI

Both the image file and the params file can be made automatically if you use our expTools software for stimulus presentation. For example, if you run 'ret', our code for presenting retinotopic stimuli, you can check the box 'Save Stim Parameters', and enter a file name in the field 'Store image matrix' (see figure on right). This will result in two files, a parameters file and an image matrix file. The parameters file name is derived from the date and time you ran the program (e.g. 20080717T161217.mat)

If you do not use expTools to run your experiment in the scanner, or if you want to create arbitrary new stimuli, you can generate the image file and params files manually. The image file must have a variable named images, containing all the unique images shown in the scan. It is n * m * k, where n and m are the image dimensions and k is the number of images. The parameters file must have a variable named stimulus containing a field seq (stimulus.seq), and seqTiming (stimulus.seqTiming). Seq is a vector indexing the image matrix. seqTiming is a vector of the same length as seq indicating onset times of each image in seconds (zero-based, i.e. the first time point of the scan is zero). See makeStimFromScan for additional optional fields and futher details.

Once your image and params files are in the Stimuli directory, you can make your stimuli via the set parameters ui.

ExpTools Stimulus GUI

In the GUI, as shown on the right, select StimFromScan for Class. This will cause a number of stimulus fields to gray out; they will not be needed since the stimulus will be defined from the images themselves. The fields 'Image file', 'Params file', and 'Image filter' will be un-grayed, as these fields become necessary when defining your stimulus. The image filter is used to convert the raw images into stimulus predictors. The binary filter is most similar to the processing used in our standard make* functions, such as make8Bars, but other filters are possible.

Jittering images

If you have eye movement data you can alter your stimulus to compensate for eye position. To do this, put a file containing the string jitter inside the Stimuli directory. From the set parameters ui, you can then select a jitter file from the drop down menu. Eye movement jitter can be used with any kind of stimulus, whether made from a standard make* file or from your own image matrix.

How to build a new pRF model

Algorithms

The main reason why the pRF model runs in the "Gray" view is because a coarse-to-fine algorithm is implemented with two stages. In the first stage the data pRF parameter estimates are obtained from data that is smoothed across the cortical surface. In the second stage, these estimates are refined without any smoothing of the data. This approach allows for both a speed improvement as well as a more robust estimate of the pRF parameters. To be able to smooth across the cortical surface the data has to be in the "Gray" view.

PRF model simulation

From stored pRF models we create simulated data sets. These simulations have several uses. We can use this method to predict the response to novel stimuli. We can use the simulations to predict the consequences of eye movements, scotomas or other neurological defects.

This section reviews how use retinotopic model (RM) to calculate a simulated data set.

The key function for the simulation is rmCalculatePredictions. The simulation also uses rmEditStimParms, which calls a sequence of functions including rmLoadParameters and rmMakeStimulus. (These routines should probably be renamed, too).

The method of calculating a simulation from the GUI is

• Load an existing RM from the File | Retinotopy Model | Select and Load Model pull down
• From the Analysis | Retinotopic Model | Set Parameters pull down, set the properties of the stimulus for your scans
• Choose the Analysis | Retinotopic Model | Calculate pRF Predictions pull down; a new dataType will be created (default: Simulated).

You can view the simulated data by choosing this dataTYPE and interacting with the time series and other elements (ROIs) as if it were real data.

Programming notes

When you open a data set that has a retinotopy model calculated, you can load it from File | Retinotopy Model | Select and Load Model. This will attach a structure (rm) to the VOLUME{} cell array for that window. For example,

>> VOLUME{1}.rm ans = retinotopyModelFile: 'Y:\users\Netta\HH080213b\Gray\Averages\retModel-20080214-221919-sFit.mat' retinotopyModels: {[1x1 struct]} retinotopyParams: [1x1 struct] modelNum: 1

To simulate data, we might want to use the model and a different stimulus. For example, Jon does this when he estimates the model using retinotopy and then predicts the response to a visual brightness illusion.

THe next two steps then involve (a) creating and loading a new stimulus, and (b) combining the stimulus with the model to produce a simulated data set (prfSavePredictions).

To specify a new stimulus, you can use the function params = rmEditStimParams(getCurView);

To load a new stimulus you can either use the current method of generating stimuli in the GUI to create a new stimulus and put it in the retinotopy model structure. The pulldown is Analysis | Retinotopic Model | Set Parameters which brings up a window where you define your new stimulus. When you are done, click OK and (a) the stimulus will be generated, and (b) the parameters will be stored in the view structure. The stimulus itself will replace the current stimulus in the rm.retinotopyParams.stim, so be careful not to over-write the existing model (this would be hard).

Script

To be done

1. We need a script to do all of this

This sort of thing ...

curModel = viewGet(VOLUME{1},'rmCurrent') curModelParams = viewGet(VOLUME{1},'rmStimParams')

rmCalculatePredictions;

Notes

There is some confusing terminology

1. Repetitions: If you have 6 stimulus cycles and 4 blanks then the stimulus repeats twice. SOD accounts for this when solving the model. He averages the first and second half and fits that. This saves time. But the meaning of repetitions is confusing - cycles is the number of cycles that, say, the wedge goes around. Repetitions is the number of complete repetitions of the stimulus, including both the wedge/ring and the blanks. In the convenitonal 6 cycle 4 blank experiment the number of repetitions is 2. - Put this comment in the retinotopy model solution, too.