Methods

Data¶

We obtained functional MRI data from 7 different sources (Table 1). We included three resting state studies with healthy volunteers (study 1, study 2, study 3, $n_{total}=118$), one task-based study (study 4, $n_{total}=33$ participants, 9 runs each), an individual participant meta-analytic activation map of pain (study 5, $n_{total}=603$ from 20 different studies), 8 task-based activation patterns obtained from coordinate-based meta-analyses via Neurosynth (study 6, 14371 studies in total, see Supplementary Table 1) and a resting state dataset focusing on Autism Spectrum Disorder (ASD) from the ABIDE (Autism Brain Imaging Data Exchange, study 6, $n_{total}=1112$, Di Martino et al. (2014)).

Study 1 was used to evaluate the potential of the resting state functional brain connectome to be considered as an attractor network, to optimize the temperature (β) and noise (σ) parameters of the fcHNN model, and to evaluate the of the proposed approach to reconstruct resting state brain dynamics. Study 2 and 3 served as replications studies for these analyses. Study 1, 2 and 3 is well suited to examine replicability and generalizability; data in these three studies were acquired in 3 different centers from 2 different countries, by different research staff, with different scanners (Philips, Siemens, GE) and different imaging sequences. Further details on study 1-3 are described in Spisak et al. (2020). The ability of the proposed approach to model task-based perturbation of brain dynamics was evaluated in Study 4, which consisted of nine task-based fMRI runs for each of the 33 healthy volunteers. In all runs, participants received heat pain stimulation. Each stimulus lasted 12.5 seconds, with 3-second ramp-up and 2-second ramp-down periods and 7.5 seconds at target temperature. Six levels of temperature were administered to the participants (level 1: 44.3°C; level 2: 45.3°C; level 3: 46.3°C; level 4: 47.3°C; level 5: 48.3°C; level 6: 49.3°C). We used run 1 (passive experience), run 3 (down-regulation) and run 7 (up-regulation) from this study. In runs 3 and 7, participants were asked to cognitively “increase” (regulate-up) or “decrease” (regulate-down) pain intensity. No self-regulation instructions were provided in run 1. See Woo et al. (2015) for details. Pain control signal for our task-based trajectory analyses on data from study 4 was derived from our individual participant meta-analysis of 20 pain fMRI studies (study 5, n=603). For details, see Zunhammer et al. (2021). To obtain fMRI activation maps for other tasks, we used Neurosynth(Tor D. (2011)), a web-based platform for large-scale, automated synthesis of functional magnetic resonance imaging (fMRI) data. We performed 8 different coordinate-based meta-analyses with the terms “motor”, “auditory”, “visual”, “face”, “autobiographical”, “theory mind”, “language” and “pain” (Supplementary Table 1) and obtained the Z-score maps from a two-way ANOVA, comparing the coordinates reported for studies with and without the term of interest, and testing for the presence of a non-zero association between term use and voxel activation. In study 7 (ABIDE), we obtained preprocessed regional timeseries data from the Preprocessed Connectome Project Craddock et al., 2013, as shared (https://osf.io/hc4md) by Dadi et al. (2019). Preprocessed timeseries data were obtained with the 122-region version of the BASC (Bootstrap Analysis of Stable Clusters) brain atlas Bellec et al., 2010.

Preprocessing and timeseries extraction¶

Functional MRI data from studies 1-4 was preprocessed with our in-house analysis pipeline, called the RPN-pipeline (https://github.com/spisakt/RPN-signature). The RPN-pipeline is based on PUMI (Neuroimaging Pipelines Using Modular workflow Integration, https://github.com/pni-lab/PUMI), a nipype-based Gorgolewski et al., 2011 workflow management system. It capitalizes on tools from FSL Jenkinson et al., 2012, ANTS Avants et al., 2011 and AFNI Cox, 1996, with code partially adapted from the software tools C-PAC Craddock et al., 2013 and niworkflows Esteban et al., 2019, as well as in-house python routines.

Brain extraction from both the anatomical and the structural images, as well as tissue-segmentation from the anatomical images was performed with FSL bet and fast. Anatomical images were linearly and non-linearly co-registered to the 1mm-resolution MNI152 standard brain template brain with ANTs (see https://gist.github.com/spisakt/0caa7ec4bc18d3ed736d3a4e49da7415 for parameters). Functional images were co-registered to the anatomical images with the boundary-based registration technique of FSL flirt. All resulting transformations were saved for further use. The preprocessing of functional images happened in the native image space, without resampling. Realignment-based motion-correction was performed with FSL mcflirt. The resulting six head motion estimates (3 rotations, 3 translations), their squared versions, their derivatives and the squared derivatives (known as the Friston-24-expansion, Friston et al. (1996)) were calculated to be used as nuisance signals. Additionally, head motion was summarized as framewise displacement (FD) timeseries, according to Power’s method Power et al., 2012, to be used in data censoring and exclusion. After motion-correction, outliers (e.g. motion spikes) in timeseries data were attenuated using AFNI despike. The union of the eroded white-matter maps and ventricle masks were transformed to the native functional space and used for extracting noise-signal for anatomical CompCor correction Behzadi et al., 2007.

In a nuisance regression step, 6 CompCor parameters (the 6 first principal components of the noise-region timeseries), the Friston-24 motion parameters and the linear trend were removed from the timeseries data with a general linear model. On the residual data, temporal bandpass filtering was performed with AFNI’s 3DBandpass to retain the 0.008–0.08Hz frequency band. To further attenuate the impact of motion artifacts, potentially motion-contaminated time-frames, defined by a conservative FD>0.15mm threshold, were dropped from the data (known as scrubbing, Satterthwaite et al. (2013)). Participants were excluded from further analysis if more than 50% of frames were scrubbed.

The 122-parcel version of the BASC (Multi-level bootstrap analysis of stable clusters) multi-resolution functional brain atlas Bellec et al., 2010 was individualized; it was transformed to the native functional space of each participant (interpolation: nearest neighbour) and masked by the grey-matter mask obtained from the anatomical image, to retain individual grey-matter voxels only. Voxel-timeseries were averaged over these individualized BASC regions.

All these preprocessing steps are part of the containerized version of the RPN-pipeline (https://spisakt.github.io/RPN-signature), which we run with default parameters for all studies, as in Spisak et al. (2020).

Functional connectome¶

Regional timeseries were ordered into large-scale functional modules (defined by the 7-parcel level of the BASC atlas) for visualization purposes. Next, in all datasets, we estimated study-level mean connectivity matrices by regularized partial correlation, via the Graphical Lasso algorithm that estimates a sparse precision matrix by solving a Lasso problem and an L1 penalty for sparsity Varoquaux et al., 2010, as implemented in nilearn Abraham et al., 2014. Diagonal elements of the matrices were set to zero.

Connectome-based Hopfield networks¶

Hopfield networks, a type of artificial neural network, consist of a single layer of $m$ fully connected nodes Hopfield, 1982, with activations $\bold{a} = (a_1, \dots, a_m )$. Hopfield networks assign an energy to any arbitrary activity configurations:

where $W$ is the weight matrix with element $w_{i,j}$ denoting the weight between nodes i and j and $\bold{b}$ is the bias, which is set to $\bold{b} = 0$ for all experiments.

During the so-called relaxation process, the activities of the nodes are iteratively updated until the network converges to a stable state, known as an attractor state. The dynamics of the network are governed by the equation referenced as eq. (1) of the main text, or in matrix form:

where $\bold{a'} = ({a'}_1, \dots, {a'}_m)$ is the activity in the next iteration and $S(.)$ is the sigmoidal activation function ($S(a) = tanh(a)$ in our implementation) and β is the temperature parameter. During the stochastic relaxation procedure, we add weak Gaussian noise to each node’s activity at every iteration:

where $\epsilon \sim \mathcal{N}(\mathbf{\mu}, \sigma)$, with σ regulating the amount of noise added and $\mathbf{\mu}$ set to 0, by default.

In this work we propose functional connectome-based Hopfield neural networks (fcHNNs) as a model for large-scale brain dynamics. FcHNNs are continuous-state Hopfield neural networks with each node representing a brain region and weights initialized with a group-level functional connectivity matrix. The weights are scaled to zero mean and unit standard deviation.

fcHNN convergence and attractors¶

We investigated the convergence properties of functional connectome-based HNNs in study 1 by contrasting the median number of iterations until reaching convergence to a permutation based null model. In more detail, the null model was constructed by randomly permuting the upper triangle of the original connectome and filling up the lower triangle to get a symmetric network (symmetry of the weight matrix is a general requirement for convergence). This procedure was repeated 1000 times. In each repetition, we initialized both the original and the permuted fcHNN with random input and counted the number of iterations until convergence. The whole procedure was repeated with $\beta=0.3, 0.35, 0.4, 0.45, 0.5, 0.55$ and 0.6 (providing 2-8 attractor states). In studies 1-3, we obtained the finite number of attractor states with the same values of β by repeatedly (10⁵ times) initializing the fcHNNs with random activations and relaxing them until convergence.

fcHNN projection¶

We mapped out the fcHNN state-space by initializing our fcHNN model with a random input, and applying the stochastic update step for 10⁵ iterations and storing all visited activity configurations. We performed a principal component analysis (PCA) on the state samples, and proposed the first two principal component (PCs) as the coordinate system for the fcHNN projection. Using a Multinomial Logistic Regression, we predicted to which attractor state each of the state samples converges to, using the first two PCs as features. The model’s performance was evaluated with 10-fold cross-validation. We visualized the attractor states position in the projection as well as the decision boundaries between the attractor states, based on this regression model.

Parameter optimization¶

Based on the convergence analysis, we selected the β value that provided the fastest median convergence ($\beta = 0.04$) for all subsequent analyses, to minimize computational costs. We optimized the noise parameter σ by comparing the distribution of the state-space of the fcHNN with the distribution of the real fMRI data in the fcHNN projection. Specifically, we evaluated eight different σ values, spaced logarithmically between 0.1 and 1, obtained the fcHNN projection for each model, transformed the real data into this projection and compared the distribution of the data, as well as the state occupancies. The similarity of the distribution of of the real and fcHNN-generated data was quantified as the Pearson’s correlation coefficient between the 2-dimensional histograms over a 2-dimensional grid of 100×100 uniformly distributed bins in the [-6,6] range (arbitrary units on the fcHNN projection axes) after applying a Gaussian smoothing with a σ of 5 bins.

State occupancies for each attractor state were calculated as the ratio of time spent on the basis of the attractor (both for the fcHNN-generated and the empirical data). Similarity of fcHNN-generated state occupancies to those observed in the real data was evaluated with a $\chi^2$-test statistic. Both test statistics were contrasted to two types of null models. The first null model was constructed to investigate whether fcHNNs can extract information from the functional connectome over and beyond what is already present in the co-variance structure of the data, next to normality assumptions. This was implemented by drawing random samples from a multivariate normal distribution with the functional connectome as the covariance matrix. The second null model was constructed to investigate whether the fcHNN model’s performance can be explained solely with the spatial autocorrelation properties of the data. This null model was implemented by spatially autocorrelation-preserving randomization of the real data. Timeframes were first Fourier-transformed into the frequency domain, the phases were randomized simultaneously and the data was transformed back with the inverse-Fouurier transformation. More details on the null models can be found in Supplementary Figure 5. Both null models were used to generate 1000 surrogates of the empirical data, which were projected into the fcHNN projection space and compared to the fcHNN-generated data with the above described test statistics. P-values were calculated by contrasting the real test statistics to the null-distributions As a result of this procedure, we selected $\sigma = 0.37$ for all subsequent analyses.

Replicability¶

We obtained the four attractor states in study 1, as described above. We then constructed two other fcHNNs, based on the study-mean functional connectome obtained in studies 2 and 3 and obtained all attractor states of these models, with the same parameter settings ($\beta = 0.04$ and $\sigma = 0.37$) as in study 1. In both replication studies we found four attractor states. The spatial similarity of attractor states across studies was evaluated by Pearson’s correlation coefficient.

Evaluation: resting state dynamics¶

To evaluate the explanatory power of the fcHNN projection, we performed PCA on the preprocessed fMRI time-frames from study 1 (analogously to the methodology of the fcHNN projection, but on the empirical timeseries data). Next, we fitted linear regression models which used the first two fcHNN or real data-based PCs as regressors to reconstruct the real fMRI time-frames. In-sample explained variances and the corresponding confidence intervals were calculated for both models with bootstrapping (100 samples). To evaluate the out-of-sample generalization of the PCs (fcHNN- and real data-based) from study 1, we calculated how much variance they can explain in study 2.

Similarity between state occupancy and distribution was calculated suring the parameter optimization. More detail on the associated null-models can be found in Supplementary figure 5.

To confirm that the real and fcHNN temporal sequences (from the stochastic relaxation) display similar temporal autocorrelation properties, we compared both to their randomly shuffled variant with a “flow analysis”. First we calculated the direction on the fcHNN projection plane between each successive TR (a vector on the fcHNN projection plane for each TR transition), both for the empirical and the shuffled data. Next, we obtained a two-dimensional binned means for both the x and y coordinates of these transition vectors (pooled across all participants), calculated over a 2-dimensional grid of 100×100 uniformly distributed bins in the [-6,6] range (arbitrary units) and applied a Gaussian smoothing with a σ of 5 bins (same approach as in described in the Parameter optimization section). Finally, we visualized the difference between the binned-mean trajectories of the empirical and the shuffled data as a “streamplot”, with the Python package matplotlib. The same approach was repeated with the fcHNN-generated data. The similarity of the real and fcHNN-generated flow analysis was quantified with Pearson’s correlation coefficient, p-values were obtained with permutation testing.

Evaluation: task-based dynamics¶

We used study 4 to evaluate the ability of the fcHNN approach to capture and predict task-induced alterations in large-scale brain dynamics. First, runs 1, 3 and 7, investigating the passive experience and the down- and up-regulation of pain, respectively, were preprocessed with the same workflow used to preprocess studies 1-3 (Preprocessing and timeseries extraction). Regional timeseries data was grouped to “pain” and “rest” blocks, with a 6-second delay to adjust for the hemodynamic response time. All activation timeframes were transformed to the fcHNN projection plane obtained from study 1. Within-participant differences of the average location on the fcHNN projection was calculated and visualized with radial plots, showing the participant-level mean trajectory on the projection plane from rest to pain, denoted with circles, as well as the group level trajectory (arrow). The significance of the position difference and energy difference of the participant-level mean activations in the projection plane was tested with a permutation test. We used the L2 norm of the two-dimensional position difference and the absolute energy difference, respectively, as test statistics. The permutation tests were run with 1000 permutations, randomly swapping the conditions within each participant.

To further highlight the difference between the task and rest conditions, a “flow analysis” was performed to investigate the dynamic trajectory differences between the conditions rest and pain. The analysis method was identical to the flow analysis of the resting sate data (Evaluation: resting state dynamics). First we calculated the direction in the projection plane between each successive TR during the rest conditions (a vector on the fcHNN projection plane for each TR transition). Next, we obtained a two-dimensional binned means for both the x and y coordinates of these transition vectors (pooled across all participants), calculated over a 2-dimensional grid of 100×100 uniformly distributed bins in the [-6,6] range (arbitrary units) and applied Gaussian smoothing with a σ 5 bins. The same procedure was repeated for the pain condition and the difference in the mean directions between the two conditions was visualized as “streamplots” (using Python’s matplotlib). We used the same approach to quantify the difference in characteristic state transition trajectories between the up- and downregulation conditions. The empirically estimated trajectory differences (from real fMRI data) were contrasted to the trajectory differences predicted by the fcHNN model from study 1.

To obtain fcHNN-simulated state transitions in resting conditions, we used the stochastic relaxation procedure ((3)), with $\mathbf{\mu}$ set zero. To simulate the effect of pain-related activation on large-scale brain dynamics, we set $\mu_i$ during the stochastic relaxation procedure to a value representing pain-elicited activity in region i. The region-wise activations were obtained calculating the parcel-level mean activations from the meta-analytic pain activation map from Zunhammer et al., 2021, which contained Hedges’ g effect sizes from an individual participant-level meta-analysis of 20 pain studies, encompassing a total of n=603 participants. The whole activation map was scaled with five different values ranging from 10^-3 to 10^-1, spaced logarithmically, to investigate various signal-to-noise scenarios. We obtained the activity patterns of 10⁵ iterations from this stochastic relaxation procedure and calculated the state transition trajectories with the same approach used with the empirical data. Next we calculated the fcHNN-generated difference between the rest and pain conditions and compared it to the actual difference through a permutation test with 1000 permutations, randomly swapping the conditions within each participant in the real data and using Pearson’s correlation coefficient between the real (permuted) and fcHNN-generated flow-maps as test statistic. From the five investigated signal-to-noise values, we chose the one that provided the highest similarity to the real pain vs. rest trajectory difference.

When comparing the simulated and real trajectory differences between pain up- and downregulation, we used the same procedure, with two differences. First, when calculating the simulated state transition vectors for the self-regulation conditions, we used the same procedure as for the pain condition, but introduced and additional signal in the nucleus accumbens, with a negative and positive sign, for up- and downregulation, respectively. We did not optimize the signal-to-noise ratio for the nucleus accumbens signal but, instead, simply used the value optimized for the pain vs. rest contrast (For a robustness analysis, see Supplementary figure 7).

Clinical data¶

To demonstrate the sensitivity of the fcHNN approach to clinically relevant alterations of large-scale brain dynamics in Autism Spectrum Disorder (ASD), we obtained data from n=172 individuals, acquired at the New York University Langone Medical Center, New York, NY, USA (NYU) as shared in the Autism Brain Imaging Data Exchange dataset (study 7: ABIDE, Di Martino et al., 2014. We focused on the largest ABIDE imaging center to ensure that our results are not biased by center effects. We excluded high motion cases similarly to our approach in studies 1-4, i.e. by ignoring (“scrubbing”) volumes with FD>0.15 and excluding participants with more than 50% of data scrubbed. Timeseries data was pooled and visualized on the fcHNN projection of study 1, separately for ASD and control participants. Next, for each participant, we grouped the timeframes from the regional timeseries data according to the corresponding attractor states (obtained with the fcHNN model from study 1) and averaged timeframes corresponding to the same attractor state to calculated participant-level mean attractor activations. We assessed mean attractor activity differences between the patient groups with a permutation test, randomly re-assigning the group labels 50000 times. We adjusted the significance threshold with a Bonferroni-correction, accounting for tests across 4 states and 122 regions, resulting in $\alpha = 0.0001$. Finally, we have calculated the trajectory differences between the two groups, as predicted by the group-specific fcHNNs (initialized with the ASD and TCD connectomes), and - similarly to the approach used in study 4 - we contrasted the fcHNN predictions to the trajectory differences observed in the real rsfMRI data. As in the previous flow analyses, we tested the significance of the similarity between the predicted and observed trajectory differences with a permutation test (1000 permutations), by shuffling group labels.