Beta-Series Modeling for Task-Based Functional Connectivity

This example shows how to run beta-series GLM models, which are a common modeling approach for a variety of analyses of task-based fMRI data with an event-related task design, including functional connectivity, decoding, and Representational Similarity Analysis (RSA).

First, we compare the standard GLM modeling approach to two common beta-series modeling approaches, Least Squares All (LSA) and Least Squares Separate (LSS). Then, we show how to use the resulting beta-series for a simple task-based functional connectivity analysis.

Prepare data and analysis parameters

Download data in BIDS format and event information for one subject, and create a standard FirstLevelModel.

For more information see the dataset description.

from nilearn.datasets import fetch_language_localizer_demo_dataset
from nilearn.glm.first_level import FirstLevelModel, first_level_from_bids
from nilearn.plotting import plot_design_matrix, plot_stat_map, show

data = fetch_language_localizer_demo_dataset()

[fetch_language_localizer_demo_dataset] Dataset directory found: /home/runner/work/nilearn/nilearn/nilearn_data/fMRI-language-localizer-demo-dataset

Find the first subject’s functional run, confounds, and events file and use them to instantiate a FirstLevelModel.

models, models_run_imgs, events_dfs, models_confounds = first_level_from_bids(
    dataset_path=data.data_dir,
    task_label="languagelocalizer",
    space_label="",
    sub_labels=["01"],
    img_filters=[("desc", "preproc")],
    n_jobs=2,
)

# Grab the first subject's FirstLevelModel object
for _idx, model in enumerate(models):
    standard_glm = model

    # Multiple fMRI runs (and accompanying event files) may
    # support a single subject's FirstLevelModel,
    # so we loop through. Note that for this dataset,
    # we have only one run per subject.
    for run_imgs, events in zip(
        models_run_imgs[_idx], events_dfs[_idx], strict=False
    ):
        fmri_file = run_imgs
        events_df = events

# We can verify that only one fMRI run was found for
# this subject as expected.
print("fmri_file:", fmri_file)

/home/runner/work/nilearn/nilearn/examples/07_advanced/plot_beta_series.py:37: RuntimeWarning:

'StartTime' not found in file /home/runner/work/nilearn/nilearn/nilearn_data/fMRI-language-localizer-demo-dataset/derivatives/sub-01/func/sub-01_task-languagelocalizer_desc-preproc_bold.json.

/home/runner/work/nilearn/nilearn/examples/07_advanced/plot_beta_series.py:37: UserWarning:

'slice_time_ref' not provided and cannot be inferred from metadata.
It will be assumed that the slice timing reference is 0.0 percent of the repetition time.
If it is not the case it will need to be set manually in the generated list of models.

fmri_file: /home/runner/work/nilearn/nilearn/nilearn_data/fMRI-language-localizer-demo-dataset/derivatives/sub-01/func/sub-01_task-languagelocalizer_desc-preproc_bold.nii.gz

Define the standard model

Here, we create a basic GLM for this one run, which we can use to highlight differences between standard event-modeling approach and beta-series models.

We will just use the one created by first_level_from_bids.

import matplotlib.pyplot as plt

# Fit the model.
standard_glm.fit(fmri_file, events_df)

# The standard design matrix has one column for each condition, along with
# columns for the confound regressors and drifts.
fig, ax = plt.subplots(figsize=(5, 10))
plot_design_matrix(standard_glm.design_matrices_[0], axes=ax)
show()

We will reuse the parameters supplied above to first_level_from_bids for all other models, so we can extract these directly from the standard_glm object.

glm_parameters = standard_glm.get_params()

# We need to override one parameter (``signal_scaling``).
glm_parameters["signal_scaling"] = standard_glm.signal_scaling

Define the Least Squares-All (LSA) model

We will now create a Least Squares-All (LSA) model. This involves creating a new condition in the design matrix for each trial of interest. It’s important to ensure that the original conditions can be inferred from the new trial-wise conditions, in order to collect the resulting beta maps into condition-wise beta-series. Here, we will do this using a unique delimiter (__), which should not be present in the original condition names.

# Transform the data frame for LSA
lsa_events_df = events_df.copy()
conditions = lsa_events_df["trial_type"].unique()
condition_counter = dict.fromkeys(conditions, 0)
for i_trial, trial in lsa_events_df.iterrows():
    trial_condition = trial["trial_type"]
    condition_counter[trial_condition] += 1
    # We use a unique delimiter here (``__``) which shouldn't be in the
    # original condition names.
    trial_name = f"{trial_condition}__{condition_counter[trial_condition]:03d}"
    lsa_events_df.loc[i_trial, "trial_type"] = trial_name

lsa_glm = FirstLevelModel(**glm_parameters)
lsa_glm.fit(fmri_file, lsa_events_df)

fig, ax = plt.subplots(figsize=(10, 10))
plot_design_matrix(lsa_glm.design_matrices_[0], axes=ax)
show()

Aggregate beta maps from the LSA model based on condition

Collect the Parameter Estimate maps.

from nilearn.image import concat_imgs

lsa_beta_maps = {cond: [] for cond in events_df["trial_type"].unique()}
trialwise_conditions = lsa_events_df["trial_type"].unique()
for condition in trialwise_conditions:
    beta_map = lsa_glm.compute_contrast(condition, output_type="effect_size")
    # Drop the trial number from the condition name to get the original name,
    # splitting on our delimiter (``__``).
    condition_name = condition.split("__")[0]
    lsa_beta_maps[condition_name].append(beta_map)

# We concatenate the lists of 3D maps into a single 4D beta-series for
# each condition.
lsa_beta_maps = {
    name: concat_imgs(maps) for name, maps in lsa_beta_maps.items()
}

Define the Least Squares-Separate (LSS) models

We next create a separate Least Squares-Separate (LSS) model, one for each trial in the conditions of interest. This is much like the LSA approach, except that we only relabel one trial in the events_df data frame. We loop through the trials, create a version of the data frame where the targeted trial has a unique trial type, fit the model to that data frame, and finally collect the targeted trial’s beta map for the beta-series.

def lss_transformer(events_df, row_number):
    """Label one trial for one LSS model.

    Parameters
    ----------
    df : pandas.DataFrame
        BIDS-compliant events file information.
    row_number : int
        Row number in the data frame.
        This indexes the trial that will be isolated.

    Returns
    -------
    df : pandas.DataFrame
        Update events information,
        with the selected trial's trial type isolated.
    trial_name : :obj:`str`
        Name of the isolated trial's trial type.
    """
    events_df = events_df.copy()

    # Determine which number trial it is *within the condition of interest*.
    trial_condition = events_df.loc[row_number, "trial_type"]
    trial_type_series = events_df["trial_type"]
    trial_type_series = trial_type_series.loc[
        trial_type_series == trial_condition
    ]
    trial_type_list = trial_type_series.index.tolist()
    trial_number = trial_type_list.index(row_number)

    # We again use a unique delimiter here (``__``) that shouldn't be in the
    # original condition names.
    # Technically, all you need is for the requested trial to have a unique
    # 'trial_type' *within* the data frame, rather than across models.
    # However, we may want to have meaningful 'trial_type's (e.g., 'Left_001')
    # across models, so that we could track individual trials across models.
    trial_name = f"{trial_condition}__{trial_number:03d}"
    events_df.loc[row_number, "trial_type"] = trial_name
    return events_df, trial_name


# Loop through the trials of interest and transform the data frame for LSS.
lss_beta_maps = {cond: [] for cond in events_df["trial_type"].unique()}
lss_design_matrices = []

for i_trial in range(events_df.shape[0]):
    lss_events_df, trial_condition = lss_transformer(events_df, i_trial)

    # Compute and collect beta maps.
    lss_glm = FirstLevelModel(**glm_parameters)
    lss_glm.fit(fmri_file, lss_events_df)

    # We save the design matrices across trials to visualize them later.
    lss_design_matrices.append(lss_glm.design_matrices_[0])

    beta_map = lss_glm.compute_contrast(
        trial_condition,
        output_type="effect_size",
    )

    # Drop the trial number from the condition name to get the original name,
    # splitting on our delimiter (``__``)
    condition_name = trial_condition.split("__")[0]
    lss_beta_maps[condition_name].append(beta_map)

# We again concatenate the lists of 3D maps into a single 4D beta-series for
# each condition.
lss_beta_maps = {
    name: concat_imgs(maps) for name, maps in lss_beta_maps.items()
}

Show the LSS design matrices for the first few trials

fig, axes = plt.subplots(ncols=3, figsize=(40, 20))
for i_trial in range(3):
    plot_design_matrix(
        lss_design_matrices[i_trial],
        axes=axes[i_trial],
    )
    axes[i_trial].set_title(f"Trial {i_trial + 1}")

show()

Compare design matrices from the three modeling approaches

DM_TITLES = ["Standard GLM", "LSA Model", "LSS Model (Trial 1)"]
DESIGN_MATRICES = [
    standard_glm.design_matrices_[0],
    lsa_glm.design_matrices_[0],
    lss_design_matrices[0],
]

fig, axes = plt.subplots(
    ncols=3,
    figsize=(40, 20),
    gridspec_kw={"width_ratios": [1, 2, 1]},
)

for i_ax, _ in enumerate(axes):
    plot_design_matrix(DESIGN_MATRICES[i_ax], axes=axes[i_ax])
    axes[i_ax].set_title(DM_TITLES[i_ax])

show()

Applications of beta-series

Beta-series can be used much like resting-state data, though generally with vastly reduced degrees of freedom compared to a typical resting-state run, given that the number of trials should always be less than the number of volumes in a fMRI run.

Two common applications of beta-series are to functional connectivity and decoding analyses. For an example of a beta-series applied to decoding, see Decoding of a dataset after GLM fit for signal extraction. Here, we show how the beta-series can be applied to functional connectivity analysis.

In the following section, we perform a task-based functional connectivity analysis using the two task conditions (“language” and “string”) from the LSS beta-series. The seed coordinate is chosen based on a previous Neurosynth meta-analysis. This section is based on Producing single subject maps of seed-to-voxel correlation, which goes into more detail about seed-to-voxel functional connectivity analyses.

import numpy as np

from nilearn.maskers import NiftiMasker, NiftiSpheresMasker

# Use coordinate taken from Neurosynth's "language" meta-analysis.
coords = [(-54, -42, 3)]

# Initialize masker for the Neurosynth seed.
seed_masker = NiftiSpheresMasker(
    coords,
    radius=8,
    detrend=True,
    standardize="zscore_sample",
    memory="nilearn_cache",
    memory_level=1,
    verbose=1,
)

# Initialize a separate masker for the whole brain.
brain_masker = NiftiMasker(
    smoothing_fwhm=6,
    detrend=True,
    standardize="zscore_sample",
    memory="nilearn_cache",
    memory_level=1,
    verbose=1,
)

Perform the seed-to-voxel correlation for the LSS beta-series.

Using the defined seed_masker, we extract the signal from our Neurosynth coordinate for the LSS beta-series maps defined for each trial, separately for each task type (“language” and “string”). We then extract the whole-brain signal using the brain_masker from the same LSS beta-series maps. Finally, we take the dot-product of these two matrices and normalize it by the number of samples; i.e., the number of trials.

We perform this analysis for each task condition separately.

language_seed_beta_series = seed_masker.fit_transform(
    lss_beta_maps["language"]
)
language_beta_series = brain_masker.fit_transform(lss_beta_maps["language"])
language_corrs = (
    np.dot(
        language_beta_series.T,
        language_seed_beta_series,
    )
    / language_seed_beta_series.shape[0]
)
language_connectivity_img = brain_masker.inverse_transform(language_corrs.T)

# Perform the same seed-to-voxel correlation for the LSS 'string' beta-series
string_seed_beta_series = seed_masker.fit_transform(lss_beta_maps["string"])
string_beta_series = brain_masker.fit_transform(lss_beta_maps["string"])
string_corrs = (
    np.dot(
        string_beta_series.T,
        string_seed_beta_series,
    )
    / string_seed_beta_series.shape[0]
)
string_connectivity_img = brain_masker.inverse_transform(string_corrs.T)

[NiftiSpheresMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7f5fa76b0550>
[NiftiSpheresMasker.wrapped] Finished fit
________________________________________________________________________________
[Memory] Calling nilearn.maskers.base_masker.filter_and_extract...
filter_and_extract(<nibabel.nifti1.Nifti1Image object at 0x7f5fa76b0550>, <nilearn.maskers.nifti_spheres_masker._ExtractionFunctor object at 0x7f5f8866b730>,
{ 'allow_overlap': False,
  'clean_args': None,
  'clean_kwargs': {},
  'detrend': True,
  'dtype': None,
  'high_pass': None,
  'high_variance_confounds': False,
  'low_pass': None,
  'mask_img': None,
  'radius': 8,
  'reports': True,
  'seeds': [(-54, -42, 3)],
  'smoothing_fwhm': None,
  'standardize': 'zscore_sample',
  'standardize_confounds': True,
  't_r': None}, confounds=None, sample_mask=None, dtype=None, sklearn_output_config=None, memory=Memory(location=nilearn_cache/joblib), memory_level=1, verbose=1)
[NiftiSpheresMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7f5fa76b0550>
[NiftiSpheresMasker.wrapped] Extracting region signals
[NiftiSpheresMasker.wrapped] Cleaning extracted signals
_______________________________________________filter_and_extract - 1.5s, 0.0min
[NiftiMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7f5fa76b0550>
[NiftiMasker.wrapped] Computing mask
________________________________________________________________________________
[Memory] Calling nilearn.masking.compute_background_mask...
compute_background_mask(<nibabel.nifti1.Nifti1Image object at 0x7f5fa76b0550>, verbose=0)
__________________________________________compute_background_mask - 0.0s, 0.0min
[NiftiMasker.wrapped] Resampling mask
________________________________________________________________________________
[Memory] Calling nilearn.image.resampling.resample_img...
resample_img(<nibabel.nifti1.Nifti1Image object at 0x7f5fd0778c10>, target_affine=None, target_shape=None, copy=False, interpolation='nearest')
_____________________________________________________resample_img - 0.0s, 0.0min
[NiftiMasker.wrapped] Finished fit
________________________________________________________________________________
[Memory] Calling nilearn.maskers.nifti_masker.filter_and_mask...
filter_and_mask(<nibabel.nifti1.Nifti1Image object at 0x7f5fa76b0550>, <nibabel.nifti1.Nifti1Image object at 0x7f5fd0778c10>, { 'clean_args': None,
  'clean_kwargs': {},
  'cmap': 'gray',
  'detrend': True,
  'dtype': None,
  'high_pass': None,
  'high_variance_confounds': False,
  'low_pass': None,
  'reports': True,
  'runs': None,
  'smoothing_fwhm': 6,
  'standardize': 'zscore_sample',
  'standardize_confounds': True,
  't_r': None,
  'target_affine': None,
  'target_shape': None}, memory_level=1, memory=Memory(location=nilearn_cache/joblib), verbose=1, confounds=None, sample_mask=None, copy=True, sklearn_output_config=None)
[NiftiMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7f5fa76b0550>
[NiftiMasker.wrapped] Smoothing images
[NiftiMasker.wrapped] Extracting region signals
[NiftiMasker.wrapped] Cleaning extracted signals
__________________________________________________filter_and_mask - 0.1s, 0.0min
[NiftiMasker.inverse_transform] Computing image from signals
________________________________________________________________________________
[Memory] Calling nilearn.masking.unmask...
unmask(array([[ 0.232559, ..., -0.331564]], dtype=float32), <nibabel.nifti1.Nifti1Image object at 0x7f5fd0778c10>)
___________________________________________________________unmask - 0.0s, 0.0min
[NiftiSpheresMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7f5fc6efb7f0>
[NiftiSpheresMasker.wrapped] Finished fit
________________________________________________________________________________
[Memory] Calling nilearn.maskers.base_masker.filter_and_extract...
filter_and_extract(<nibabel.nifti1.Nifti1Image object at 0x7f5fc6efb7f0>, <nilearn.maskers.nifti_spheres_masker._ExtractionFunctor object at 0x7f5f8a857ee0>,
{ 'allow_overlap': False,
  'clean_args': None,
  'clean_kwargs': {},
  'detrend': True,
  'dtype': None,
  'high_pass': None,
  'high_variance_confounds': False,
  'low_pass': None,
  'mask_img': None,
  'radius': 8,
  'reports': True,
  'seeds': [(-54, -42, 3)],
  'smoothing_fwhm': None,
  'standardize': 'zscore_sample',
  'standardize_confounds': True,
  't_r': None}, confounds=None, sample_mask=None, dtype=None, sklearn_output_config=None, memory=Memory(location=nilearn_cache/joblib), memory_level=1, verbose=1)
[NiftiSpheresMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7f5fc6efb7f0>
[NiftiSpheresMasker.wrapped] Extracting region signals
[NiftiSpheresMasker.wrapped] Cleaning extracted signals
_______________________________________________filter_and_extract - 1.0s, 0.0min
[NiftiMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7f5fc6efb7f0>
[NiftiMasker.wrapped] Computing mask
________________________________________________________________________________
[Memory] Calling nilearn.masking.compute_background_mask...
compute_background_mask(<nibabel.nifti1.Nifti1Image object at 0x7f5fc6efb7f0>, verbose=0)
__________________________________________compute_background_mask - 0.0s, 0.0min
[NiftiMasker.wrapped] Resampling mask
________________________________________________________________________________
[Memory] Calling nilearn.image.resampling.resample_img...
resample_img(<nibabel.nifti1.Nifti1Image object at 0x7f5fd0778c10>, target_affine=None, target_shape=None, copy=False, interpolation='nearest')
_____________________________________________________resample_img - 0.0s, 0.0min
[NiftiMasker.wrapped] Finished fit
________________________________________________________________________________
[Memory] Calling nilearn.maskers.nifti_masker.filter_and_mask...
filter_and_mask(<nibabel.nifti1.Nifti1Image object at 0x7f5fc6efb7f0>, <nibabel.nifti1.Nifti1Image object at 0x7f5fd0778c10>, { 'clean_args': None,
  'clean_kwargs': {},
  'cmap': 'gray',
  'detrend': True,
  'dtype': None,
  'high_pass': None,
  'high_variance_confounds': False,
  'low_pass': None,
  'reports': True,
  'runs': None,
  'smoothing_fwhm': 6,
  'standardize': 'zscore_sample',
  'standardize_confounds': True,
  't_r': None,
  'target_affine': None,
  'target_shape': None}, memory_level=1, memory=Memory(location=nilearn_cache/joblib), verbose=1, confounds=None, sample_mask=None, copy=True, sklearn_output_config=None)
[NiftiMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7f5fc6efb7f0>
[NiftiMasker.wrapped] Smoothing images
[NiftiMasker.wrapped] Extracting region signals
[NiftiMasker.wrapped] Cleaning extracted signals
__________________________________________________filter_and_mask - 0.1s, 0.0min
[NiftiMasker.inverse_transform] Computing image from signals
________________________________________________________________________________
[Memory] Calling nilearn.masking.unmask...
unmask(array([[ 0.380079, ..., -0.364701]], dtype=float32), <nibabel.nifti1.Nifti1Image object at 0x7f5fd0778c10>)
___________________________________________________________unmask - 0.0s, 0.0min

Visualize both correlation maps.

fig, axes = plt.subplots(figsize=(10, 8), nrows=2)
conn_imgs = [language_connectivity_img, string_connectivity_img]
conn_img_lbls = ["language", "string"]

for img, lbl, ax in zip(conn_imgs, conn_img_lbls, axes, strict=False):
    display = plot_stat_map(
        img,
        threshold=0.5,
        vmax=1,
        cut_coords=coords[0],
        title=lbl,
        figure=fig,
        axes=ax,
    )
    display.add_markers(
        marker_coords=coords,
        marker_color="g",
        marker_size=200,
    )

fig.suptitle(
    "Least Squares-Separate (LSS) Beta-Series Functional Connectivity"
)

show()

Estimated memory usage: 317 MB