I am looking for a way to extract the mean BOLD response over all voxels in an ROI, in order to create a mean time course answer to repeated events. I would use some BIDS data, preprocessed using fMRIprep, something like sub-01_task-01_run01.nii.gz, events.tsv, and mask.nii.gz. The objective is to compare it with canonical HRF. Ideally, I would like to generalize over all runs, maybe even all subjects. And I can only use Python. This would be the perfect plot:
PyHRF seems to do exactly what I’m looking for, but the repo hasn’t been updated for ~10 years, it uses Python2, and the website is down - the manual is still online.
from nilearn import image, plotting
from nilearn.maskers import NiftiLabelsMasker
import matplotlib.pyplot as plt
import numpy as np
# Load your functional image and ROI mask
func_img = "path_to_func.nii.gz" # 4D fMRI data
roi_mask = "path_to_roi.nii.gz" # ROI mask (must be the same space as func_img)
# Extract mean time series from ROI
masker = NiftiLabelsMasker(labels_img=roi_mask, standardize=True)
timeseries = masker.fit_transform(func_img)
# Compute mean across all voxels in the ROI
mean_bold = np.mean(timeseries, axis=1)
# Plot the extracted BOLD response
plt.figure(figsize=(6, 4))
plt.plot(mean_bold, label="Mean BOLD Signal", color="b")
plt.xlabel("Time (scans)")
plt.ylabel("BOLD Signal (a.u.)")
plt.title("Mean BOLD Time Series in ROI")
plt.legend()
plt.show()
Example plot showing different HRFs:
from nilearn.glm.first_level import compute_regressor
# Example stimulus onset times (adjust based on your paradigm)
onset_times = np.array([10, 30, 50, 70]) # Onsets in seconds
duration = 5 # Stimulus duration in seconds
tr = 2 # Repetition time (TR) in seconds
n_scans = len(mean_bold) # Number of time points
# Generate different HRF models
hrfs = ["glover", "spm", "fir"]
plt.figure(figsize=(6, 4))
for hrf_model in hrfs:
hrf_regressor, _ = compute_regressor(
event_onsets=[onset_times], event_durations=[duration],
frame_times=np.arange(n_scans) * tr, hrf_model=hrf_model
)
plt.plot(hrf_regressor, label=f"HRF: {hrf_model}")
plt.plot(mean_bold, label="Mean BOLD", color="k", linestyle="dashed")
plt.xlabel("Time (scans)")
plt.ylabel("Amplitude")
plt.title("HRF Model Comparison")
plt.legend()
plt.show()
