6. Functional Connectivity and Resting State

Why functional connectivity?

A task-based GLM tells you where the brain responds when a stimulus appears. It is silent about how regions talk to each other — and silent about brains in the absence of any task at all. Functional connectivity (FC) plugs that gap. The starting observation, due to Biswal et al. (1995), is that even with a subject lying still in the scanner, the BOLD signal in distant grey-matter voxels fluctuates coherently in the low-frequency band (~0.01–0.1 Hz). Sensorimotor cortices on the two hemispheres rise and fall together; the posterior cingulate hums in step with the medial prefrontal cortex. These spontaneous co-fluctuations turn out to organise themselves into the same large-scale networks that show up during tasks — the default-mode network being the canonical example.

A functional connectome is the resulting graph: nodes are brain regions, edges are the statistical dependencies between their time courses. Resting-state fMRI is the most common substrate, but the same machinery applies to task fMRI (cross-link 7. Task-fMRI Connectivity (PPI and RSA)).

what FC is and is not

FC is a statistical dependency between time series; it is not causal and not directional. “Functional” here just contrasts with “structural” (DTI tractography). Inferring direction needs dynamic models (DCM, Granger) — out of scope for this chapter.

---

0 Two ways to build an FC map

Before the machinery, the big-picture choice. There are two complementary ways to turn a resting-state run into connectivity, plus a data-driven third:

• Seed-based (hypothesis-driven): pick one region of interest — the seed — extract its mean time course, and correlate it against every other voxel in the brain. The result is a single whole-brain FC map: “what is coupled to this region?” Simple, interpretable, and the classic way the default-mode network was characterised.
• Whole-connectome (parcellation-based): reduce the brain to N nodes and correlate all pairs, giving an N × N connectivity matrix — “how is everything coupled to everything?” This is the bulk of this chapter (§1–§6).
• Data-driven (ICA/dictionary learning): let the data decompose itself into networks with no seed and no atlas (§7).

Seed-based analysis is the smallest, most transparent unit and the conceptual bridge to task connectivity — a PPI (7. Task-fMRI Connectivity (PPI and RSA)) is a seed-based analysis with a task-interaction term added. Here is the whole method:

from nilearn.maskers import NiftiSpheresMasker, NiftiMasker
from nilearn.interfaces.fmriprep import load_confounds
import numpy as np
 
# 1. Seed = a 6 mm sphere at the PCC (a default-mode hub), MNI coords
pcc = [(0, -52, 26)]
confounds, sample_mask = load_confounds(
    bold, strategy=("motion", "high_pass", "wm_csf", "global_signal"))
 
seed_masker = NiftiSpheresMasker(pcc, radius=6, standardize="zscore_sample",
                                 t_r=t_r, low_pass=0.1, high_pass=0.01)
seed_ts = seed_masker.fit_transform(bold, confounds=confounds,
                                    sample_mask=sample_mask)          # (T, 1)
 
# 2. Every voxel's time course, cleaned the same way
brain_masker = NiftiMasker(standardize="zscore_sample", t_r=t_r,
                           low_pass=0.1, high_pass=0.01)
brain_ts = brain_masker.fit_transform(bold, confounds=confounds,
                                      sample_mask=sample_mask)        # (T, n_voxels)
 
# 3. Correlate seed with every voxel -> a seed FC map
seed_to_vox = (brain_ts.T @ seed_ts).ravel() / seed_ts.shape[0]      # Pearson r per voxel
fc_map = brain_masker.inverse_transform(np.arctanh(seed_to_vox))     # Fisher-z, back to a NIfTI

Plot fc_map and you get the familiar default-mode network. The rest of this chapter scales this idea up from one seed to an all-pairs connectome.

---

1 From voxels to nodes

A connectome lives or dies by its parcellation. You need to reduce ~200 000 grey-matter voxels to a tractable set of nodes (typically 100–500) before correlating anything. Three families of nodes:

1.1 Deterministic atlases

Each voxel belongs to exactly one labelled region. These are the workhorses for a first project because the regions have names you can put in a paper.

Atlas	Regions	Native MNI space	Nilearn fetcher
Harvard-Oxford (cort + sub)	48 + 21	MNI152NLin6Asym (FSL)	fetch_atlas_harvard_oxford
AAL	116	MNI152	fetch_atlas_aal
Schaefer 2018	100 / 200 / 400 / 1000	MNI152NLin6Asym	fetch_atlas_schaefer_2018
Yeo 7 / 17 networks	7 / 17	MNI152NLin6Asym	fetch_atlas_yeo_2011

match the MNI template your data lives in

fMRIPrep’s default normalisation target is MNI152NLin2009cAsym, while most older atlases (Harvard-Oxford, Yeo, Schaefer) ship in MNI152NLin6Asym. The two templates differ by ~3 mm at the cortical edge — enough to mislabel insular voxels as STG. Either request the matching space in fMRIPrep (--output-spaces MNI152NLin6Asym:res-2, see 2.3 fmriPrep with Docker) or resample the atlas with image.resample_to_img(atlas, ref) before use.

1.2 Probabilistic / coordinate atlases

• MSDL (39 regions, Varoquaux 2011) — probabilistic, optimised for resting-state. Use with NiftiMapsMasker.
• Power 264 (Power et al. 2011) — 264 spherical 10 mm ROIs at MNI coordinates. Pre-assigned to 13 networks; standard for graph-theoretic work.
• DiFuMo — multi-scale (64–1024 components) probabilistic dictionary.

1.3 Data-driven parcels

When no atlas matches your hypothesis you can let the data parcellate itself: ICA (CanICA), dictionary learning (DictLearning), or Ward/K-means clustering (nilearn.regions.Parcellations). The trade-off is interpretability: data-driven components do not come with anatomical labels.

sensible defaults

For a master’s project, start with Schaefer-400 if you care about cortex only, or Harvard-Oxford + Harvard-Oxford-sub if you want subcortical nodes too. Both give you labelled regions, both run in seconds.

---

2 Extracting time series

Once nodes are defined, average the preprocessed BOLD within each region. Nilearn’s maskers do this in one call, and crucially they let you regress out confounds on the same data scaling, which is what you want for FC.

import os
from nilearn import datasets
from nilearn.maskers import NiftiLabelsMasker
from nilearn.interfaces.fmriprep import load_confounds
 
bids_dir     = os.environ["BIDS_DIR"]
fmriprep_dir = os.environ["FMRIPREP_DIR"]
sub, task    = "01", "rest"
 
bold = (f"{fmriprep_dir}/sub-{sub}/func/"
        f"sub-{sub}_task-{task}_space-MNI152NLin6Asym_desc-preproc_bold.nii.gz")
 
# Schaefer-400, 7-network version, in MNI152NLin6Asym at 2 mm
atlas = datasets.fetch_atlas_schaefer_2018(n_rois=400, yeo_networks=7, resolution_mm=2)
 
# Read TR from the BIDS sidecar — never hard-code it
import json
with open(f"{os.environ['BIDS_DIR']}/sub-{sub}/func/"
          f"sub-{sub}_task-{task}_bold.json") as f:
    t_r = json.load(f)["RepetitionTime"]
 
# Confounds: the Parkes 2018 sweet spot for rest-state
confounds, sample_mask = load_confounds(
    bold,
    strategy=("motion", "wm_csf", "global_signal", "scrub", "high_pass"),
    motion="full",               # 24 HMP: 6 + derivatives + squares
    wm_csf="basic",
    global_signal="basic",       # GSR — controversial, see below
    scrub=5, fd_threshold=0.5, std_dvars_threshold=1.5,
)
 
masker = NiftiLabelsMasker(
    labels_img=atlas.maps,
    labels=atlas.labels,
    standardize="zscore_sample",
    detrend=True,
    low_pass=0.10, high_pass=0.01, t_r=t_r,
    memory="nilearn_cache", verbose=0,
)
 
ts = masker.fit_transform(bold, confounds=confounds, sample_mask=sample_mask)
# ts shape: (n_timepoints_kept, 400)

A few things are doing real work here:

• standardize="zscore_sample" z-scores each region’s time course. Required for the Pearson correlation that follows to be on a comparable scale.
• detrend=True removes linear drift (low-pass scanner drift the high-pass filter does not catch).
• low_pass=0.10, high_pass=0.01 restrict the analysis to the canonical rest-state band. Apply filtering at the same step as confound regression, otherwise the filter reintroduces frequencies the regression had killed (Lindquist’s “don’t filter twice” rule).
• confounds=... regresses motion, white-matter, CSF, and (optionally) global signal out of the regions’ time series. This is the most consequential choice in the pipeline — see 5.2 Confound Selection for the rationale, the 24HMP + aCompCor + GSR + scrubbing strategy, and the per-subject mean-FD screening code.
• sample_mask drops the time points flagged for motion (scrub), instead of zero-ing them. Important when you compute correlations on the remaining points.

GSR is a choice, not a default

Global-signal regression sharpens network structure and reduces motion artefact, but it forces the connectome to have negative correlations somewhere by construction (the rows of the design must sum to zero). Parkes et al. (2018) found GSR helps on motion benchmarks; many groups still avoid it. Decide before you look at the data, and report it.

---

3 Computing the connectome

Once you have a (T, N) time-series matrix, the connectome is just a similarity matrix between columns.

from nilearn.connectome import ConnectivityMeasure
 
cm = ConnectivityMeasure(kind="correlation", standardize="zscore_sample")
conn = cm.fit_transform([ts])[0]      # shape (400, 400)

kind accepts:

• "correlation" — Pearson r. The default; what most papers report.
• "partial correlation" — r between two regions after regressing out all other regions. Estimates direct (rather than mediated) coupling. Sensitive to estimation noise when N is large.
• "covariance" / "precision" — the matrix-valued versions.
• "tangent" — see §5.

Fisher z before you average

Pearson r is bounded in [-1, 1] and skewed near the tails. Apply Fisher’s z transform $z=\frac{1}{2}\ln \left(\frac{1+r}{1-r}\right)=\mathrm{arctanh}(r)$ before averaging across subjects or running parametric stats; invert with np.tanh for display.

---

4 Visualising the connectome

import numpy as np
from nilearn import plotting
from nilearn.plotting import find_parcellation_cut_coords
 
coords = find_parcellation_cut_coords(labels_img=atlas.maps)
 
# 4a. heatmap
plotting.plot_matrix(conn, vmin=-0.8, vmax=0.8,
                     labels=[l.decode() for l in atlas.labels],
                     reorder=True, figure=(10, 10))
 
# 4b. 3D glass brain, showing only the strongest 20% of edges by |weight|
plotting.plot_connectome(conn, coords, edge_threshold="80%",
                         node_size=20, title="Schaefer-400 connectome")
 
# 4c. interactive HTML — drop into a Quartz page or open in browser
view = plotting.view_connectome(conn, coords, edge_threshold="95%")
view.save_as_html("connectome.html")

reorder=True runs hierarchical clustering on the rows so block structure (i.e. networks) becomes visible. For probabilistic atlases use find_probabilistic_atlas_cut_coords instead.

---

5 Group-level connectomes

The naive approach: compute each subject’s correlation matrix, Fisher-z it, average. Test edges with one-sample t-tests across subjects.

subjects = ["01", "02", "03", "04", "05"]
mats = []
for s in subjects:
    ts_s = masker.fit_transform(bold_for(s), confounds=conf_for(s))
    mats.append(ConnectivityMeasure(kind="correlation").fit_transform([ts_s])[0])
 
z_mats = np.arctanh(np.clip(np.stack(mats), -0.999, 0.999))
group  = np.tanh(z_mats.mean(axis=0))

This works, but each subject’s covariance is noisy and you lose information about between-subject geometry. The tangent-space projection (Varoquaux et al. 2010) fixes both: it represents each subject as a deviation from a group mean covariance computed on the Riemannian manifold of positive-definite matrices. Tangent-space features routinely outperform raw correlations as inputs to classifiers and group statistics.

cm_tan = ConnectivityMeasure(kind="tangent")
tan_mats = cm_tan.fit_transform(list_of_ts)   # list of (T_i, N) arrays
group_ref = cm_tan.mean_                       # group mean covariance

Use tangent when you plan to do classification or group comparison; use correlation (+ Fisher z) when you want per-edge p-values you can map back to anatomy.

---

6 Resting-state networks (RSNs)

The canonical Yeo et al. (2011) parcellation groups the cortex into seven networks: Visual, Somatomotor, Dorsal Attention, Ventral Attention, Limbic, Frontoparietal, and Default Mode. The Schaefer atlas you fetched above is pre-assigned to these networks via its label strings.

# Schaefer labels look like b'7Networks_LH_Vis_1', b'7Networks_RH_Default_PCC_1', ...
nets = np.array([l.decode().split("_")[2] for l in atlas.labels])
 
# Within- and between-network FC
from itertools import product
unique_nets = np.unique(nets)
M = np.zeros((len(unique_nets), len(unique_nets)))
for i, a in enumerate(unique_nets):
    for j, b in enumerate(unique_nets):
        mask = np.outer(nets == a, nets == b)
        if a == b:
            mask &= ~np.eye(len(nets), dtype=bool)   # drop self-edges
        M[i, j] = np.arctanh(conn[mask]).mean()
M = np.tanh(M)

M is a 7 × 7 matrix whose diagonal is within-network coherence and whose off-diagonals are between-network coupling — a much more compact summary than 400 × 400, and the right granularity for behavioural correlations in small samples.

---

7 Data-driven networks: ICA and dictionary learning

When you have no atlas or want to discover networks from the data themselves:

from nilearn.decomposition import CanICA, DictLearning
 
ica = CanICA(n_components=20, mask_strategy="whole-brain-template",
             smoothing_fwhm=6.0, threshold=3.0,
             standardize="zscore_sample", random_state=0,
             memory="nilearn_cache").fit(list_of_bolds)
components_img = ica.components_img_     # 4D image, one volume per component
 
dl = DictLearning(n_components=20, smoothing_fwhm=6.0, alpha=10,
                  standardize="zscore_sample", random_state=0).fit(list_of_bolds)

CanICA recovers spatially independent networks (visual, default mode, salience, …); DictLearning produces sparser, often more interpretable maps. Either output can be fed into RegionExtractor to split each component into separate ROIs, then back into NiftiMapsMasker to compute a fully data-driven connectome.

picking k

20 components is a good starting point for whole-brain rest. 40–70 captures finer sub-networks; above 100 you start splitting noise.

---

8 Common pitfalls

checklist before you publish

1. Confound strategy applied only once. Do not regress confounds in NiftiLabelsMasker and in a separate step — you will partial out the same variance twice and inflate FC.
2. Always Fisher-z Pearson r before averaging or t-testing.
3. Motion-screen subjects — exclude any with mean FD > 0.2 mm or >20% scrubbed volumes (see the screening snippet in 5.2 Confound Selection). Without this, group FC differences become group motion differences.
4. Match MNI templates. Atlas and BOLD must live in the same space. Resample one or rerun fMRIPrep with the right --output-spaces.
5. Report the band-pass. 0.01–0.10 Hz is standard; 0.008–0.09 Hz is common; report what you used.
6. Mind the TR. With multi-band TR ≈ 0.8 s you have ~6× the samples per minute — both the filter cut-offs and the scrubbing thresholds should be re-examined.

---

9 Where to go next

• Task-based FC, PPI, RSA → 7. Task-fMRI Connectivity (PPI and RSA).
• Graph-theoretic measures (clustering coefficient, modularity, rich-club) — nilearn does not ship these. Use bctpy (Python port of the Brain Connectivity Toolbox) or networkx on the thresholded matrix.
• Dynamic FC — sliding-window correlations, HMMs (see hmmlearn, pydfc).
• Machine learning on connectomes — feed cm_tan.fit_transform(...) straight into sklearn; predict age, diagnosis, behaviour. Always cross-validate at the subject level.

---

References

• Biswal, B., Yetkin, F. Z., Haughton, V. M., & Hyde, J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magnetic Resonance in Medicine, 34(4), 537–541. doi:10.1002/mrm.1910340409
• Yeo, B. T. T., Krienen, F. M., Sepulcre, J., Sabuncu, M. R., Lashkari, D., Hollinshead, M., … Buckner, R. L. (2011). The organization of the human cerebral cortex estimated by intrinsic functional connectivity. Journal of Neurophysiology, 106(3), 1125–1165. doi:10.1152/jn.00338.2011
• Schaefer, A., Kong, R., Gordon, E. M., Laumann, T. O., Zuo, X.-N., Holmes, A. J., Eickhoff, S. B., & Yeo, B. T. T. (2018). Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI. Cerebral Cortex, 28(9), 3095–3114. doi:10.1093/cercor/bhx179
• Power, J. D., Cohen, A. L., Nelson, S. M., Wig, G. S., Barnes, K. A., Church, J. A., … Petersen, S. E. (2011). Functional network organization of the human brain. Neuron, 72(4), 665–678. doi:10.1016/j.neuron.2011.09.006
• Varoquaux, G., Baronnet, F., Kleinschmidt, A., Fillard, P., & Thirion, B. (2010). Detection of brain functional-connectivity difference in post-stroke patients using group-level covariance modeling. MICCAI 2010, LNCS 6361, 200–208. doi:10.1007/978-3-642-15705-9_25
• Parkes, L., Fulcher, B., Yücel, M., & Fornito, A. (2018). An evaluation of the efficacy, reliability, and sensitivity of motion correction strategies for resting-state functional MRI. NeuroImage, 171, 415–436. doi:10.1016/j.neuroimage.2017.12.073