pRF Tutorial

Maybe this should be deprecated for the one Nathan drafted, based on this. pRF_Tutorial, I think (BW).

This page contains a tutorial on the population receptive field (pRF) method first described in [http://white.stanford.edu/~brian/papers/mri/2007-Dumoulin-NI.pdf Dumoulin and Wandell (2008)]. It walks you through the analysis of a sample data set available online.

• pRF_Model

== (1) Download the sample data set ==

The sample data set for this tutorial is downloadable at [http://white.stanford.edu/software/downloadsData/mrVistaPRFSession.zip '''http://white.stanford.edu/software/downloadsData/mrVistaPRFSession.zip'''].

This .zip file contains a compressed session directory, with data already initialized in mrVista format. The experiment for this session was a set of rotating-wedge and contracting-ring retinotopy scans. For each experimental run, the stimulus traversed the visual field (the wedge rotated around the screen, or the ring contracted from peripheral to foveal) 6 cycles. There was also an intermittent 12-second mean luminance or "blank" period, where the stimulus turned off, at a different frequency of 4 cycles/run.

The average wedge and ring runs have already been computed, and are stored in the data type "Averages". A segmented T1-weighted whole brain anatomy is also included in the 3DAnatomy/ subdirectory. This includes pre-built gray matter surface Mesh.

This session can be analyzed using either the Retinotopy_Tutorial or pRF_Model method. This tutorial describes how to use the pRF method.

You should place this directory somewhere in your file system to analyze. You should also have a running version of MATLAB (r2006a or newer), and the Software repository downloaded and present in your MATLAB path.

== (2) Start up mrVista ==

Open MATLAB and navigate to the directory with the example session.

Type mrVista . This will open an Inplane window. You should be familiar with navigating around this window for this tutorial.

mrVista will load the pre-computed Retinotopy_Tutorial file, computed as part of the traveling wave analyses. You should be viewing the phase map from this analysis. If you don't see this map, select the menu '''View''' | '''Phase Map'''.

{| |Image:pRF Tutorial Inplane with CorAnal.jpg |}

Although these analyses are not part of the pRF analysis, it is useful to look at the data, to ensure the data are sound. For the wedge and ring experimental design here, it is possible to run both the pRF analysis, as well as the earlier traveling wave analysis. Other designs, such as the "8 bar" stimulus in which bars sweep across visual space in different directions, can only be analyzed with the pRF method.

Now let's open up a gray 3-view window. Type mrVista 3 at the MATLAB prompt. A new Volume window will open.

Make sure the "Averages" data type is selected.

{| |Image:pRF Tutorial Gray with CorAnal.jpg |}

Make sure the time series are transformed from the Inplane view to the Volume view. This has already been done for you in the example data set, but if you want to re-do it, select the menu '''Xform''' | '''Inplane -> Volume''' | '''tSeries (all scans)''' | '''Trilinear Interpolation'''.

• This step takes the functional data in the Inplane view (which reflects the format in which they were originally collected), applied the Alignment which has alrady been computed for this data set, and interpolates the data to sample it once for every anatomical voxel in the whole-brain T1-weighted anatomy.
• We usually use linear interpolation, in which the values at each point in the anatomy are derived by interpolating spatially to the nearest functional voxels. But we could also use nearest-neighbor interpolation, where the values exactly match the nearest functional voxel.

We will perform our subsequent analyses in this Volume view, so you can close the Inplane window if you want.

== (3) Define the stimulus parameters ==

In the Volume window, select the menu '''Analysis''' | '''Retinotopy Model''' | '''Set Stimulus Parameters'''.

You will get a dialog allowing you to defined what stimulus was used for each scan in this data type. (This dialog actually varies with the version of MATLAB. For MATLAB Version 2007a and newer, a new interface will pop up. Because this new design crashes for MATLAB versions 2006b and older, you will get an older version of the dialog on these versions. The screenshots below are of the older version.)

'''IMPORTANT NOTE''': The pRF analysis always fits all scans in a data type. For this example session, we are analyzing the Averages data, which have two scans: a rotating-wedge scan and a contracting-ring scan. If you are analyzing other data sets, and only want to analyze a subset of your scans, use the option '''Edit''' | '''Data Type''' | '''Group Scans''' to copy the scans you want to analyze together into their own data type.

• It is also important to make sure you set every scan in the data type. Note the slider at the top of the GUI which takes you across the scans.

Set the stimulus parameters for each scan.

• At the dialog top, there is a scan slider which lets you step through scans in the data type. Make sure you set appropriate parameters for each scan in the data type.
• '''Stimulus Type''': Defines the type of stimulus. For the first scan, we select "Wedge", indicating rotating wedges. For scan 2, this is "Ring", indicating expanding/contracting ring.
• ''''Stimulus Radius (deg)''':
• ''''Stimulus Width(deg)''':
• ''''Stimulus Starting Phase(deg)''':
• ''''Stimulus Direction''':
• ''''Stimulus Cycles (#)''':
• ''''Mean-Luminance Blocks''':
• ''''Stimulus Repetitions(#)''':
• ''''Removed Time Frames with Stimuli(deg)''':
• ''''DCT Frequency Max for Detrending(cycles/scan)''':
• ''''HRF Type''':
• ''''HRF Parameters''':
• ''''Frame Interval (sec)''':
• ''''Time Frames (#)''':
• ''''Flip Stimuli Left/Right''':
• ''''Image File''':
• ''''Params File''':
• ''''Image Filter''':

The parameters for the sample data set should already be set appropriately. (The DCT Frequency Max may be set to zero; this is fine, and turns off all detrending except for the constant-offset DC componenet. We actually usually use detrending up to 3 cycles/scan, so you may want to set 3 as this parameter.) {| |Image:pRF Tutorial Stimulus Params Scan 1.jpg | |Image:pRF Tutorial Stimulus Params Scan 2.jpg |}

Press 'OK' on the dialog. The code should now generate the stimulus representation.

• For those who want to get their hands dirty, both the stimulus representation, and the pRF analysis parameters, are stored in the variable '''VOLUME{1}.rm.retinotopyParams'''. You can access these in an object-ish fashion using the syntax params = viewGet(VOLUME{1}, 'rmparams'); .
• params has two sub-fields: params.analysis has the fields describing how you would like to analyze the data, such as the sampling grid in visual space and flags indicating whether you want to do a coarse-to-fine fit. The params.stim field contains a description of the stimulus, as well as some things which are tightly related to the stimulus. This includes a description of the hemodynamic response function (HRF) you use to map from the stimulus and pRF to a predicted BOLD response, and details about the temporal frame rate of the functional data you are analyzing. There is one entry in params.stim for each scan in the data type. So for this session, params.stim(1) describes the wedge scan, and params.stim(2) describes the ring scan.
• When you make the stimuli (which was done automatically when you pressed 'OK' on the dialog), the final resulting images get stored in params.analysis.allstimages. This contains a description of the images concatenated across all scans.

Now, let's look at a movie of the resulting stimulus, to make sure it matches what we showed during the session.

• Select the menu '''Plots''' | '''Retinotopy Model''' | '''View Stimulus Aperture'''.

{| |Image:pRF Tutorial stimulus example.png |}

You can page through the time dimension (in units of MR frames) using the slice slider, or play through the stimulus using the menu option '''DisplayVol Options''' | '''Play through slices'''. The 'Time Course' button lets you browse the temporal pattern of stimulation for different pixels.

== (4) Estimate the pRF parameters for the data ==

Solving the pRF model is generally the longest step in the analysis. On our modern process servers, solving the model can take around 3-5 hours. It can take longer on slower machines.

A set of pRF models have already been run and saved for this example session, so if you're eager to learn the methods and don't want to wait for it to run, you can skip this step and load the included model files.

You have many choices to solve the model. We'll use the default choice, which is to do a full coarse-to-fine fitting for each voxel, and run it within our current instance of MATLAB. The default model we will solve is a circular 2D Gaussian for each voxel.

• To do this, select the menu '''Analysis''' | '''Retinotopy Model''' | '''Run'''.
• If you are working on a linux machine, an alternate way to apply the same sort of model would be to select '''Analysis''' | '''Retinotopy Model''' | '''Run Process in Background (spawn)'''. This option "spawns" a new example of MATLAB as a background process, and runs the same sort of model we're now solving in the current MATLAB window. You don't get the detailed feedback on how things are going, but you can close this MATLAB session, and even log out, and the pRF solution will still be running in the background.

When the model finishes, it will produce three '''.mat''' files containing three iterations of the model solution. These files will all be located in the data directory ('Gray/Averages/'), will start with the text 'retModel-' followed by a time stamp indicating when they were solved, and will have the following suffixes at the end:

• '''retModel-*-gFit.mat''' represents the grid-fit stage. This stage take a discrete "grid" of pre-set pRF parameters (x0, y0, sigma) and fits them to each voxel. There are approximately 100,000 different sets of parameters it generates. It selects the best-fitting set of parameters for each voxel. ** Because we're using a coarse-to-fine fit (the default mode in the Gray/Volume views), this fitting is actually on a version of the time series that is spatially-smoothed along the cortical surface. The motivation is to have voxels which are nearby in cortical space to have similar initial pRF estimates. In the refinement stage, we will not use this spatial smoothing, and will fit the original data.
• '''retModel-*-sFit.mat''' represents a search fit. This is a constrained optimization that is run for each voxel, where the set of pRF parameters (x0, y0, sigma) for each voxel can assume any value to best fit the data (within certain constraints of the search space--it avoids solutions which are wildly different from the coarse fit).
• '''retModel-*-fFit.mat''' represents the final fit of the model. It contains the pRF estimaes from the search-fit, but recomputes the measurement of proportion variance explained for each voxel. Whereas the search fit file estimated this based on a refinement of the grid fit, the '''-fFit''' file uses the final pRF estimate for each voxel to explain the original data.

It is the '''-fFit''' file which we usually analyze, and which we will look at in subsequent sections.

== (5) Load a solved model ==

Select the menu item '''File''' | '''Retinotopy Model''' | '''Load and Select Model File'''.

You will get a file dialog. Select the file '''retModel-20091118-192911-fFit.mat'''. This is the final-fit retinotopy model for this data set.

As the code loads the retinotopy model, it loads four data fields from the model into the four data slots used by mrVista. (See the Retinotopy_Tutorial to understand how these fields are normally used.). By default, the following data are loaded into the following fields:

• The '''co''' slot has the ''variance explained'' for each voxel. This is the proportion of that voxel's variance in the time series explained by the best-fit pRF model (R^2).
• The '''ph''' slot has the ''polar angle'' for each voxel. This is the angle, in radians clockwise from the positive X axis ("3-o-clock"), of each voxel's pRF center. ** '''Note''': There is a (very annoying) Y-axis flip that has consistently existed in all pRF models. Unfortunately, the accessor code for the pRF model has a lot of post-hoc corrections for this, so undoing the bug requires a lot of concomittant changes in the code. For the time being, we have not checked in the full set of fixes, since the existing code is internally consistent. However, one consequence of this flip is that the polar angle moves CLOCKWISE from 3-o-clock, instead of counter-clockwise, as it normally should.
• The '''amp''' slot has the ''pRF size'' for each voxel. This is expressed as the standard deviation of the 2D circular Gaussian function used for this pRF model. ** For some more complex models, we might have multiple Gaussians in different directions. In this case, the default code loads the sigma of the Gaussian in the larger direction.
• The '''map''' slot has the ''eccentricity'' of each voxel. This is the distance, in visual degrees, of the pRF center from the fovea.

{| |Image:pRF Tutorial Gray View Eccentricity.jpg | |Image:pRF Tutorial Gray View Polar Angle.jpg |- |Image:pRF Tutorial Gray View pRF Size.jpg | |Image:pRF Tutorial Gray View Variance Explained.jpg |}

=== Loading individual data fields ===

The pRF has many data fields, and you can load these in whichever of the '''co''', '''ph'''', '''amp''' or '''map''' slots you like. To do this, select the menu '''File''' | '''Retinotopy Model''' | '''Load Model Parameter'''. You will get a dialog allowing you to choose a data field, and select which data field to load it in.

{| |Image:rmLoadParameter dialog.jpg |}

=== Accessing the solved model ===

You can access the saved parameters using the '''viewGet''' and '''viewSet''' formalisms. To get the parameters for the currently-loaded retinotopy model, use the command params = viewGet(VOLUME{1}, 'RMParams');

The params structure is described above. We can also get the model itself, by using:

model = viewGet(VOLUME{1}, 'RMModel');

This produces a cell-array of models, each of which has a model structure. This design allows a given model file to have several different models applied (for instance, single 2-D Gaussian, or a 2-D Gaussian with a second Gaussian surround), although in the models we analyze here, we only ever have a single, 2D circular Gaussian model.

The model structure has several key subfields:

• '''model.x0''': x-position of the pRF center for each voxel;
• '''model.y0''': y-position of the pRF center for each voxel;
• '''model.sigma.major''': size (standard deviation, in visual degrees) of the Gaussian in the larger direction;
• '''model.sigma.minor''': size (standard deviation, in visual degrees) of the Gaussian in the smallerdirection;
• '''model.sigma.theta''': angle between the major and minor axes of the Gaussian

Because the model is a 2D Gaussian, sigma major and sigma minor are the same, and theta is 0 (although this value doesn't matter in the circular case).

To access these, it is recommended to use the '''rmGet''' and '''rmSet''' formalisms; for instance:

x0 = rmGet(model{1}, 'x0'); sigma = rmGet(model{1}, 'sigmamajor'); sigma_angle = rmGet(model{1}, 'theta');

Using this accessor, you can also easily get derived parameters such as:

ecc = rmGet(model{1}, 'eccentricity'); pol = rmGet(model{1}, 'polar-angle'); % angle in radians CW from 3-o-clock varexp = rmGet(model{1}, 'variance explained'); % proportion of each voxel's time course variance explained by the model

== (6) View data on the gray matter surface ==

To view the data maps, we'll load some Mesh. Some of these are pre-made and included with the sample data set.

Viewing these meshes works on most Windows and linux MATLAB installations, but may not work on Macs (the VTK code which renders the meshes has not yet been compiled for Mac).

• Let's start with the right hemisphere. Click the crosshairs in the right hemisphere of cortex.
• Type '''mrmStart''' . A mrMesh server window should appear; you can minimize this window, but keep it open.
• Select the menu '''Gray''' | '''Surface Mesh''' | '''Load and Display''' (or press Control-Comma).
• The mesh files are .mat files, saved separately for each hemisphere (in the folders 3DAnatomy/Left/3DMeshes/, and 3DAnatomy/Right/3DMeshes/). The code starts you off in the hemisphere corresponding to the crosshairs position. By default, this should be 3DAnatomy/Right/3DMeshes/, but if not, you can navigate there.
• Select the file '''ras_right_new_football.mat'''.

You should then get a window opened which looks like this:

{| |Image:pRF Tutorial RH Mesh Example 1.png |}

You should be able to rotate and zoom the mesh using the left and right mouse buttons, respectively; and move the mesh using Ctrl-left click. See the Mesh page for more background.

Now we will project some data maps onto the mesh. Let's start with the eccentricity map.

• Make sure you're viewing the eccentricity map: if not, select '''View''' | '''Parameter Map''' in the Volume window.
• Select '''Gray''' | '''Update Mesh''' (or press Control-0).
• '''NOTE''': If you get an error saying ??? Index exceeds matrix dimensions you can fix it by running: '''Gray / Mesh Display Options / Recompute vertex/gray map'''. You can just call this once from the drop down menu. Then you should not have to call it again.

{| |Image:pRF Tutorial RH Mesh Example Ecc.png |}

Notce how we have a nice progression from the red foveal representation, near the occipital pole, to the peripheral representation more anteriorly?

Now we'll view the polar angle.

• Select '''View''' | '''PhaseMap''' in the Volume window.
• Select '''Gray''' | '''Update Mesh''' (or press Control-0).

{| |Image:pRF Tutorial RH Mesh Example Pol 1.png |}

The boundaries between visual field maps should be fairly clear in this hemisphere. The color map that is loaded by default emphasizes these reversals, at representations of the upper vertical meridian (UVM), lower vertical meridian (LVM), and contralateral horizontal meridian (CHM) -- in this case, the contralateral visual field is the left visual field. This color map lacks the nice red/green, blue/yellow opponency of maps used by packages like Freesurfer (although that color map is also available). But it does have this easy-to-remember mnemonic:

• '''Blue''' is '''up''' like the sky;
• '''Green''' is '''horizontal''', like grass;
• '''Red''' is '''down''', like the fires of Hell.

We'll go ahead and select V1 based on these reversals. V1 runs along the calcarine sulcus, from a dorsal LVM representation, to a ventral UVM representation. To select V1, first select the mrMesh window with the mesh in it, then:

• Rotate the mesh so you can get a clear shot of the reversals (see the images here);
• Press the '''d''' key to enter '''d'''raw mode. The mouse pointer will turn into a pencil icon.
• Click on vertices of the V1 boundaries. For now, we'll start at the peripheral (anterior) part of the blue UVM representation, work our way back towards the back of the brain, where it meets the red LVM representation, then move along the LVM forward back towards the periphery. (Note that the LVM is a little messy in the mid-eccentricities.)
• When you've finished along the LVM, press the '''c''' button to '''c'''lose the polygon. The last point will be connected to the first point; this should make a line along the anterior end of V1.
• We need to fill in the mesh ROI; right now we've only selected the perimeter. Press the '''f''' key to enter '''f'''ill mode. The mouse pointer will turn into a paintbox icon. Click inside the V1 polygon. ** Note that if you click outside the V1 polygon, all parts of the mesh ''except'' those inside V1 will be selected.
• Grab this polygon from the mesh. Go back to the Volume window. Select the menu '''Gray''' | '''Mesh ROIs''' | '''Get ROI from mesh (drawn with "d" key), all Layers'''.
• This will create a new ROI in the Volume window. ** If we'd selected the 'layer 1 only' option, the mesh ROI would only include the first gray matter layer above the gray/white interface. For some situations and segmentations, this might be desirable. For now, we'll take all the gray layers, acknowledging that if our segmentation is incorrect, we may accidentally be including voxels which extend outside the gray matter.

{| |Image:pRF Tutorial RH Mesh Example Pol 2.png | |Image:pRF Tutorial RH Mesh Example Pol 3.png |}

== (7) View ROI data ==

Having loaded a model, and having defined ROIs, we can perform several analyses on the pRF estimates in each ROI. Most of these options are accessible through the menu '''Analysis''' | '''Retinotopy Model'''.

=== Browse the pRFs / time courses within an ROI ===

Select '''Plot Receptive Field and Model Fit (all time points)'''. A separate window will open up, with the '''rmPlotGUI''' tool for the selected ROI.

{| |Image:rmPlotGUI Example.png |}

('''NOTE:''' If you were to select the other option, '''Plot Receptive Field and Model Fit (single cycle)''', you would get the same interface. This is because, when we set the stimulus parameters above, we defined the stimulus as having a single repetition, with 6 rotation cycles. If we didn't have the 4 mean-luminance periods, we could have set the number of repetitions to 6, to match the stimulus rotation cycles. In that case, selecting the single-cycle option will show the mean single-cycle time series from each scan.)

There are several parts to this interface:

• '''Voxel Slider''': selects the current voxel whose data you want to view.
• '''Receptive Field''': Shows the best-fit pRF for this voxel, as saved in the currently-loaded retinotopy model file.
• '''Time Series''': Shows three time series traces: (1) The observed time series for the voxel (black), after trend removal and conversion to percent signal change; (2) the predicted response produced by the estimated pRF (blue); (3) the residuals (observed - prediction). Each of these traces can be toggled using the checkboxes to the right of the plot.
• '''pRF stats''': Shows the parameters for the selected voxel's pRF, and the coordinates of the selected voxel. The coordinates are specified (for this Gray/Volume ROI) as [axi cor sag], specifying the axial (S -> I), coronal (A -> P), and sagittal (L -> R) slice of the selected voxel.