Skip to content

Fast & real time pRF mapping (Python)

mario.senden edited this page May 11, 2022 · 4 revisions

Fast / real time pRF mapping is a regression-based approach. Specifically, the BOLD signal is predicted from a linear combination of a set of basis functions encoding the effective stimulus. These basis functions are hashed Gaussians. This means that it essentially utilizes randomly placed isotropic Gaussians for encoding the stimulus but rather than using individual Gaussians, several of those are hashed together to form a single feature.

To find the weight of each feature, the tool can either use ridge regression (most useful for an offline scenario) or gradient descent (most useful for an online scenario).

Instantiation

The tool needs to be instantiated with a number of parameters:

  • f_sampling - sampling frequency of data acquisition (1 / TR)
  • n_voxels - number of voxels in data set
  • n_features - number of features
  • n_gaussians - number of Gaussian hashed together for each feature
  • fwhm - full width at half maximum of Gaussian
  • eta - learning rate (inverse of ridge regularization parameter)
parameters = {'f_sampling': sampling_frequency,
           'r_stimulus': stimulus_resolution,
           'n_voxels': number_of_voxels,
           'n_features': number_of_features,
           'n_gaussians': number_of_Gaussians,
           'fwhm': full_width_half_max,
           'eta': learning_rate}

hgr = HGR(parameters)

Usage

Ridge regression & gradient descent To find the weights for each feature in each voxel, it is possible to use ridge regression or gradient descent. The former can quite simply be carried out by calling the ridge function and passing all data and the effective stimulus hgr.ridge(data, stimulus). For gradient descent, the update function needs to be called repeatedly for each time point hgr.update(data[t,:] stimulus[t,:]). Note that in contrast to the pRF tool, the stimulus needs to have dimensions time-by-pixels rather than pixels-by-time. Also note, that in order to use the same instantiation of the class for online analysis of several runs, the instance needs to be reset prior to each run (hgr.reset())

After weight estimation, receptive fields can be obtained by the dot product between a set of features and their weights. Subsequently, it is recommended to normalize the raw receptive field of each voxel and shrink it (using a shrinkage exponent alpha).

features = hgr.get_features()
weights = hgr.get_weights()
RFs = np.dot(features, weights)

mx = np.max(RFs, axis=0)
mn = np.min(RFs, axis=0)
RFs = (RFs - mn) / (mx - mn)
RFs = RFs**alpha

RFs = np.reshape(RFs, (stimulus_resolution, stimulus_resolution, number_of_voxels))

Voxel selection The tool has an in-built voxel selection function which performs blocked cross-validation for time-series data. By default, the function splits the data into 4 splits and returns a binary mask selecting those voxels whose average fit over split exceeds the 95th percentile of the fit distribution over all voxels. Alternatively, the function can select a pre-specified number of voxels or use a pre-specified fit threshold for voxel selection.

mask, corr_fit = hgr.get_best_voxels(data,stimulus,
                                       n_splits = number_of_splits,   # number of splits for cross-validation
                                       type = type_of_cutoff          # cutoff type ('percentile', 'number', 'threshold')
                                       cutoff = cutoff_value)         # cutoff value

pRF parameter estimation After optimizing regression weights, the tool can also be used to obtain pRF parameters using results = hgr.get_parameters(). To speed up this process, the tool can be passed a binary mask as optional input. Specifically, the mapping function takes four optional arguments:

  • n_batch - batch size; i.e. number of voxel estimated in parallel (default = 10000)
  • max_radius - radius of measured visual field (default = 10)
  • alpha - shrinkage parameter (default = 1)
  • mask - a binary mask specifying for which voxels the analysis should be carried out

The function returns a structure (results) with five fields:

  • 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

The most efficient way to use this tool for obtaining pRF parameters is to first perform voxel selection, then perform ridge regression and lastly obtain parameters limited to selected voxels.