Geometry of neural dynamics along the cortical attractor landscape reflects changes in attention

Domain: neuroscience
Published: 2026-02-12
Editorial outcome: minor_revision

Abstract
Large-scale brain activity reflects changes in attention. To understand how, we tested a hypothesis that the geometry of neural dynamics on the cortical attractor landscape, or movement along its “hills and valleys”, reflects attentional states. A dynamical systems model separating intrinsic dynamics from stimulus-driven influences was fit to fMRI data collected during rest, tasks, and movie-watching. Model simulations revealed a set of attractors aligned with canonical functional brain networks. The speed and direction of neural trajectories toward these attractors varied systematically with attentional states over time and across contexts. When participants were paying attention to effortful tasks, neural dynamics converged fast and directly toward a task-relevant attractor. In contrast, when participants were engaged in sitcom episodes, neural dynamics occupied a flatter region of the landscape, directed away from attractors. These findings demonstrate that while attractor locations are largely determined by cortical organization, the geometry of neural dynamics changes systematically across attentional states and contexts.

Introduction
The brain is a multiscale system that dynamically evolves over time to produce behavior. Dynamical systems modeling, which formalizes how neural activity changes over time using differential equations, has played a central role in neuroscience1,2. Such models have uncovered numerous biophysical mechanisms of the brain, from the early work of Hodgkin and Huxley3 who described how ion currents generate action potentials, to Wilson and Cowan4 who described how excitatory-inhibitory neuronal interactions shape population-level dynamics. More recent work has extended this approach beyond single cells and neuronal populations to characterize whole-brain dynamics5–10.

From a dynamical systems perspective, large-scale brain activity unfolds over time as a trajectory within a high-dimensional state space, where each dimension typically represents the activity of a brain region. This state space is shaped by an attractor landscape, or hills and valleys that guide the trajectory of neural activity (Fig. 1). Research has shown that trajectories visit a small number of recurring brain states, each defined by a unique pattern of regional activity and interactions11–15. These brain states function as attractors, or deep valleys in the landscape, where neural activity is most likely to converge16–19. In turn, neural dynamics are characterized by transitions between or perturbations around these attractors.

While attractor landscapes govern neural dynamics, it is important to understand how they relate to our internal states, including physiological (e.g., hunger, stress, arousal) and cognitive states (e.g., attention, emotion, motivation)20. Prior studies have provided initial evidence of this connection, showing that certain brain states are more or less likely to occur when a person is focused versus unfocused21, performing better or worse on a task22,23, or engaged versus disengaged in a movie24. However, the computational mechanism of this connection remains unknown25. One way to characterize such mechanisms is through the geometry of neural dynamics—the structure and flow of trajectories that describe how brain activity evolves over time along the attractor landscape (Fig. 1). This geometry can be quantified by the speed and direction with which brain activity moves toward or away from attractors. We hypothesized that the geometry of neural dynamics changes over time in relation to internal state fluctuations. Supporting this idea, Munn et al.26 showed evidence of flattened attractor landscape during phasic bursts in noradrenergic locus coeruleus activity (related to arousal)27 and steepened attractor landscape during phasic bursts in cholinergic basal forebrain activity (related to vigilance or attentional focus)28. This suggests a possibility that one’s internal state at a given moment may correspond to whether the brain occupies a flatter or steeper region of the attractor landscape, which in turn shapes the trajectory of neural dynamics.

Here, we propose that the geometry of neural dynamics on the attractor landscape characterizes moment-to-moment and context-to-context variations in internal states. In this study, we specifically test this in relation to measures of sustained attention. Dynamical systems models were fit to whole-brain fMRI data collected during rest, tasks, and movie-watching. Our model separates neural activity into two components: one that is intrinsic, driven internally by regional activity and interactions, and the other extrinsic, driven externally from the stimulus. By simulating neural trajectories from the model, we identify attractors toward which the neural activity converges. At each moment, we estimate the speed and direction of these simulated trajectories, using them to infer the steepness of the local attractor landscape. Importantly, we relate the change in geometries to behavioral measures of attention which were collected with the fMRI data. These measures include participants’ continuous ratings of engagement while watching comedy sitcoms and button response times during controlled sustained attention tasks.

Prior work has typically fit a set of static model parameters from neural activity dynamics, meaning variations across time have often been reduced to a single model that is agnostic to temporal change. Here, we reproduce evolving neural geometry from a model of fixed parameters and relate these dynamics to fluctuations in internal states. With this approach, we test a hypothesis that the geometry of large-scale cortical dynamics along the attractor landscape systematically reflects changes in attentional states during task and movie-watching contexts.

Results
A dynamical systems model of large-scale cortical activity
We applied a large-scale parametric dynamical systems model, developed and validated by Singh et al.9 and Chen et al.19, to fit the time series of BOLD activity measured in human cortex with fMRI. The model defines the rate of change in neural activity x_ t, as a function of the neural activity x t (x t 2Rn, n= number of neural units) and time-aligned experimental variables u t (u t 2Rp, p= number of experimental variables), such as task, stimulus, or behavior. In our model, X corresponds to the BOLD activity timeseries of 200 parcels covering the cortex29. U corresponds to audiovisual and semantic features extracted from the stimuli that participants watched and heard inside the scanner, which were reduced to 100 principal component dimensions. Mathematically, the model can be described as x_ t=F(x t)+G(u t), where F and G are deterministic functions that represent, respectively, the impact of intrinsic activity and input-driven perturbations on the evolution of neural activity.

The followings are the specifics of our model, where F(x t)=W\psi \alpha (x t)-D\odot x t and G(u t)=\beta u t with W, D, \alpha, and \beta being model parameters (Fig. 2).

x_ t = (x(t+\Delta t)-x(t))/\Delta t = W\psi \alpha (x t)-D\odot x t + \beta u t (1)

\psi \alpha (x t) = \sqrt{\alpha^2+(bx t+0.5)^2} - \sqrt{\alpha^2+(bx t-0.5)^2} (2)

Given that \Delta t equals to 1 TR in our data, the equation can be simplified as follows.

x^ t+1 = x t + W\psi \alpha (x t)-D\odot x t + \beta u t (3)

The goal of the model is to predict neural activity pattern of the consecutive timestep (x^ t+1), by minimizing the prediction error (i.e., difference between the predicted x^ t+1 and observed x t+1) using algorithmic optimization (here, stochastic gradient descent). The prediction of the next timestep is based on the neural activity x t plus the signals received by connections from other parcels (W\psi \alpha (x t)) minus the self-decay (D\odot x t) plus the neural activity driven by the external inputs (\beta u t). In turn, fitting this model corresponds to decomposing intrinsic and extrinsic (i.e., input-driven) neural dynamics.

The weight matrix (W 2Rn×n) represents directional interaction between neural units, namely the effective connectivity. Nonlinearity is introduced by a parametrized sigmoidal transfer function (\psi a), which maps neural activity to a bounded output at a range from −1 to 1. The slope of the transfer function differs for every parcel, parameterized by \alpha (\alpha 2Rn). b is fixed at 20/3 following prior studies9,19. A self-decay (D 2Rn) captures a return to baseline in absence of external inputs or interactions, with \odot denoting element-wise product. Having high decay rate corresponds to having low temporal autocorrelation, indicating less persistence in neural activity. \beta represents linear transformation from the stimulus embedding space to neural activity space. Note that this is an individualized model where a set of parameters are estimated for each fMRI run of each participant.

We analyzed two openly available fMRI datasets, the SONG dataset24 (N=27) and the HCP dataset30–32 (N=119). In both datasets, each participant underwent multiple sessions of rest, task, and movie-watching conditions. We extracted timeseries of stimulus features for the movie-watching runs, including low-level visual and visuo-semantic features of video frames, low-level audio of the sounds, and audio-semantic features of the speech and dialogue in the movies. These features were projected onto 100 principal components that explained the largest variance across runs. Only low-level visual and visuo-semantic features were extracted from the SONG dataset task runs because the stimuli were purely visual. Given that no stimulus was provided for resting-state runs in either dataset, we fit a simplified model that only considers intrinsic but not extrinsic neural dynamics: x_ t=W\psi \alpha (x t)-D\odot x t (removing the input term \beta u t from the equation).

Model parameters recapitulate functional brain connectivity and stimulus encoding
We first tested whether our model, optimized to predict neural activity of successive timesteps, captured parameters sensitive to individual and cognitive state differences. These properties are critical, as their emergence indicates the model’s biological plausibility. We analyzed the movie-watching runs of the SONG and HCP datasets because they were fit on the full model that includes the input term (see Supplementary Fig. S1 for the range of model parameter estimates).

Model performance was estimated based on how well the model predicted neural activity of the next time steps, specifically r2(x^ 2:T, x 2:T). The total explained variance was, on average, r2=0.452±0.043 across a total of 80 runs in the SONG dataset and r2=0.561±0.053 across 476 runs in the HCP dataset (Fig. 3A). This indicates that our model explained approximately half of the variance in neural activity. We then decomposed the explained variance into three components. The interareal interaction explained r2=0.041±0.009 (SONG), 0.055±0.015 (HCP), local recurrence explained r2=0.284±0.051, 0.237±0.060, and external input-driven activity explained r2=0.017±0.003, 0.015±0.004 of the variances.

We validated model performance across runs. Model parameters estimated from one run of a participant were better at predicting that same participant’s brain activity in other runs than predicting other participants’ brain activity (SONG: within: r2=0.418±0.037, across: r2=0.402±0.011; paired t-test t(79)=4.049, p=0.0001; HCP: within: r2=0.526±0.052, across: r2=0.508±0.015; t(475)=9.375, p=2.8e-19; Supplementary Fig. S2). This indicates that the model captured individual-specific patterns of brain dynamics and explained more variance within individuals than across individuals.

We further validated the model by comparing model parameters, W and \beta, to their analogue descriptive statistics that are commonly used in the field. W was compared to an undirected functional connectivity (FC), estimated as the parcel-by-parcel Fisher’s transformed Pearson’s correlation coefficients. \beta was compared to regression coefficients estimated from a linear encoding model, that predicts neural activity from stimulus time series: x t=(encoding coefficient)×u t+(residual) (encoding coefficient 2Rn×p, residual 2Rn). We found that the estimated W was highly comparable to FC for both the SONG (cosine similarity=0.929±0.009; z=1178.36, p<0.0001 compared to shuffled chance distribution) and HCP datasets (cosine similarity=0.913±0.014; z=2830.05, p<0.0001) (Fig. 3B, D). Likewise, \beta was highly comparable to encoding coefficients for both the SONG (cosine similarity=0.324±0.100; z=448.58, p<0.0001) and HCP datasets (cosine similarity=0.243±0.036; z=1022.18, p<0.0001) (Fig. 3C).

Do these parameters capture differences between individuals as well as differences in cognitive states? Models were considered sensitive to individual differences if the parameters estimated from runs of the same person were more similar compared to parameters of different individuals (Fig. 3E). Both W and \beta were sensitive to individual differences, and W more strongly reflected individual differences than \beta (Wilcoxon rank-sum test between same vs. different individual pairs; W: z=14.957, p=1.4e-50 for SONG, z=44.667, p=0.0 for HCP; \beta: z=2.714, p=0.007 for SONG, z=13.828, p=1.7e-43 for HCP) (Fig. 3F). This aligns with prior findings that functional connectivity is stable within an individual and distinctive across individuals, thus driven more by trait than state differences33,34.

Models were considered sensitive to cognitive state differences if parameters estimated from runs of the same movie stimulus were more similar compared to runs of different stimuli (Fig. 3E). Both W and \beta were sensitive to cognitive state differences, and \beta more strongly reflected cognitive state differences than W (W: z=21.414, p=9.9e-102 for SONG, z=75.509, p=0.0 for HCP; \beta: z=45.189, p=0.0 for SONG, z=242.537, p=0.0 for HCP) (Fig. 3F). The corresponding descriptive statistics—the FC and encoding coefficients—closely followed this trend. These results indicate that the parameters estimated from our dynamical systems model are comparable to validated descriptive statistics and hold biological plausibility.

The model reveals stable attractors organized along the cortical hierarchy
In dynamical systems, the differential equation determines and predicts how the system’s trajectory would evolve over time from a given initial state, assuming that no internal or external perturbation exists beyond what is parameterized. Depending on the nature of the system, the trajectory may eventually converge to a set of stable fixed point attractors, settle onto stable periodic orbits known as limit cycles, bounce between saddle points when attractors are absent, or exhibit chaotic behavior19.

We hypothesized that our data-driven models would reveal a finite set of attractors, or stable patterns of neural activity toward which neural trajectories are likely to converge. To find attractors, we ran forward simulations of the model starting from multiple initial states, each corresponding to the observed 200-parcel neural activity pattern at every time step of the fMRI time series (Table 1). The observed neural activity at each time step was simulated 5000 time-steps forward in time, which we considered was sufficient for convergence. This method simulates the intrinsic drift, or the likely path that neural activity would follow in the absence of external perturbations beyond those parameterized by the model. If the forward simulation reached a state where no further change occurred over the last ten iteration steps, we considered it to have converged to a fixed point, representing the bottom of a valley in the inferred attractor landscape (Fig. 1). Fixed points estimated across all time steps were grouped as the same attractor if their pairwise Euclidean distance was less than 0.1 (Fig. 4A).

We focused this analysis on the HCP dataset, given that it contained a larger number of total runs (1666 runs) compared to the SONG dataset (188 runs). We included rest, task, and movie-watching runs in the analysis. When conducting forward simulations on all runs respectively, a majority of runs converged onto a set of point attractors: many converged onto 2 (41.06%) or 4 attractors (51.08%), and some converged onto 6 (7.02%) or 8 attractors (0.66%) (Fig. 4A). The stability of the attractors in all of these runs was confirmed by computing the Jacobian matrix of the model and verifying that all eigenvalues occupy the unit circle, indicating that small perturbations around an attractor would return the system to that attractor. Only 3 out of 1666 runs (0.18%) did not converge to fixed point attractors but exhibited oscillatory limit cycles. This indicates that when large-scale cortical activity evolves according to the dynamics specified by the model equations, it is most likely to fall into a set of attractors. Note that the model algorithm is designed to identify an even number of attractors, each representing opposing patterns of neural activity (Supplementary Text S1).

We predicted that these attractors would tile the core gradients of cortical organization that span between the transmodal default mode network (DMN) areas and the unimodal sensory and motor areas35,36. Furthermore, we expected that these attractors would correspond to canonical brain states—distinct and recurring patterns of neural activity—that have been replicated in multiple studies24,37. To test these hypotheses, we combined the activity patterns of all fixed point attractors estimated in every run. We applied a k-means clustering to find 4 attractor clusters, because a majority of runs exhibited either four or less attractors (Fig. 4B).

The mean activity patterns of these attractor clusters (Fig. 4C) were compared to the known cortical gradients estimated by Margulies et al.36 (Fig. 4D). This previous work applied a nonlinear dimensionality reduction algorithm on hundreds of participants’ resting-state functional connectivity data to find top gradients that explained the largest variances. The primary gradient distinguished unimodal from transmodal areas, and the secondary gradient distinguished sensory from motor areas, which were argued to be an “intrinsic coordinate system” of the human brain38 and replicated by multiple research groups39.

Our two attractor clusters were highly comparable to the top cortical gradients (\rho values=0.850 and 0.896, p values<0.0001). Specifically, the attractor clusters served as axes that separated (i) the DMN from the unimodal (UNI) sensory and motor areas and (ii) the visual network (VIS) from the sensorimotor network (SM). We did not find meaningful clusters of attractors nor resemblance to gradients when the attractors were estimated from the randomly circular-shifted fMRI data (Supplementary Fig. S4). We also found no significant participant-level differences in attractor positions, indicating that the attractor locations were largely consistent across participants (Supplementary Text S2). These results provide a dynamical systems explanation for the recurrence of a small number of brain states over time: attractors correspond to brain states, the positions of these attractors are constrained by the brain’s functional network organization, and neural dynamics unfold along this attractor landscape.

Geometry of neural dynamics differs across fluctuating attentional states during task and movie-watching
We characterized the attractor landscapes of large-scale neural dynamics, or the likely path that neural activity will take when following the model’s equations (Fig. 4A). Would the geometry of these paths vary depending on a person’s attentional state at a given moment? In other words, would one’s attentional state—from a state of being focused on a task versus zoning out to being engaged in a movie versus being bored—relate to where the brain state is positioned within the hills and valleys?

Revisiting the forward simulation in Table 1, we identified two vectors at each time step based on the observed neural activity x t: one defining the intrinsic drift, or the flow governed by regional connections and self-decays (called the “intrinsic vector” defined by W\psi \alpha (x t)-D\odot x t, shortened as W~D), and the other defining the flow nudged by the external inputs at the corresponding moment (called the “extrinsic vector” defined by \beta u t, shortened as ~B). To quantify this, for every time step, we extracted the intrinsic and extrinsic vectors and calculated their angles and magnitudes (Fig. 5A). The angle indicates the direction of the vector with respect to the position of the attractor it eventually converged, with high angle indicating the vector directing away from the attractor. The magnitude represents the degree of change from the initial state, with high magnitude indicating a fast-moving vector. Because these measures were estimated at each time step, we were able to extract their time series over the course of the fMRI run. In essence, they provide a geometric characterization of neural dynamics on the attractor landscape.

We asked whether neural trajectories toward attractors—characterized by the angle and magnitude of the intrinsic and extrinsic activity—systematically varied based on the attentional state of a participant. We focused our analysis on the SONG dataset, which contained both naturalistic movie-watching (i.e., two runs of comedy sitcom episode watching) and controlled attention tasks (i.e., two runs of gradual-onset continuous performance task; gradCPT)24. Importantly, these runs included time-varying behavioral measures of attention. During sitcom episodes, participants continuously rated how engaging they found the episode by adjusting the scale bar40. This measure was collected after the fMRI scan, as participants re-watched the same episode. (Narrative engagement has been characterized as a state of heightened emotional arousal and attentional focus40–43. Because narrative engagement accompanies changes in both arousal and attentional states, denoting them as “attentional state” is a simplification made in this article.) During gradCPT, participants pressed a button at every second whenever target images appeared inside the scanner. The inverse of response time variability served as a proxy of sustained attention, with moments of stable response times indicating high attention and variable response times indicating low attention44.

We correlated the behavioral time series of each participant with the four geometric measures’ time series, which was compared to a respective chance distribution where each geometric measure was correlated with circular-shifted behavioral time courses (Table 2 and Fig. 5B).

Both the angle and magnitude of intrinsic activity (W~D) were significantly correlated with attention dynamics (Table 2 and Fig. 5B). The relationships were comparable between repeated runs of the movie-watching and attention tasks, highlighting the reliability of our results. Interestingly, an opposite relationship to attention dynamics was found between the two contexts. The angle of the intrinsic vector was large when participants reported high engagement toward episodes, whereas the angle was small when participants performed stably in gradCPT. The magnitude of the intrinsic vector was small when participants reported high engagement toward episodes, whereas the magnitude was large when participants performed stably in gradCPT. This means that neural dynamics toward attractors—largely determined by the landscape—not only varied across attentional state dynamics but in a manner different across contexts.

On the other hand, the angle and magnitude of the extrinsic activity (~B) did not meaningfully relate to attention (Table 2). There were weak trends of correlations, which were seemingly derived from emergent correlations amongst the four geometric measures (Supplementary Fig. S5). This highlights that attention relates to changes in intrinsic neural dynamics, but not stimulus-driven neural dynamics.

In a study by Song et al.24 in which this original data was collected, different brain states were associated with high attention in movie-watching and attention task contexts. Participants reported high engagement to sitcom episodes during the “base” state, a state where no functional network exhibited dominant activity and was positioned at the center of the latent manifold. On the contrary, participants exhibited stable task performance to gradCPT during the DMN state, a state defined by high activity in the DMN (consistent with findings by refs. 45–47). Motivated by these results, we categorized attractors into either of the four ends of gradients 1 (DMN and UNI) and 2 (VIS and SM), based on the neural pattern similarity between attractors and the gradients. We assessed the relationships between the geometric measures and attention dynamics, separately for these four attractors. The goal was to see if the relationship, shown in Fig. 5B and Table 2, differs depending on where the neural trajectories converge.

We primarily found that the likelihood of attractors differed across the two contexts (Fig. 5C). The neural activity was nearly equally likely to fall into one of the four attractors during movie-watching, whereas the ends of gradient 1, the DMN and UNI, were much more likely to serve as attractors during sustained attention tasks, compared to the ends of gradient 2. This finding showing different likelihood of attractors across contexts indicates that the attractor landscapes differed across contexts.

During movie-watching, we found that the main relationship between the geometric measures and attention dynamics remained consistent, irrespective of which attractor the neural activity fell into. When participants reported high engagement to movies, the intrinsic drift directed away from the attractors with decreased magnitude (Fig. 5D). This implies that the brain was in a flattened or shallow region of the landscape when people were attentive, such that the brain activity was more likely to lie at the center of the manifold and gravitated less toward the attractors (Fig. 6A). On the contrary, during attention tasks, the main effect we found in Table 2 was specific to when the neural activity converged onto the DMN attractor (Fig. 5D). When attentive, neural activity approached the DMN attractor faster and more directly. This implies that the brain was on a steeper region of the attractor landscape near the DMN attractor when attentive (Fig. 6B).

Together, the results suggest that the neural dynamics toward attractors differed across high and low attentional states, in a context-dependent manner. This indicates that the attractor landscape of large-scale cortical dynamics systematically varies depending on attentional states and task demands.

Methods
Model description
The dynamical systems model used in this study is adopted from the neural mass model called the mesoscale individualized neurodynamic (MINDy) model9,19. The model is designed to fit the neural activity time series—whichever units or scales the neural activities are sampled from—and the time-aligned experimental variables, such as task, stimulus, or behavior. For our use, the neural activity corresponded to the BOLD activity of 200 cortical parcels collected from human fMRI. The experimental variable corresponded to audiovisual and semantic feature embeddings of the movies that participants watched inside the scanner and visual feature embeddings of the images that were presented as task stimuli. However, the choice of neural units and experimental variables can vary depending on the study and research question. No tailoring specific to the fMRI data (e.g., deconvolution of the hemodynamic response function) was made.

To reiterate, the model is defined as the following equation, with the neural activity at time t represented as x t (x t 2Rn, n= 200 parcels) and stimuli at time t represented as u t (u t 2Rp, p= 100 PCs of the feature embeddings), with \Delta t set to 1 TR.

x^ t+1 = x t + W\psi \alpha (x t)-D\odot x t + \beta u t (4)

\psi \alpha (x t) = \sqrt{\alpha^2+(bx t+0.5)^2} - \sqrt{\alpha^2+(bx t-0.5)^2} (5)

The weight matrix (W 2Rn×n) represents directional interaction between neural units, which conceptually corresponds to effective connectivity. W\psi \alpha (x t) represents the weighted sum of the neural unit’s nonlinearly transformed activity at time t multiplied by its directed connections from every other neural unit. To constrain the estimation of W as the sum of sparse random component and a low-rank structured component60, we decomposed the parameter into W=W S+W L where W S (W S 2Rn×n) represents a sparse matrix after L 1 regularization and W L=W 1 W T 2 (W 1,2 2Rn×r) is given by low-rank approximation (r=n=3) with sparsity also given to W 1 and W 2 with L 1 regularization. This formalization allows the model to capture global low-dimensional motifs of connectivity superimposed on a sparse network, which gives rise to the known low-dimensional and modular structure of large-scale functional organization.

Nonlinearity follows a parametrized sigmoidal transfer function (\psi a), which maps neural activity to a bounded output at a range of −1 to 1. The slope of the transfer function is determined by the estimated \alpha. b is fixed at 20/3 in our model.

\beta 2Rn×p represents a linear coupling between neural activity and incoming stimuli that are time-aligned with one another. The time was aligned by convolving the stimulus time series with the canonical hemodynamic response function. \beta conceptually corresponds to a linear mapping from the stimulus space to the neural space.

A decay term (D 2Rn) represents convergence to baseline activity at the absence of external inputs or interactions, with \odot denoting element-wise product. D is an algorithmically important parameter, because D is initially set to a value significantly higher than 1, which flips the right-hand side of the equation toward a large negative factor of x t. This accentuates the difference in neural activity between units, thus allowing effective estimation of the directed connectivity W. From our empirical tests, if D is set to a biologically plausible value of <1, the estimated W does not recapitulate the descriptive functional connectivity measure.

Model fitting
We fit the neural activity time series acquired at each run, in batches of size 300 TRs. Specifically, we predicted the neural activity of consecutive timesteps x^ 2:T based on the observed neural activity x 1:T-1 (where T=300). Model parameters were optimized across 2500 iterations using stochastic gradient descent, specifically the Nesterov-accelerated Adaptive Moment Estimation optimizer as chosen in previous studies9,19. The loss L was calculated as follows. \Lambda represents the regularization term, which enforces sparsity in connections. Regularization terms were fixed to \lambda 1=0.075, \lambda 2=0.2, \lambda 3=0.05.

L = 1/2 ||x^ 2:T - x 2:T||2 2 + \Lambda (6)

\Lambda = \lambda 1||W S|| 1 + \lambda 2 Tr(|W S|) + \lambda 3 ||W 1|| 1 + ||W 2|| 1 (7)

Because D was initialized at a value higher than 1, the prediction accuracy (i.e., rank correlation between the predicted and observed neural activity patterns) started at a value near −1 at the start of the training. The prediction accuracy increased gradually across iterations. In contrast, W quickly became comparable to a functional connectivity matrix in the initial phase of training, but across more iterations, it gradually approached a diagonally dominant matrix, which is conceptually similar to the first-order autoregressive model. To prioritize biologically meaningful parameter optimization rather than brain activity prediction, we stopped the training at 2500 iterations (in batches of 300 consecutive TRs) at which point W was similar to functional connectivity, but the model’s accuracy still remained negative. To boost prediction accuracy, we fit an ordinary least squares linear regression model to estimate pW, pD, and pB, which are scalar values (pW, pD, pB 2R). Parameters W, D, and \beta were scaled by these value estimates, respectively.

x_ t = pW*W\psi \alpha (x t) - pD*D\odot x t + pB*\beta u t (8)

Model fitting and ensuing analyses were conducted in Python (v3.12.2) and Pytorch (v2.4.1).

Hyperparameters
Hyperparameter selection largely followed the original implementation of the MINDy model9, but with some simplifications made (Table 3).

fMRI datasets
Two openly available fMRI datasets were analyzed: the SONG dataset (participants recruited in South Korea; N=27) and the HCP dataset (participants recruited in the USA; N=119). Both datasets were approved by the institutional review boards and de-identified for sharing. No information on participants’ age, sex, or gender was collected in the present study. Preprocessing steps followed Song et al.24. We applied a 200-parcel cortical atlas by Schaefer et al.29 where the BOLD activities of voxels corresponding to each parcel were averaged to represent the parcel activity9,19. The SONG dataset includes two runs of resting-state, two runs of controlled sustained attention task called the gradual-onset continuous performance task (gradCPT), two runs of comedy sitcom-watching, and one run of educational documentary-watching (3T scans with TR=1 s). The HCP dataset includes four runs of resting-state in 3T (TR=0.72 s) and four in 7T (TR=1 s), two runs of working memory tasks in 3T (6 different types of cognitive task runs), and four runs of movie-watching in 7T. The total duration, ranging from 405 to 1486 TRs.

Input features
Low-level visual features were characterized by hue, saturation, and pixel intensity, which were estimated per frame and averaged across frames within an event (rgb2hsv function in MATLAB R2024a). In the gradCPT run with face images, hue and saturation were excluded from the analyses because the images were presented in grayscale. Low-level audio features were represented with amplitude and pitch of left and right stereos, which were estimated per frame and averaged across frames within an event (audioread and pitch functions in MATLAB).

Visuo-semantic features were represented by 512-dimensional embeddings of OpenAI’s pretrained Contrastive Language-Image Pre-training model61 (huggingface; clip-vit-base-patch32). For audio-semantic features, we first applied OpenAI’s WhisperX model to the audio file to transcribe speech and dialogues, along with their time stamps62. Transcripts were sampled at a TR resolution and were transformed using 512-dimensional embeddings of Google’s Universal Sentence Encoder63 (tensorflow v2.19.0; https://tfhub.dev/google/universal-sentence-encoder-multilingual/3). A mean value of the respective dimension was assigned to moments when speech or dialogue did not exist.

Embedding time series were convolved with a canonical hemodynamic response function and normalized across time per dimension. Principal component analyses were conducted on the concatenated embedding time series of the 3 movie-watching runs in SONG and 4 runs in HCP, respectively. We selected the top 100 principal component time series, which explained 74.97% and 77.69% of variance, respectively. Principal component analysis was conducted on the gradCPT face run in isolation, given that the same images were presented in the same sequence for all participants (100.00% of explained variance for 100 principal components). For gradCPT scene runs, because presented images and sequences differed for all participants, participant-specific embeddings were concatenated for a principal component analysis (85.34% of explained variance). Resulting principal component time series were again normalized across time per dimension to be used as stimulus input u in the model.

Validation of model parameters
Parameters were estimated through training the model on the data collected from an individual fMRI run (2500 training iterations). After parameters were fixed, we predicted the next-time-step neural activity from the observed neural activity time series to calculate explained variance (r2). Note that the data used for training and testing were the same. When comparing parameter W with an undirected functional connectivity matrix (parcel-by-parcel Fisher’s transformed Pearson’s correlation coefficients), we took the average of the upper and lower triangles of W (region i→j and region j→i in a directed graph) and took the cosine similarity with the edge strengths of the functional connectivity matrix. Due to the low temporal resolution of the fMRI data, our W estimates are largely symmetrical. Encoding coefficients were estimated using an ordinary least squares regression with a residual term. Cosine similarity between all values in \beta and values in encoding coefficients were computed. For both metrics, values in descriptive statistics were randomly shuffled 10,000 times, which served as the respective chance distribution. Z statistics and two-tailed p values were calculated with respect to chance distributions. Cosine similarities were calculated between parameter estimate of run pairs. The cosine similarity values were grouped into whether they correspond to the same vs. different individuals or the same vs. different movie stimuli, which were compared using Wilcoxon rank sum tests.

Attractor clusters and gradients
Attractors were estimated by forward simulating the observed neural activity pattern at each time step based on the estimated model parameters (Table 1). For each time step, we performed 5000 simulations. If there was no longer a change in the last 10 iterations of the 5000th forward simulations (||x_|| 1<1e-6), we considered the system to have converged to a fixed point. If the Euclidean distance between the pair of estimated fixed points obtained from different time steps was less than 0.1, those were grouped as the same attractor. Using this approach, we found that most fMRI runs converged to 2 or 4 attractors. Principal component analysis was applied to each individual fMRI run only for visualization purpose in a 3D space (Fig. 4A).

We then aggregated attractors (each defined as a neural activity pattern across 200 parcels) across all fMRI runs and applied k-means clustering to identify 4 attractor clusters. The choice of k=4 was motivated by the observation that most runs exhibited either 2 or 4 attractors, with higher numbers occurring less frequently. The mean activity pattern of each cluster represented the attractor cluster (Fig. 4B, C). Two pairs of clusters exhibited a pattern correlation of −1. Principal component analysis was also applied to the aggregated attractor patterns, only for visualization purpose (Fig. 4B).

Cortical voxel gradients estimated by Margulies et al.36 were obtained from NeuroVault (https://identifiers.org/neurovault.collection:1598). Gradient values were averaged within each parcel to generate gradient estimates of the 200 parcels (Fig. 4D). The neural activity patterns of the attractor clusters were then compared to the parcel-level gradient values using Spearman’s rank correlation.

Relating speed and direction of neural dynamics toward attractors with measures of attention
With the neural activity pattern at each time step as an initial state (x t), we considered one-step forward simulation of x^ t+1 = x t + W\psi \alpha (x t)-D\odot x t + \beta u t and separated vectors representing internally-driven change (W\psi \alpha (x t)-D\odot x t, shortened to W~D) and externally-driven change (\beta u t, shortened to ~B). The angle between the respective vector and the eventual end state was calculated with the following equation,

\theta = arccos((v1*v2)/(||v1||*||v2||)) (9)

where v1*v2 represents the dot product of the two vectors, and ||v1|| and ||v2|| are the Euclidean norms of the vectors. Before applying the inverse cosine, the cosine similarity was clipped to the range of [−1, 1]. The magnitude of the vector was calculated as the Euclidean norm. Because these measures were estimated at each time step, repeating this across all time steps within a run resulted in their time series.

To probe attentional state changes during movie-watching, Song et al.24 asked participants to continuously rate their engagement, on a scale of 1 to 9, as they re-watched the movies after the fMRI scan. Each participant’s engagement rating was normalized across time and convolved with a hemodynamic response function. To probe attentional state changes during gradCPT, we analyzed participants’ button response times as in the previous study24. After linear interpolation of no button response trials and regressing out the linear trend, we calculated response time variability by taking the deviance from the mean response time at every TR. Because studies showed that moments of low response time variability correspond to high sustained attention or optimal task performance, the inverse response time variability time course was used as a proxy for attention dynamics44. Again, the inverse response time variability was normalized across time and convolved with a hemodynamic response function. See Song et al.24 for details of behavioral experiments and data analyses.

Angle and magnitude time series were correlated with attention measures sample data TR resolution using Spearman’s rank correlation, for each individual run. The mean of all participants’ Fisher’s r-to-z transformed \rho values were compared to permuted null distributions to calculate z statistics, where behavioral time series were circular-shifted across time with a random multiplication of either 1 or −1 (10,000 iterations). FDR correction was applied across 16 significance tests (Table 2).

The relationship between geometric measures and attention was analyzed separately depending on the position of the attractor at each time step, or where the eventual end state lies at the final round of forward simulation. Because we found that attractor clusters lie at the ends of the known primary and secondary gradients (Fig. 4), we categorized the attractor of each time step to one of the four ends of the two gradients. The attractor was labeled as either the DMN, UNI, VIS, or SM, based on the highest correlation coefficient. Correlation between the geometric measure and attention was calculated from a subset of time series corresponding to the respective attractor cluster category. In a majority of participants’ gradCPT runs, VIS and SM attractors did not appear. Because fewer than 10 fMRI participants’ gradCPT face or gradCPT scene runs included VIS and SM attractors, statistical analysis was not conducted for these cases. Significance was again tested by comparing the mean of Fisher’s r-to-z transformed \rho values to null distributions created from circular-shifted behavioral time series. Twelve comparisons were corrected for in Fig. 5D, and 24 comparisons were corrected for in Supplementary Fig. S6B using FDR correction.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors have addressed all my concerns.", "reviewer_2": "The authors have sufficiently addressed all my concerns. I support the manuscript for\npublication and congratulate the authors on their work.", "reviewer_3": "I sincerely appreciate the authors’ very thorough responses to the reviewer comments. In\nparticular, the inclusion of the Jacobian analysis and null-model testing, together with the\nmany clarifications and discussions added to the main text, subsequently strengthens the\ncredibility of the findings. Also, I believe these additions meaningfully expand the scope and\ndepth of the manuscript’s discussion.\nI have one minor comment regarding the choice of k in the k-means analysis. I appreciate the\nauthors for performing the three quantitative validation tests. As is often the case, these\nmetrics favor the smallest nontrivial solution (k = 2) in terms of absolute scores, and are\ntherefore typically interpreted by examining elbows/knees or local extrema. From the results, k\n= 6 appears suboptimal, whereas k = 4 (and maybe k = 8) can be plausible choices. Yet the\ncurves do not exhibit a well-defined elbow/knee, and the authors’ a priori justification for k = 4\nis reasonable."}, "author_response": "We are glad that our response has addressed both reviewers’ prior concerns. The two reviewers’\nsuggestions were extremely helpful in improving our manuscript.\nReviewer #3 (Remarks to the Author):\nI sincerely appreciate the authors’ very thorough responses to the reviewer comments. In\nparticular, the inclusion of the Jacobian analysis and null-model testing, together with the\nmany clarifications and discussions added to the main text, subsequently strengthens the\ncredibility of the findings. Also, I believe these additions meaningfully expand the scope and\ndepth of the manuscript’s discussion.\nWe agree that the Jacobian analysis and null-model testing, along with others, significantly improved\nthe paper. We thank the reviewer for helpful suggestions.\nI have one minor comment regarding the choice of k in the k-means analysis. I appreciate the\nauthors for performing the three quantitative validation tests. As is often the case, these\nmetrics favor the smallest nontrivial solution (k = 2) in terms of absolute scores, and are\ntherefore typically interpreted by examining elbows/knees or local extrema. From the results, k\n= 6 appears suboptimal, whereas k = 4 (and maybe k = 8) can be plausible choices. Yet the\ncurves do not exhibit a well-defined elbow/knee, and the authors’ a priori justification for k = 4\nis reasonable.\nWe appreciate the reviewer’s understanding our a priori justification for using k=4. We had a clear\ndata-driven motivation to examine k=4, so did not conduct these tests in our original analyses.\nRunning the validation tests showed that the metrics favored the smallest solution with some hint of\nelbows/knees for k=4 and k=8. We have chosen the original k=4 in the manuscript and did not further\nchange the manuscript."}]

Decision letter