-
Notifications
You must be signed in to change notification settings - Fork 5
Population Receptive Field Mapping (Python)
In the population receptive field mapping approach to retinotopy, receptive fields are formalized as parametric models (e.g. isotropic Gaussian). The parameters for each voxel are fit such that a BOLD signal predicted from moving a stimulus (according to the stimulation protocol of the retinotopy experiment) across the receptive field (RF) model matches the BOLD observed for that voxel. In the pRF
tool provided in the CNI toolbox, a grid search is used to find the best parameters for all voxels in an isotropic 2D Gaussian model. After linear encoding of the stimulus using this model, compressive spatial summation may be applied.
The three parameters fit by the tool are the location of the RF (in Cartesian coordinates) and its size. By default, the grid does not explore a fixed set of size values. Rather it exploits that RF sizes are linearly related to eccentricity and explores a range of slopes for the size-eccentricity relationship. This effectively allows for exploration of a greater range of receptive field sizes. It is, however, also possible to simply explore a fixed set of size values. In terms of RF location, the visual field is by default split into a circular grid of points, whose density decays exponentially with eccentricity. Close to the fovea, the grid is thus denser than in the periphery, thus taking cortical magnification into account. It is, however, also possible with the pRF
tool to generate a linear grid of RF locations.
The tool needs to be instantiated with a number of parameters:
- f_sampling - sampling frequency of data acquisition (1 / TR)
- h_stimulus - stimulus height
- w_stimulus - stimulus width
- n_samples - number of samples (functional volumes)
- n_rows - number of rows (1st volumetric dimension of 4D data tensor)
- n_cols - number of columns (2nd volumetric dimension of 4D data tensor)
- n_slices - number of slices (3rd volumetric dimension of 4D data tensor)
if data.ndim==4:
n_samples, n_rows, n_cols, n_slices = data.shape
elif data.ndim==2:
n_samples, n_rows = data.shape
n_cols, n_slices = 1, 1
parameters = {'f_sampling': sampling_frequency,
'h_stimulus': stimulus_height,
'w_stimulus': stimulus_width,
'n_samples': n_samples,
'n_rows': n_rows,
'n_cols': n_cols,
'n_slices': n_slices}
Before RF parameters can be estimated, it is necessary to generate predicted BOLD signals based on a specific grid. This can be done by first importing a stimulus and then using the generate_timecourses
function.
Importing a stimulus
The tool requires a rank 3 stimulus tensor whose rows and columns reflect stimulus resolution and whose slices reflect time. If such a tensor has already been loaded it can be provided to the tool prf.set_stimulus(stimulus)
Alternatively, the tool can create such a tensor from stimulus image (.png) files by calling prf.import_stimulus()
and selecting the folder containing the stimulus images.
Creating timecourses
By default, the generate_timecourses
function generates a grid of thirty locations in each the x and y dimension as well as 10 slopes. Locations range from -10 to +10 degrees of visual angle (maximum radius of stimulated visual field). Slopes range from 0.1 to 1.2. Instead of exploring slopes, it is also possible to explore sigmas directly. In that case sigmas range from min_sigma
(default = 0.1) to max_sigma * max_radius
, where max_sigma = 1
by default and max_radius
is the maximum radius of the stimulated visual field. Compressive spatial summation is not applied by default.
Using default parameters, creating predicted BOLD signals boils down to prf.create_timecourses();
However, all these parameters can be changed.
prf.create_timecourses(use_slope = true, # flag to indicate whether to explore slopes (true) or sigmas directly (false)
num_xy = number_xy, # number of points in x and y dimension
num_size = number_sizes, # number of points in the size dimension (either slopes or sigmas)
max_radius = max_radius, # maximum radius of stimulated visual field
min_slope = min_slope, # lower bound of slope (only needed when use_slope = true)
max_slope = max_slope, # upper bound of slope (only needed when use_slope = true)
min_sigma = min_sigma, # lower bound of sigma (only needed when use_slope = false)
max_sigma = max_sigma, # upper bound of sigma (only needed when use_slope = False)
css_exponent = alpha, # exponent for compressive spatial summation
sampling = 'log', # location sampling type ('log' or 'linear')
Mapping
The prf
instance can now be used to map pRF parameters for data from runs sharing the same parameters and stimulation protocol specified before: results = prf.mapping(data)
During mapping, voxels whose mean signal intensity falls below a threshold will be skipped.
This threshold can be adjusted. Specifically, the mapping
function takes two optional arguments:
- threshold - a mean signal intensity threshold below which a voxel is skipped (default = 100)
- mask - a binary mask specifying for which voxels the analysis should be carried out (if a mask is used, the threshold is omitted)
The function returns a dictionary (results
) with six keys:
- corr_fit - correlation between observed and best fitting predicted BOLD signal
- mu_x - x-coordinate of RF
- mu_y - y-coordinate of RF
- sigma - RF size
- eccentricity - eccentricity of RF
- polar_angle - polar angle of RF
Values corresponding to each key retain the volumetric dimensions of the data.