Understanding individual differences in cognitive control is a central goal in psychology and neuroscience. Reliably measuring these differences, however, has proven extremely challenging, at least when using standard measures in cognitive neuroscience such as response times or task-based fMRI activity. While prior work has pinpointed the source of the issue—the vast amount of cross-trial variability within these measures—solutions remain elusive. Here, we propose one potential way forward: an analytic framework that combines hierarchical Bayesian modeling with multivariate decoding of trial-level fMRI data. Using this framework and longitudinal data from the Dual Mechanisms of Cognitive Control project, we estimated individuals’ neural responses associated with cognitive control within a color-word Stroop task, then assessed the reliability of these individuals’ responses across a time interval of several months. We show that in many prefrontal and parietal brain regions, test–retest reliability was near maximal, and that only hierarchical models were able to reveal this state of affairs. Further, when compared to traditional univariate contrasts, multivariate decoding enabled individual-level correlations to be estimated with significantly greater precision. We specifically link these improvements in precision to the optimized suppression of cross-trial variability in decoding. Together, these findings not only indicate that cognitive control-related neural responses individuate people in a highly stable manner across time, but also suggest that integrating hierarchical and multivariate models provides a powerful approach for investigating individual differences in cognitive control, one that can effectively address the issue of high-variability measures.

Analyses of functional connectivity (FC) in resting-state brain networks (RSNs) have generated many insights into cognition. However, the mechanistic underpinnings of FC and RSNs are still not well-understood. It remains debated whether resting-state activity is best characterized as reflecting noise-driven fluctuations around a single stable state, or instead, as a nonlinear dynamical system with nontrivial attractors embedded in the RSNs. Here, we provide evidence for the latter, by constructing whole-brain dynamical systems models from individual resting-state fMRI (rfMRI) recordings, using the Mesoscale Individualized NeuroDynamic (MINDy) framework. The MINDy models consist of hundreds of neural masses representing brain parcels, connected by fully trainable, individualized weights. We found that our models manifested a diverse taxonomy of nontrivial attractor landscapes including multiple equilibria and limit cycles. However, when projected into anatomical space, these attractors mapped onto a limited set of canonical RSNs, including the default mode network (DMN) and frontoparietal control network (FPN), which were reliable at the individual level. Further, by creating convex combinations of models, bifurcations were induced that recapitulated the full spectrum of dynamics found via fitting. These findings suggest that the resting brain traverses a diverse set of dynamics, which generate several distinct but anatomically overlapping attractor landscapes. Treating rfMRI as a unimodal stationary process (i.e., conventional FC) may miss critical attractor properties and structure within the resting brain. Instead, these may be better captured through neural dynamical modeling and analytic approaches. The results provide new insights into the generative mechanisms and intrinsic spatiotemporal organization of brain networks.