State-Space Modeling for fMRI Data

Reference Files

Read these before generating code or providing detailed guidance:

Quick-Reference Decision Tree

Step 1: What is the paradigm?

Resting state → States represent intrinsic brain dynamics. No external timing. HRF is a nuisance — you are modeling BOLD-level dynamics unless you explicitly deconvolve. Use HMM or HMM-MAR on parcellated or ICA data. → Read references/paradigm_guide.md § Resting State

Task-based (event/block design) → External task structure provides ground truth for validation. Key question: are you modeling task-evoked states or residual dynamics beyond the task? HRF alignment is critical — task events are neural-time but BOLD is delayed 4-6s. → Read references/paradigm_guide.md § Task-Based

Naturalistic (movie, TV, gaming) → Continuous stimulation without discrete trial structure. States reflect stimulus-driven + endogenous dynamics. Shared stimulus enables inter-subject state alignment. Typically longer runs, which helps estimation. → Read references/paradigm_guide.md § Naturalistic

Step 2: What data format?

→ See references/preprocessing.md § Dimensionality Reduction for CIFTI and voxel data.

Step 3: Which model family?

"I want discrete brain states and their transitions" → HMM family. Gaussian HMM if states differ in mean/FC patterns. HMM-MAR if states differ in temporal dynamics (autoregressive structure).

"I want continuous latent dynamics that switch between regimes" → SLDS or rSLDS. States govern a linear dynamical system; latent space is continuous.

"I want external inputs (task events, stimuli) to drive state transitions" → Input-output HMM or input-driven SLDS.

"I want hierarchical structure (fast states within slow states)" → Hierarchical HMM or hierarchical SLDS. Useful for naturalistic paradigms.

"Deep data on few subjects" → Per-subject models with more complex structure (rSLDS, HMM-MAR with many lags). Cross-validate within-subject using held-out runs or time segments.

"Many subjects with short scans" → Group-level HMM or concatenation approaches. Simpler models (Gaussian HMM) are more robust. → See references/group_inference.md

→ Full model specs, equations, and breakdown conditions: references/model_catalog.md

Step 4: HRF strategy

The BOLD signal is not a direct readout of neural activity — it's the neural signal convolved with the HRF, introducing a ~4-6s delay and temporal smoothing.

Why this matters: If neural states switch rapidly (<5s), the HRF blurs transitions in BOLD. A 2-second neural state may appear as a gradual 8-10s transition.

→ Full treatment with code: references/hrf_modeling.md

Step 5: Key parameter decisions

Number of states (K): No single right answer. Use:

• BIC (tends to prefer fewer states with fMRI)
• Cross-validated log-likelihood (held-out time segments or runs)
• Free energy (variational methods)
• Stability: fit multiple K values, check state reproducibility across initializations
• Sanity check: states should have interpretable spatial patterns and plausible dwell times

Covariance structure (Gaussian emissions):

• Full: captures FC per state, many parameters — needs substantial data or regularization
• Diagonal: fewer parameters, more robust with limited data
• Factor analysis / low-rank: good for high-dimensional data

Autoregressive order (HMM-MAR): Typical range 1–5 for TR=0.7–2s. Cross-validate.

Initialization: K-means or GMM (much better than random). 20–50 random restarts.

Stickiness: If model infers rapid switching (every TR), check for motion artifacts or add a sticky prior. Sticky HMM adds self-transition bias, discouraging unrealistic rapid switching.

Step 6: What to do after fitting

Connect states to behavior or cognition: → Read references/downstream_analysis.md § Behavioral Correlates

Write a methods/results section or prepare figures: → Read references/downstream_analysis.md § Reporting Checklist and § Required Figures

Sanity check — does the pipeline even work? → Run simulate_and_recover() from references/downstream_analysis.md § Simulation Testing before trusting real-data results.

Covariates and Confounds

Regress out before SSM fitting:

• Head motion parameters (6 or 24-parameter expansion)
• CSF and WM signals (or aCompCor components)
• Use fMRIPrep/XCP-D recommended confound strategy
• Global signal regression: controversial — can introduce anticorrelations

Covariates that can enter the SSM:

• Task timing (onsets, durations) → transition model or emission model
• Stimulus features (naturalistic) → emission model covariates
• Subject-level covariates (age, group) → hierarchical prior

→ Full confound strategy: references/preprocessing.md

Examples

Example 1: Resting-state HMM with hmmlearn

User: "I have fMRIPrep-preprocessed resting-state data parcellated to Schaefer-200. I want to find recurring brain states."

Steps:

1. Read references/paradigm_guide.md § Resting State and references/code_templates.md
2. Recommend Gaussian HMM (states differ in FC patterns, not dynamics)
3. Suggest K=6-12 with BIC/cross-validation
4. Sticky HMM to avoid TR-level switching
5. Use full covariance if data permits (>300 TRs per state), else diagonal

Key code pattern (from references/code_templates.md):

from hmmlearn import hmm

model = hmm.GaussianHMM(
    n_components=K,
    covariance_type="full",  # or "diag" for limited data
    n_iter=100,
    random_state=42
)
# Use K-means init; multiple restarts; pass lengths array for multi-run
model.fit(X, lengths=run_lengths)

Example 2: Task-based SSM with HRF-aware IO-HMM (N-back paradigm)

User: "I have N-back fMRI (0-back, 2-back). I want to model how cognitive load shifts brain states over time."

Steps:

1. Read references/paradigm_guide.md § Task-Based and references/hrf_modeling.md
2. Recommend IO-HMM — external task input drives state transitions
3. HRF strategy: convolve task regressors with canonical HRF before using as inputs
4. Read references/code_templates.md for IO-HMM template

Key decisions:

• Task input: HRF-convolved indicator for 0-back vs. 2-back condition
• Model: ssm.HMM with InputDrivenTransitions
• Covariates: regress out motion + CSF/WM before fitting

Example 3: Naturalistic fMRI with rSLDS (movie-watching)

User: "I have movie-watching fMRI from 50 subjects. I want continuous latent dynamics that switch between regimes, not just discrete states."

Steps:

1. Read references/paradigm_guide.md § Naturalistic and references/model_catalog.md § rSLDS
2. Recommend rSLDS — latent space with state-dependent switching boundaries
3. Start with SLDS first; upgrade if SLDS fits poorly
4. Group-level inference: references/group_inference.md

Key considerations:

• Fit on BOLD directly (naturalistic = slow states, HRF less critical)
• Use ssm.SLDS(recurrent=True) with Laplace-EM
• Shared stimulus across subjects enables inter-subject state alignment

Common Pitfalls

1. Fitting SSMs to raw data. Always preprocess fully first (fMRIPrep + XCP-D with appropriate denoising).

2. Ignoring HRF when interpreting fast states. States lasting 2-3 TRs in BOLD ≠ 2-3 second neural events. Acknowledge the BOLD timescale.

3. Too many states for the data. Rule of thumb: ≥50-100 time points per state for stable full-covariance Gaussian HMM. A 5-min scan (~300 TRs) cannot reliably support a 12-state model with 100 ROIs.

4. Naive multi-run concatenation. The state at end of run 1 has no temporal relation to start of run 2. Pass lengths array to reset the forward algorithm at run boundaries.

5. State label switching across subjects. HMMs are invariant to state relabeling. Use Hungarian algorithm on emission parameters or fit a group model to align labels.

6. Motion-driven states. High-motion time points create artifactual states. Scrub/censor before fitting; verify states don't correlate with framewise displacement.

7. Circular analysis. Don't select K on full data then refit and report. Use proper cross-validation or a held-out set.

Python Library Landscape

→ Working code examples for each library: references/code_templates.md

Language & Compute Environment

This skill is Python-first. All code targets Python ≥ 3.10 with the library versions in each template. The ssm library (Linderman lab) provides a JAX backend that runs on GPU transparently — the same API works on CPU or GPU.

GPU acceleration — recommended, not required:

All code in references/code_templates.md runs on CPU. GPU is a drop-in speedup. See references/code_templates.md §12 (JAX/GPU) and §13 (dynamax Lego models).

JAX note: Install JAX before importing ssm for GPU support:

# CPU-only JAX (default)
pip install jax jaxlib ssm

# GPU (CUDA 12)
pip install "jax[cuda12]" -f https://storage.googleapis.com/jax-releases/jax_cuda_releases.html
pip install ssm

Julia: For Bayesian HMM inference with full MCMC uncertainty, Turing.jl and StateSpaceModels.jl are excellent alternatives with richer posterior sampling than Python's pyhsmm. Recommend if the user is already in the Julia ecosystem or needs calibrated credible intervals on transition probabilities.