Fitting a 2D population receptive field model to fMRI data from a visual experiment

Author: Malte Lüken (m.luken@esciencecenter.nl)

This examples shows how to fit a population receptive field (pRF) model to blood oxygenation level-dependent (BOLD) functional magnetic resonance imaging (fMRI) data.

A pRF model maps neural activity in a brain region of interest (ROI; e.g., V1 in the human visual cortex) to an experimental stimulus (e.g., a bar moving through the visual field). Here, we use the visual domain as an example, where the pRF is the part of the visual field that stimulates activity in the region of interest. Because the visual field is two-dimensional, the pRF model also has two dimensions.

Because prfmodel uses Keras for model fitting, we need to make sure that a backend is installed before we begin. In this example, we use the TensorFlow backend.

import os
from importlib.util import find_spec

# Set keras backend to 'tensorflow' (this is normally the default)
os.environ["KERAS_BACKEND"] = "tensorflow"
# Hide tensorflow info messages
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "1"

if find_spec("tensorflow") is None:
    msg = "Could not find the tensorflow package. Please install tensorflow with 'pip install .[tensorflow]'"
    raise ImportError(msg)

Loading the surface

Before we start modelling, we need to make sure that we have all the requirements to visualize the BOLD response data. To visualize fMRI data, we can project it onto an anatomical surface of the brain. In this tutorial, we use a template surface based on the surface of the Human Connectome Project (see https://doi.org/10.6084/m9.figshare.13372958). For visualization, we use the nilearn package. We download the surface and load the flat meshes for both hemispheres.

from prfmodel.examples import load_surface_mesh

# Download surface mesh and load as nilearn.surface.PolyMesh object
mesh = load_surface_mesh(dest_dir="data", surface_type="flat")

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1773141557.993844    4067 port.cc:153] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.
I0000 00:00:1773141557.994370    4067 cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used.
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1773141559.169837    4067 port.cc:153] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.
I0000 00:00:1773141559.170145    4067 cudart_stub.cc:31] Could not find cuda drivers on your machine, GPU will not be used.

We can visualize the flat surface with nilearn.

%matplotlib inline
import matplotlib.pyplot as plt
from nilearn.plotting import plot_surf

SURF_VIEW = (90, 270)

fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(8, 6))

plot_surf(mesh, hemi="both", view=SURF_VIEW, axes=ax, figure=fig)

fig.suptitle("Flat surface (HCP template)");

Loading the BOLD response

Now that we have the flat surface, we load the raw BOLD response data from a single subject that can be represented on the surface mesh. We load data from both hemispheres of the subject by setting hemisphere="both".

from prfmodel.examples import load_single_subject_fmri_data

response_raw = load_single_subject_fmri_data(dest_dir="data", hemisphere="both")
response_raw.shape  # shape (num_vertices, num_frames)

(118584, 120)

The BOLD response time courses have 120 time frames. We convert the raw BOLD response signal to percent signal change (PSC).

response_psc = ((response_raw.T / response_raw.mean(axis=1)).T - 1.0) * 100.0

The response object contains the BOLD response timecourse for each voxel on the surface mesh. It has shape (num_vertices, num_frames) where num_vertices is the number of vertices on the surface mesh and num_frames the number of time frames of the recording. The timecourses from the left and right hemispheres are concatenated (left is first). Because the data has been converted to PSC, each timecourse should have a mean of zero.

First, we define a helper function to plot a statistic on the surface.

import numpy as np
from nilearn.plotting import plot_surf_stat_map


def plot_surf_stat_map_helper(
        stat_map: np.ndarray,
        vmin: float | None = None,
        vmax: float | None = None,
        title: str | None = None,
    ) -> tuple[plt.Figure, plt.Axes]:
    """Helper function to plot a surface with a stat map."""
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(8, 6))

    plot_surf_stat_map(
        mesh,
        stat_map,
        vmin=vmin,
        vmax=vmax,
        hemi="both",
        view=SURF_VIEW,
        cmap="inferno",
        axes=ax,
        figure=fig,
        title=title,
    )

    # Expand the 3D axes to fill the space to the left of the colorbar
    surf_axes = [a for a in fig.axes if a.name == "3d"]
    cbar_axes = [a for a in fig.axes if a.name != "3d"]
    cbar_x0 = min(a.get_position().x0 for a in cbar_axes)
    for a in surf_axes:
        a.set_position([0.0, 0.0, cbar_x0 - 0.01, 0.97])

    return fig, ax

We can then plot the standard deviation of each timecourse on the surface.

# Calculate standard deviation of each timecourse
response_sd = response_psc.std(axis=1)

plot_surf_stat_map_helper(response_sd, vmax=5.0, title="Response standard deviation");

We can see that there is high variation in the signal in the visual areas (e.g., V1-V3) as we would expect for a visual stimulus.

Creating the stimulus

For pRF modeling, we also need the experimental stimulus that the subject has seen during the fMRI recording. The stimulus was created using MATLAB so, we first load the raw design matrix that is located in the data directory.

Note - stimulus preprocessing: In this tutorial, we preprocess the raw design matrix before applying the pRF model to it. However, while we use one possible preprocessing approach, we want to mention that there are alternative approaches resting on different assumptions that can lead to similar modelling results.

from scipy.io import loadmat

design_path = "data/vis_design.mat"
design = np.transpose(loadmat(design_path)["stim"])
design.shape  # (num_frames, num_x, num_y)

(120, 500, 281)

The first dimension is time which matches the time frames (second dimension) of our BOLD timecourses. The other dimensions correspond to the x- and y-coordinates of the pixels on the screen on which the stimulus was presented. We can see that the screen was rectangular. For pRF modelling, we prefer a square design, so we transform the design matrix into square shape by padding the y-axis.

# Calculate difference in number of pixels between x and y dimension
pixel_offset = int((design.shape[1] - design.shape[2]) / 2)

# Pad the y-dimension with zeros to match x dimension
design_square = np.pad(
    design,
    [
        (0, 0),  # We ignore the time frame and x dimension
        (0, 0),
        (pixel_offset, pixel_offset+1),
    ],
    constant_values=0,  # Pad with zeros
)
design_square.shape

(120, 500, 500)

Let’s take a look at the range of values in the design.

design_square.min(), design_square.max()

(0.0, 255.0)

We can see that the design still holds raw pixel color intensity values ranging from 0 to 255. For pRF modelling, we first binarize them to 0 (stimulus absent) and 1 (color intensity > 0; stimulus present). Then, we apply a Gaussian filter to smooth the transitions between stimulus absence and presence. Note that we also swap the x- and y-dimension to match the order that prfmodel expects.

from scipy.ndimage import gaussian_filter

design_smooth = np.moveaxis(
    gaussian_filter(
        np.clip(design_square, 0, 1),  # Clip values to interval 0 and 1
        sigma=10,
        axes=(1, 2),  # Apply filter to each time frame separately
    ),
    1, 2,  # Swap x and y dimension because prfmodel expects y dimension to be first!
)
design_smooth.min(), design_smooth.max()

(0.0, 1.0)

Now, the design values range between 0 and 1. However, it still has 500 pixels in each screen dimension. This resolution is unecessarily high for our pRF model, so we only select every 5th pixel to get a lower resolution.

design_smooth_low_res = design_smooth[:, ::5, ::5]
design_smooth_low_res.shape

(120, 100, 100)

The pRF is defined on the visual field and not on the screen pixels. Therefore, we also need to define the coordinates of each pixel in the visual field. The visual field of our stimulus ranges from -10 to 10 degrees of visual angle.

num_grid_points = design_smooth_low_res.shape[1]

x_min, y_min, x_max, y_max = -10.0, -10.0, 10.0, 10.0

coordinates = np.stack(
    # Get all combinations of x- and y-coordinates
    np.meshgrid(
        np.linspace(x_min, x_max, num_grid_points),
        np.linspace(y_min, y_max, num_grid_points),
    ),
    axis=-1,
)
coordinates.shape

(100, 100, 2)

The coordinates object contains the x- and y-coordinate for every screen pixel in the visual field. We combine the final design and the coordinates in a prfmodel.stimuli.PRFStimulus container that we will later use as input for our pRF model.

from prfmodel.stimuli import PRFStimulus

stimulus = PRFStimulus(design_smooth_low_res, coordinates, ["y", "x"])  # y comes before x!

print(stimulus)

PRFStimulus(design=array[120, 100, 100], grid=array[100, 100, 2], dimension_labels=['y', 'x'])

When printing the stimulus object, we can see the shapes of its attributes.

We can visualize the stimulus using prfmodel.stimuli.animate_2d_prf_stimulus().

from IPython.display import HTML
from prfmodel.stimuli import animate_2d_prf_stimulus

ani = animate_2d_prf_stimulus(stimulus, interval=50)  # Pause 50 ms between time frames

HTML(ani.to_html5_video())

We can see that the stimulus consists of four rectangular bars that move vertically and horizontally across the screen.

Inspecting the data

Before we model the BOLD response, we look at the timecourses and inspect the quality of the data. To do this we select a subset of vertices and plot their BOLD response over time. Note that the unit of time frames is repetition time (TR) which is 1.5 seconds for this dataset.

import plotly.io as pio
import plotly.express as px

pio.renderers.default = "notebook_connected"  # Requires internet connection to work
pio.templates.default = "simple_white"

fig = px.line(
    response_psc[::500, :].T,
    animation_frame="variable",
    range_x=(0, 120),
    range_y=(-5, 5),
    labels={
        "index": "Time frame (in TR)",
        "value": "BOLD response (in PSC)",
        "variable": "Vertex",
    },
    title="Vertex time courses",
)
fig.update_layout(showlegend=False, height=450)
fig.show()

While the timecourses are quite noisy, we can see that, for some vertices, there are four regular peaks in the signal. These peaks correspond to the bars moving through the visual field. The goal of our pRF model is to predict these peaks as closely as possible. By comparing how similar the pRF model predictions are to the observed timecourses, we can identify vertices and areas of the brain that respond to our visual stimulus. This allows us to create a stimulus-specific pRF map of the brain.

We can get an even better overview by plotting all timecourses at once in a heatmap.

aspect_ratio = response_psc.shape[1] / response_psc.shape[0]

fig, ax = plt.subplots(1, 1, figsize=(6, 6))

# We use matplotlib because plotly cannot handle this many vertices
im = ax.imshow(
    response_psc,
    aspect=aspect_ratio,
    cmap="inferno",
    vmin=-2,
    vmax=5,
)

ax.set_xlabel("Time frame (in TR)")
ax.set_ylabel("Vertex index")
fig.colorbar(im, ax=ax, label="BOLD response (in PSC)");

Again, we can see the four peaks in the BOLD response of some vertices in both hemispheres. However, the location of the peaks in the timecourse differs between vertices and hemispheres.

Defining the pRF model

Now that we have our BOLD response data and stimulus in place, we can create a pRF model to predict a response to this stimulus. We use the most popular pRF model that is based on the seminal paper by Dumoulin and Wandell (2008): It assumes that the stimulus (our moving bar) elicits a response that follows a Gaussian shape in two-dimensional visual space. This response is convolved with an impulse response that follows the shape of the hemodynamic response in the brain. Finally, a baseline and amplitude parameter shift and scale our predicted response to match the observed BOLD response.

The prfmodel.models.gaussian.Gaussian2DPRFModel class performs all these steps to make a combined prediction. However, we need to add a custom impulse response model to account for the fact that each time frame is one TR (1.5 seconds; the default in prfmodel is 1.0 seconds). Thus, we set the resolution of our predicted impulse response to the TR.

from prfmodel.models.gaussian import Gaussian2DPRFModel
from prfmodel.models.impulse import DerivativeTwoGammaImpulse

# Define repetition time (TR)
tr = 1.5

# Create custom impulse response model
impulse_model = DerivativeTwoGammaImpulse(
    resolution=tr,
    offset=tr / 2.0,
)

We can visualize the predicted impulse response for a set of default parameters.

import pandas as pd

# Define default parameters for impulse response
impulse_default_params = pd.DataFrame({
    "delay": [6.0],
    "dispersion": [0.9],
    "undershoot": [12.0],
    "u_dispersion": [0.9],
    "ratio": [0.48],
    "weight_deriv": [-0.5],
})

# Predict impulse response with default parameters
impulse_response = np.asarray(impulse_model(impulse_default_params))

fig = px.line(
    pd.DataFrame({
        "Time frame (in TR)": np.arange(impulse_response.shape[1]),
        "Impulse response": impulse_response[0]
    }),
    x="Time frame (in TR)",
    y="Impulse response",
    title="Predicted impulse response",
)
fig.update_layout(height=450)
fig.show()

E0000 00:00:1773141648.895891    4067 cuda_platform.cc:52] failed call to cuInit: INTERNAL: CUDA error: Failed call to cuInit: UNKNOWN ERROR (303)

We insert the impulse model into the prfmodel.models.gaussian.Gaussian2DPRFModel that makes combined model predictions.

# Define pRF model with custom impulse response submodel
prf_model = Gaussian2DPRFModel(
    impulse_model=impulse_model,
)

We define a set of starting parameters to make a combined prediction with our pRF model.

# Combine pRF starting parameters with impulse response default parameters
start_params = pd.concat([pd.DataFrame(
    {
        "mu_x": [0.0],
        "mu_y": [0.0],
        "sigma": [1.0],
        "baseline": [0.0],
        "amplitude": [1.0],
    },
), impulse_default_params], axis=1)

# Make prediction with pRF model
simulated_response = np.asarray(prf_model(stimulus, start_params))

fig = px.line(
    pd.DataFrame({
        "Time frame (in TR)": np.arange(simulated_response.shape[1]),
        "BOLD response (in PSC)": simulated_response[0]
    }),
    x="Time frame (in TR)",
    y="BOLD response (in PSC)",
    title="Predicted timecourse",
)
fig.update_layout(height=450)
fig.show()

Some vertices do not have valid timecourses, that is, they contain NaN values. We filter out all vertices that have missing BOLD response measurements.

response_is_valid = np.all(np.isfinite(response_psc), axis=1)
response_valid = response_psc[response_is_valid]

Fitting the pRF model

We will fit the pRF model to our BOLD response data using two stages. We begin with a grid search to find good values for our parameters of interest (mu_x, mu_y, and sigma). Then, we use least squares to estimate the baseline and amplitude of our model.

Note: Typically, we would also use stochastic gradient descent (SGD) to finetune our model fits after the least-squares stage. To keep the example computationally simple, we omit this step here but encourage readers to explore this approach at the end.

Let’s start with the grid search by defining ranges of mu_x, mu_y, and sigma that we want to construct a grid of parameter values from. For baseline and amplitude, we only provide a single value so that they will stay constant across the entire grid. The parameters of the HRF model also remain constant at their default value.

param_ranges = {
    "mu_x": np.linspace(-10, 10, 20),  # Range of visual field in experiment
    "mu_y": np.linspace(-10, 10, 20),
    "sigma": np.linspace(0.005, 10.0, 20),
    "delay": [6.0],
    "dispersion": [0.9],
    "undershoot": [12.0],
    "u_dispersion": [0.9],
    "ratio": [0.48],
    "weight_deriv": [-0.5],
    "baseline": [0.0],
    "amplitude": [1.0],
}

For all three parameters, we defined ranges of 20 values that will be used to construct the grid. That is, the grid search will evaluate all possible combinations of these values and return the combination that fits the simulated data best. This will result in a grid containing \(20 ^ 3 = 8000\) parameter combinations. This is still a relatively small grid and we recommend specifying finer grids in practice.

Let’s construct the prfmodel.fitters.grid.GridFitter and perform the grid search. Note that we set batch_size=20 to let the prfmodel.fitters.grid.GridFitter evaluate 20 parameter combinations at the same time (which saves us some memory). As the loss (i.e., the metric to minimize between model predictions and data), we use cosine similarity which ignores differences in scale between model predictions and observed data.

from keras.losses import CosineSimilarity
from prfmodel.fitters.grid import GridFitter

# Create grid fitter object
grid_fitter = GridFitter(
    model=prf_model,
    stimulus=stimulus,
    loss=CosineSimilarity(reduction="none"),  # Grid fitter needs no aggregation of loss
)

# Run grid search
grid_history, grid_params = grid_fitter.fit(
    data=response_valid,
    parameter_values=param_ranges,
    batch_size=20,
)

We can print the parameters estimated by the grid search.

grid_params

	mu_x	mu_y	sigma	delay	dispersion	undershoot	u_dispersion	ratio	weight_deriv	baseline	amplitude
0	-10.000000	-10.000000	0.005000	6.0	0.9	12.0	0.9	0.48	-0.5	0.0	1.0
1	-5.789474	10.000000	0.005000	6.0	0.9	12.0	0.9	0.48	-0.5	0.0	1.0
2	7.894737	-1.578947	10.000000	6.0	0.9	12.0	0.9	0.48	-0.5	0.0	1.0
3	8.947369	0.526316	0.531053	6.0	0.9	12.0	0.9	0.48	-0.5	0.0	1.0
4	5.789474	10.000000	2.635263	6.0	0.9	12.0	0.9	0.48	-0.5	0.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
93529	8.947369	-10.000000	0.531053	6.0	0.9	12.0	0.9	0.48	-0.5	0.0	1.0
93530	8.947369	-5.789474	0.531053	6.0	0.9	12.0	0.9	0.48	-0.5	0.0	1.0
93531	5.789474	10.000000	0.005000	6.0	0.9	12.0	0.9	0.48	-0.5	0.0	1.0
93532	6.842105	-10.000000	0.531053	6.0	0.9	12.0	0.9	0.48	-0.5	0.0	1.0
93533	5.789474	-10.000000	0.005000	6.0	0.9	12.0	0.9	0.48	-0.5	0.0	1.0

93534 rows × 11 columns

We can see that the estimates for mu_x, mu_y, and sigma are one combination in our grid. In the second step, we also optimize the amplitude of the pRF model together with the baseline using least squares. This adjusts the scale of our model predictions to scale the observed data. We set batch_size=200 to estimate least-squares fits for batches of vertices sequentially and save memory.

from prfmodel.fitters.linear import LeastSquaresFitter

# Create least-squares fitter
ls_fitter = LeastSquaresFitter(
    model=prf_model,
    stimulus=stimulus,
)

# Run least squares fit
ls_history, ls_params = ls_fitter.fit(
    data=response_valid,
    parameters=grid_params,
    slope_name="amplitude",  # Names of parameters to be optimized with least squares
    intercept_name="baseline",
    batch_size=200,
)

ls_params

	mu_x	mu_y	sigma	delay	dispersion	undershoot	u_dispersion	ratio	weight_deriv	baseline	amplitude
0	-10.000000	-10.000000	0.005000	6.0	0.9	12.0	0.9	0.48	-0.5	0.000000e+00	0.000000e+00
1	-5.789474	10.000000	0.005000	6.0	0.9	12.0	0.9	0.48	-0.5	6.801418e-10	1.717191e-16
2	7.894737	-1.578947	10.000000	6.0	0.9	12.0	0.9	0.48	-0.5	-2.090634e-01	4.589847e-01
3	8.947369	0.526316	0.531053	6.0	0.9	12.0	0.9	0.48	-0.5	1.360284e-09	1.788481e-15
4	5.789474	10.000000	2.635263	6.0	0.9	12.0	0.9	0.48	-0.5	-6.851990e-01	1.330586e+00
...	...	...	...	...	...	...	...	...	...	...	...
93529	8.947369	-10.000000	0.531053	6.0	0.9	12.0	0.9	0.48	-0.5	-2.720567e-09	-1.958322e-15
93530	8.947369	-5.789474	0.531053	6.0	0.9	12.0	0.9	0.48	-0.5	2.720567e-09	3.602110e-15
93531	5.789474	10.000000	0.005000	6.0	0.9	12.0	0.9	0.48	-0.5	0.000000e+00	0.000000e+00
93532	6.842105	-10.000000	0.531053	6.0	0.9	12.0	0.9	0.48	-0.5	-1.896294e-02	3.569337e-01
93533	5.789474	-10.000000	0.005000	6.0	0.9	12.0	0.9	0.48	-0.5	-5.441135e-09	-1.813444e-15

93534 rows × 11 columns

We can see that the amplitudes are substantially lower compared to the starting value (and our initial guess) for many vertices.

Now that the core pRF parameters and the auxiliary baseline and amplitude parameters are optimized, we can compare the model predictions against the observed responses. Because we want to make predictions for all vertices in the brain, we wrap our prf_model in the prfmodel.utils.batched() modifier function. The modifier changes the behavior of the model to make predictions for batches of vertices sequentially. This saves us a lot of memory at the expense of minimal runtime overhead.

from prfmodel.utils import batched

predict_batched = batched(prf_model)

# Make predictions with optimized parameters
pred_response = np.asarray(predict_batched(stimulus, ls_params, batch_size=200))

We can quantify how well the predictions align with the observed timecourses using the R-squared metric. This metric indicates the proportion of variance in the observed data explained by our model predictions.

from keras.metrics import R2Score

r2_metric = R2Score(class_aggregation=None)  # Don't aggregate score over vertices

r_squared = np.asarray(r2_metric(response_valid.T, pred_response.T))  # Transpose to compute score across time frames
r_squared.shape

(93534,)

We can look at the distribution of R-squared values across vertices.

fig = px.histogram(
    x=r_squared,
    nbins=15,
    labels={"x": "R-squared"},
).update_layout(yaxis_title="Frequency", height=450)
fig.show()

We can see that many vertices have a score close to zero meaning that the pRF model does not predict the observed response well. This is expected since not the entire brain responds to our relatively simple visual stimulus. However, a substantial amount of vertices also have higher scores, suggesting that some brain areas do respond to it.

Let’s plot the predicted and the observed timecourses for a subsample of vertices.

import plotly.graph_objects as go

each_k = 500

# Sort vertices that are not in V1 according to R-squared
best_vertices = np.flip(np.argsort(r_squared))[::each_k]

response_valid_best_vertices = response_valid[best_vertices]
pred_response_best_vertices = pred_response[best_vertices]
r_squared_valid_best_vertices = r_squared[best_vertices]
sigma_best_vertices = ls_params["sigma"].values[best_vertices]

df_valid = pd.DataFrame(response_valid_best_vertices.T)
df_pred = pd.DataFrame(pred_response_best_vertices.T)

df_valid["source"] = "Observed"
df_pred["source"] = "Predicted"

df = pd.concat([df_valid, df_pred], axis=0)
df["time"] = np.tile(np.arange(df_valid.shape[0]), 2)

df_melted = df.melt(id_vars=["source", "time"], var_name="vertex", value_name="response")

fig = px.line(
    df_melted,
    x="time",
    y="response",
    color="source",
    animation_frame="vertex",
    range_y=[-5, 5],
    labels={
        "time": "Time frame (in TR)",
        "response": "BOLD response (in PSC)",
        "vertex": "Vertex",
        "source": "",
    },
    title=f"Observed and predicted vertex responses (subsampled vertices)",
)

# Add a text trace to display per-vertex stats; this will be updated in each animation frame
fig.add_trace(go.Scatter(
    x=[60],
    y=[4.7],
    mode="text",
    text=[f"R-squared = {r_squared_valid_best_vertices[0]:.3f}, sigma = {sigma_best_vertices[0]:.3f}"],
    showlegend=False,
    hoverinfo="skip",
    textfont=dict(size=13),
))
stats_trace_idx = len(fig.data) - 1

# Append a stats text update to each animation frame's trace data
for i, frame in enumerate(fig.frames):
    r2 = r_squared_valid_best_vertices[i]
    sigma = sigma_best_vertices[i]
    frame.data = list(frame.data) + [go.Scatter(
        text=[f"R-squared = {r2:.3f}, sigma = {sigma:.3f}"],
    )]

fig.update_layout(showlegend=True, height=450)
fig.show()

We can see that the model predicts the observed timecourses for these vertices very well.

We can also visualize the R-squared score on the flat surface mesh.

r_squared_full = np.full((response_psc.shape[0],), fill_value=np.nan)
r_squared_full[response_is_valid] = r_squared

plot_surf_stat_map_helper(r_squared_full, vmin=0.0, vmax=1.0, title="Variance explained (R-squared)");

We can see that the vertices with the highest scores are located in the visual areas of both hemispheres.

Analyzing the pRF results

To analyze and interpret the pRF parameters, we will zoom in on vertices with R-squared > 0.6 which is roughly 7% of the vertices in the brain (note that this threshold is somewhat arbitrary).

# Create mask for vertices above R-squared threshold
is_above_threshold = r_squared > 0.6

# Compute proportion of vertices above threshold
is_above_threshold.mean()

0.07056257617550837

For these vertices, we can look at different quantities to interpret the pRFs estimated by our model.

First, we look at the pRF size indicated by sigma and plot it on the surface.

# Allocate size vector for the entire surface
size = np.full((response_psc.shape[0],), fill_value=np.nan)

# Fill valid vertices with estimated size parameters
size[response_is_valid] = np.where(is_above_threshold, ls_params["sigma"], np.nan)

plot_surf_stat_map_helper(size, vmin=0.0, vmax=10.0, title="pRF size (sigma)");

We can see that vertices in the early visual pathway (e.g., V1) have smaller sizes than those in the higher areas.

Besides the pRF size, we can also look at the position of the pRF relative to the center of the screen. First, we compute the angle of the center of the pRF relative to the center (i.e., the polar angle).

def calc_angle(mu_x: float, mu_y: float) -> float:
    """Compute the polar angle of a pRF from the x- and y-coordinate of its center."""
    return np.angle(mu_x + mu_y * 1j)

angle = np.full((response_psc.shape[0],), fill_value=np.nan)

angle_valid = calc_angle(ls_params["mu_x"], ls_params["mu_y"])
angle[response_is_valid] = np.where(is_above_threshold, angle_valid, np.nan)

plot_surf_stat_map_helper(angle, vmin=-np.pi, vmax=np.pi, title="pRF center polar angle");

The surface plot shows that vertices in the left hemisphere have a positive pRF angle, meaning that their pRF center is on the right side of the screen. In contrast, vertices in the right hemisphere have negative pRF angle, thus, their pRF center is located on the left side of the screen.

We can also look at the eccentricity, that is, the distance of the pRF center from the center of the screen.

def calc_eccentricity(mu_x: float, mu_y: float) -> float:
    """Compute the eccentricity of a pRF from the x- and y-coordinate of its center."""
    return np.abs(mu_x + mu_y * 1j)

eccentricity = np.full((response_psc.shape[0],), fill_value=np.nan)

eccentricity_valid = calc_eccentricity(ls_params["mu_x"], ls_params["mu_y"])
eccentricity[response_is_valid] = np.where(is_above_threshold, eccentricity_valid, np.nan)

plot_surf_stat_map_helper(eccentricity, vmin=0.0, vmax=10.0, title="pRF center eccentricity");

For eccentricity, we can see segments with gradual transitions in both hemispheres that correspond to pRFs that are closer or further away from the center of the screen.

Conclusion

This example showed how to fit a two-dimensional Gaussian pRF model to empirical fMRI data collected from an experiment in the visual domain. First, we plotted the raw BOLD response data on the cortical surface. Second, we created the experimental stimulus. Then, we defined a pRF model and optimized its parameters using a grid search followed by least-squares to adjust for baseline and amplitude differences. Finally, we visualized model fit, the estimated parameters, and derived measures on the cortical surface.

Next steps

The predictions by our pRF model can potentially be improved. We suggest different directions for improving the pRF model fit:

• Increasing the number of points in the parameter grid for the grid search

• Finetuning the pRF model parameters with stochastic gradient descent with prfmodel.fitters.sgd.SGDFitter

• Applying preprocessing steps before fitting the pRF model (e.g., high-pass filtering)

• Optimizing the impulse response model parameters in the grid search

• Building a more complex pRF model (e.g., compressive spatial summation, see Kay et al., 2013)

Stay Tuned

More tutorials on fitting models to empirical data and creating custom models are in the making.

For questions and issues, please make an issue on GitHub or contact Malte Lüken (m.luken@esciencecenter.nl).

References

Dumoulin, S. O., & Wandell, B. A. (2008). Population receptive field estimates in human visual cortex. NeuroImage, 39(2), 647–660. https://doi.org/10.1016/j.neuroimage.2007.09.034

Kay, K. N., Winawer, J., Mezer, A., & Wandell, B. A. (2013). Compressive spatial summation in human visual cortex. Journal of Neurophysiology, 110(2), 481–494. https://doi.org/10.1152/jn.00105.2013

Van Essen, D. C., Ugurbil, K., Auerbach, E., Barch, D., Behrens, T. E. J., Bucholz, R., Chang, A., Chen, L., Corbetta, M., Curtiss, S. W., Della Penna, S., Feinberg, D., Glasser, M. F., Harel, N., Heath, A. C., Larson-Prior, L., Marcus, D., Michalareas, G., Moeller, S., … WU-Minn HCP Consortium. (2012). The Human Connectome Project: A data acquisition perspective. NeuroImage, 62(4), 2222–2231. https://doi.org/10.1016/j.neuroimage.2012.02.018