Hi @Jon_Brooks ,
You could also just use Nilearn to code a gPPI (of course, will need you to apply the fMRIPrep-generated MNI xfms to your images first). Here is an example created by AI:
import numpy as np
import pandas as pd
from scipy.linalg import toeplitz
from nilearn.glm.first_level import (
FirstLevelModel,
compute_regressor,
make_first_level_design_matrix,
)
from nilearn.glm.first_level.hemodynamic_models import glover_hrf
from nilearn.maskers import NiftiSpheresMasker
# ------------------------------------------------------------
# 1. Inputs
# ------------------------------------------------------------
func_img = "sub-01_task-example_bold.nii.gz"
events = pd.read_csv("sub-01_task-example_events.tsv", sep="\t")
confounds = pd.read_csv("sub-01_task-example_confounds.tsv", sep="\t")
tr = 2.0
n_scans = 240
frame_times = np.arange(n_scans) * tr
seed_xyz = [(0, -52, 26)] # example PCC coordinate
condition_names = ["condition_A", "condition_B", "condition_C"]
# ------------------------------------------------------------
# 2. Extract the physiological seed signal
# ------------------------------------------------------------
seed_masker = NiftiSpheresMasker(
seed_xyz,
radius=6,
detrend=False,
standardize=False,
t_r=tr,
)
seed_bold = seed_masker.fit_transform(
func_img,
confounds=confounds,
).squeeze()
# Centering is important for interpretability and collinearity.
seed_bold = seed_bold - seed_bold.mean()
# ------------------------------------------------------------
# 3. Approximate HRF deconvolution
# ------------------------------------------------------------
def ridge_deconvolve(signal, tr, alpha=1.0):
"""Approximate neural signal using regularized HRF deconvolution."""
n = len(signal)
hrf = glover_hrf(
t_r=tr,
oversampling=1,
time_length=32.0,
onset=0.0,
)
hrf = hrf[:n]
# Convolution matrix H such that BOLD ≈ H @ neural_signal.
first_col = np.r_[hrf[0], np.zeros(n - 1)]
first_row = np.r_[hrf, np.zeros(max(0, n - len(hrf)))]
H = toeplitz(first_col, first_row[:n])
# Ridge-regularized inverse.
lhs = H.T @ H + alpha * np.eye(n)
rhs = H.T @ signal
return np.linalg.solve(lhs, rhs)
seed_neural = ridge_deconvolve(seed_bold, tr=tr, alpha=1.0)
seed_neural -= seed_neural.mean()
# ------------------------------------------------------------
# 4. Construct one neural psychological vector per condition
# ------------------------------------------------------------
def condition_boxcar(events, condition, frame_times):
"""
Create a scan-resolution, unconvolved psychological vector.
For sub-TR event timing, construct this on an oversampled time grid
instead of directly at scan resolution.
"""
vector = np.zeros(len(frame_times), dtype=float)
selected = events.loc[events["trial_type"] == condition]
for onset, duration in zip(selected["onset"], selected["duration"]):
active = (
(frame_times >= onset)
& (frame_times < onset + duration)
)
vector[active] = 1.0
# Mean-centering is standard for the interaction.
return vector - vector.mean()
def convolve_interaction(interaction_neural, frame_times, name):
"""
Convolve a scan-resolution neural interaction using Nilearn's HRF.
"""
onsets = frame_times
durations = np.zeros_like(frame_times)
amplitudes = interaction_neural
experimental_condition = np.vstack(
[onsets, durations, amplitudes]
)
regressor, _ = compute_regressor(
experimental_condition,
hrf_model="glover",
frame_times=frame_times,
con_id=name,
oversampling=1,
)
return regressor[:, 0]
ppi_regressors = {}
for condition in condition_names:
psych_neural = condition_boxcar(
events,
condition,
frame_times,
)
interaction_neural = psych_neural * seed_neural
ppi_regressors[f"ppi_{condition}"] = convolve_interaction(
interaction_neural,
frame_times,
name=f"ppi_{condition}",
)
# ------------------------------------------------------------
# 5. Build the full first-level design matrix
# ------------------------------------------------------------
# Select suitable nuisance variables for your dataset.
confound_columns = [
"trans_x", "trans_y", "trans_z",
"rot_x", "rot_y", "rot_z",
]
nuisance = (
confounds[confound_columns]
.fillna(0)
.to_numpy()
)
extra_regressors = np.column_stack(
[
seed_bold,
*[ppi_regressors[f"ppi_{c}"] for c in condition_names],
nuisance,
]
)
extra_names = (
["seed"]
+ [f"ppi_{c}" for c in condition_names]
+ confound_columns
)
design_matrix = make_first_level_design_matrix(
frame_times=frame_times,
events=events,
hrf_model="glover",
drift_model="cosine",
high_pass=0.01,
add_regs=extra_regressors,
add_reg_names=extra_names,
)
# ------------------------------------------------------------
# 6. Fit the model
# ------------------------------------------------------------
model = FirstLevelModel(
t_r=tr,
noise_model="ar1",
standardize=False,
signal_scaling=0,
minimize_memory=False,
)
model = model.fit(
func_img,
design_matrices=design_matrix,
)
# ------------------------------------------------------------
# 7. Compute gPPI contrasts
# ------------------------------------------------------------
# Condition-specific PPI:
ppi_A_map = model.compute_contrast(
"ppi_condition_A",
output_type="z_score",
)
# Differential connectivity: A > B
ppi_A_gt_B_map = model.compute_contrast(
"ppi_condition_A - ppi_condition_B",
output_type="z_score",
)
and doesn’t work with it either 