Spectral graph model for fMRI: A biophysical, connectivity-based generative model for the analysis of frequency-resolved resting-state fMRI

Understanding the governing principles underlying the brain’s resting-state activity is an area of immense importance. Since its first documentation in 1990 (Ogawa et al., 1990), the blood oxygenation level-dependent (BOLD) signal in functional magnetic resonance imaging (fMRI) has been a critical and quickly evolving method to study changes in brain activity during task and rest. More recently, the synchronized fluctuations between brain regions have been characterized as a type of functional connectivity (FC) (Fox & Raichle, 2007). This has itself launched yet another paradigm shift toward analyzing BOLD signals as a network, where connections across regions are defined as the Pearson’s correlation between the average BOLD signal time series of those respective regions (Bullmore & Sporns, 2009). This network science framework has opened lines of inquiry that seek to model and, therefore, predict empirical FC using only the brain’s wiring diagram, known as its structural connectivity (SC), as measured by diffusion-weighted MRI. Despite these successes, resting-state (rs-)fMRI analysis methods have reached saturation, and the next generation of techniques will need to address two broad limitations: the focus on the fMRI signal’s covariance (FC) without considering its power spectrum and the use of statistical rather than biophysical models to obtain features of interest. Since every limitation provides an opportunity for a paradigm shift, below we describe the opportunities available for future methodological innovation. We then introduce our proposal, which we believe could be an early but important step in this direction.

Current methods broadly seek to model, analyze, or extract features of interest from functional connectivity (FC)—that is, the second-order covariance structure inherent in the data. Whether this FC is evaluated via pairwise Pearson’s correlation between time series, or through ICA to achieve functional connectivity networks (FCNs), it essentially neglects the actual temporal structure of fMRI time series, specifically their frequency content. The BOLD signal is generally bandpass filtered from 0.01 to 0.08 Hz to avoid contamination by physiological frequency artifacts; this practice is so common that this frequency filter range is set by default in popular connectivity processing pipelines (Ciric et al., 2017; Whitfield-Gabrieli & Nieto-Castanon, 2012). Historically, this has not been a particular encumbrance due to the slow nature of BOLD response and the long TRs needed for full brain coverage. Yet consequentially, the spectral content of the fMRI BOLD signal has remained underaddressed.

There is mounting evidence, however, that higher frequencies in spontaneous BOLD fMRI activity may be meaningful. Using advanced accelerated fMRI techniques, canonical resting-state networks (RSNs) have been reported with components at higher frequencies up to 1.5 Hz (Chen & Glover, 2015; Gohel & Biswal, 2015; Kalcher et al., 2014). Further work has shown that BOLD frequency bands throughout the range 0.01–0.25 Hz were highly reproducible and had meaningfully varying network topology across RSNs (Yuen et al., 2019). While caution is still needed regarding high-frequency BOLD noise contamination (Raitamaa et al., 2021), denoising techniques that remove physiological noise sources demonstrate impressive and reproducible results (De Blasi et al., 2020; Pruim et al., 2015), thus opening the door to meaningful high-frequency BOLD analysis. It is already evident that certain frequency-specific measures, such as the fraction of low-frequency power, are useful descriptors of diseases such as schizophrenia (Fryer et al., 2015). Phase-sensitive coherence was previously proposed as a measure of FC, convincingly demonstrating the value of frequency-resolved FC (Curtis et al., 2005; Sun et al., 2004).

Although generative models of fMRI are available via dynamic causal models (DCMs), which infer the coefficient matrix of an underlying vector autoregression (VAR) process (Friston et al., 2014), they were historically limited to correlation-based FC and not power spectra. Recent advances in spectral dynamic causal modeling have expanded their ability to parameterize fMRI coherence spectra with shape and scale parameters (Razi & Friston, 2016; Razi et al., 2017). This is an important advance but one that remains phenomenological and unrelated to the underlying structural substrate or biophysical processes.

Second, current methods are to a large extent statistical—they seek to uncover the difference in summary graph theoretic statistics of FC matrices between diagnostic or task groups. A large body of work has evolved to uncover the graph theoretic features of brain FC networks. These sophisticated network analysis methods have proven value in differentiating between disease (Bassett & Bullmore, 2009; He et al., 2008) and cognitive domains (Park & Friston, 2013), and have proved especially helpful in exploring the neural correlates of developmental and psychiatric disorders (Fornito et al., 2015). Such methods are indeed highly valuable from a practical point of view, but do not typically produce new insight about the underlying biological or biophysical processes that give rise to the observed fMRI data. Despite their broad success, network methods by design remain limited to statistical rather than biological descriptors of brain activity. Indeed, the broad success of graph analysis methods suggests that such a biophysical interpretability is not necessary for many application areas. Yet, the ability to connect the underlying biophysics with observed FC is not only possible but may also lead to enhanced understanding of brain phenomena.

The availability of frequency-rich fMRI further enhances the opportunity to achieve physiologically relevant analysis, as higher frequencies may arise from a combination of the hemodynamic response, neural and cortical time constants, and the aggregate behavior of the brain’s slow oscillatory dynamics and global coupling. Hence there is an opportunity to explore biophysically grounded models of how fMRI time series and its spectral content arise from the brain’s structural substrate. Previous biophysical models of neural activity on coupled systems (Deco, Kringelbach, et al., 2017; Ritter et al., 2013) have been very successful in capturing empirical fMRI signal’s FC but were not designed to capture its spectral content—see Section 5.

Arguably, both opportunities identified above may be addressable via biophysically realistic yet appropriately parsimonious signal generation models that can interrogate BOLD frequency content and its frequency-dependent FC simultaneously. In this article, we propose a different approach of analyzing fMRI data via a powerful yet simple, parsimonious, and linear signal generation model arising from the underlying structural substrate. We show that a full mathematical exposition of this model-based approach can parsimoniously and elegantly exploit both opportunities listed above. Remarkably, this generative model is solved directly in frequency domain in closed form as a summation over graph eigen spectra. Hence, we call this model a spectral graph model (SGM) for fMRI. It can capture both spectral and second-order covariance structures (i.e., frequency-resolved cross-spectral density) simultaneously. We implemented model fitting on two large cohorts of healthy subjects via the optimization of model parameters on individual subjects.

Our approach obviates the need for large time-consuming nonlinear simulations and their attendant burden of parameter inference. The fitted model thus constitutes not only a mapping between structure and function in the brain, but also gives a novel way of analyzing spectrally rich fMRI data. Furthermore, the biologically meaningful global parameters may offer crucial and individualized context to how the subject’s fMRI signals, their spectra, and network topology are related to the underlying structural wiring substrate. We expect this model-based analysis tool will aid practitioners in analyzing fMRI data in a novel, frequency-resolved environment. An analysis overview is depicted in Figure 1.

Notation. In our notation, vectors and matrices are represented in bold, and scalars by normal font. We define a vector of ones as $1$⁠. The Fourier transform of a signal $x(t)$ is denoted as $ℱ{X}(ω)$ for angular frequency $ω=2πf$⁠, where $f$ is the frequency in Hertz. The structural connectivity matrix is denoted by $C={cl,m}$⁠, consisting of connection strength $cl,m$ between any two brain regions $l$ and $m$⁠.

For an undirected, weighted graph representation of the structural network $C={cl,m}$⁠, we model the average BOLD fMRI activation signal for the $l$-th region as $xl(t)$⁠:

whereby the fMRI signal at the $l$-th region is controlled by a characteristic time constant $τ$ (in seconds), and $α$ is a global coupling constant (unitless). Signals coming from regions $m$ connected to region $l$ are scaled by their connection strength $cl,m$⁠. The symbol $⋆$ stands for convolution, $pl(t)$ is the noise or driving function at node $l$⁠. These are global parameters and are the same for every region $l$⁠.

The above model then is the signal generation equation for fMRI. Below we describe each component of this model.

The function $f(t)$ represents the ensemble average neural impulse response, which accounts for various delays introduced by neural processes including synaptic capacitance, dendritic arbors, axonal conductance, and other local oscillatory processes that involve detailed interactions between excitatory and inhibitory neurons, etc. This impulse response is modeled by a Gamma-shape function with a single characteristic time constant, given by $τ$⁠:

This impulse response may also be considered to incorporate, in addition to neural processes, the neurovascular coupling represented conventionally by the hemodynamic response function. Note that while previous analogous rate models for MEG frequencies explicitly introduce axonal conduction delays between brain regions (Cabral et al., 2014; Nakagawa et al., 2014), here we instead incorporate those delays implicitly within $f(t)$ via its single parameter $τ$ rather than explicitly, since our specific focus is to explore much longer time lags in functional activity that arise from other sources. In the fMRI regime, these various processes are not possible to separate or identify, hence a combined model parameter via $τ$ is reasonable.

For convenience and parsimony, we use the same time constant $τ$ to also capture the self-decay, or the network diffusion, term in the signal model (first term on the right hand side; see e.g., Abdelnour et al., 2014). This choice is deliberate, since the self-decay is intimately related to the local impulse response.

We note that it would be straightforward to allow two different time constants, one for the self-decay term and another for the impulse response. The effect of this choice was empirically explored in Section 4.

The global parameter $α$ acts as a controller of weights given to long-range white-matter connections. In the related neural modeling literature, this parameter has been referred to as the coupling constant (Breakspear et al., 2003; Honey et al., 2009; Xie et al., 2021). This parameter represents aspects of global integration and segregation, and may be mediated by attention, neuromodulatory systems, and corticothalamic control signals. Again, these processes are not possible to be separately identified in the fMRI regime, hence lumping them makes sense.

Let us define the Laplacian matrix

where $I$ is the identity matrix and $C$ is the connectivity matrix as defined above. Then the above pair-wise equation readily generalized to the entire brain network, yielding a vectorized differential equation

We note first that Equation (3) is a generalization of the passive diffusion model of brain activity spread—an idea that has been successfully applied to predict steady-state zero-lag FC from the SC’s Laplacian (Abdelnour et al., 2014). Indeed, removing the impulse response $f(t)$⁠, the above system admits a closed-form solution of the free-state evolution of the expected signal with initial condition $x0$ and no external driving signal $p$⁠:

Following the framework well described in Abdelnour et al. (2014, 2015), the above signal equation readily yields the theoretical functional connectome from the Laplacian $ℒ$ of the structural connectome as follows. Assume that only a single region $i$ is experiencing activity at $t=0$⁠, hence $x0=ei$⁠, a unit vector with zeros everywhere except at the $i$-th location. Concatenating this for all regional activations in turn, we obtain:

for some constant $tcrit$ that may be empirically inferred. In our prior work (Abdelnour et al., 2018, 2021), this model was generalized to the “eigen-decomposition” model that relates the SC and FC matrices’ eigenvalues and eigenvectors, and further generalized to a series expansion of eigenvalues (Atasoy et al., 2016; Becker et al., 2018; Deslauriers-Gauthier et al., 2020; Meier et al., 2016; Tewarie et al., 2020).

Here we wish to introduce a richer signal model with local impulse responses, and impart the model with the ability to manifest rich frequency dependencies. Hence we propose the following linear systems approach that relies on Fourier transform of the signal equation. Since the above equations are linear, we can obtain a closed-form solution in the Fourier domain as demonstrated below.

Due to its Gamma shape, the ensemble average neural response function $f(t)$ admits a simple closed-form Fourier transform: $ℱ(f(t))=F(ω)=1τ2(jω+1τ)2$⁠. Taking Fourier transform of the vectorized signal Equation (3), using the notation $ℱ(x(t))=X(ω)$ and $ℱ(p(t))=P(ω)$⁠, we obtain:

Here we used the Fourier pair $ℱ(dx(t)dt)=jωX(ω)$⁠. The above equation can be rearranged as:

To evaluate this expression, we employ the eigen decomposition of the Laplacian:

where $U={uk,∀k∈[1,N]}$ are the eigenvectors and $Λ(α)=diag([λ1(α),…,λN(α)])$ contains on its diagonal the eigenvalues for a given coupling constant $α$⁠. We note that due to the definition of the Laplacian above, the eigenvectors do not have any $α$ dependence, while the eigenvalues do. We also note that the Laplacian above is assumed to be Hermitian, that is, the structural connectivity is symmetric. This is realistic with in vivo imaging in humans but not necessarily in animal models, where directed connectomes are available.

Now we exploit the orthonormality of eigenvectors, that is, $U−1=UH$ and $I=UUH$⁠, to obtain $(jωI+1τF(ω)ℒ(α))−1$$=U(jωI+1τF(ω)Λ(α))−1UH$⁠. Then the above signal equation can be rewritten in closed form as the summation over the eigenvectors:

Equation (7) is the closed-form steady-state solution of the macroscopic signals at a specific angular frequency $ω$⁠. Henceforth we drop the explicit dependency on $α$ whenever convenient. When this equation is required to be explicitly evaluated, we assume that $P(ω)=1$⁠, a vector of ones. This corresponds to either uniform stimulation or perturbation with white noise.

The above model specification requires only two parameters $α$ and $τ$⁠, the latter time constant parametrizes both the Gamma-shaped cortical response function and the network spread process. Yet there was no special reason why the two processes could not be characterized by different time constants. Hence we also evaluated an alternative, less parsimonious model by splitting the singular $τ$ into $τ1$ for network spread and $τ2$ for the cortical response time:

Equation (8) was empirically assessed in comparison with the two-parameter model Equation (7) in Section 4.

With a signal generation equation in both time and frequency domains, it is now possible to explicitly write the structure–function relationship in terms of the eigen-decomposition of the structural Laplacian $ℒ$⁠. There are several equivalent ways to achieve this; here we use the most intuitive approach. Since the cross-spectral density (CSD) is given by $CX(ω)=ℰ(X(ω)XH(ω))$⁠, we get, under the assumption that $ℰ(P(ω)PH(ω))=σ2I$⁠:

In this formulation the eigenvectors of the predicted FC are the same as those of the structural Laplacian (i.e., $U$⁠), while the eigenvalues are related via the diagonal matrix of new (frequency-dependent) eigenvalues $Γk,k(ω|τ,α)=γk(ω|τ,α)$⁠. In this manner, we have reduced the full cross-spectral density of fMRI to modeling just the diagonal eigenvalues of the structural connectome; all region-pair coherences are thus expected to be captured entirely by the eigenvectors $U$⁠. The explicit relationship between the theoretical FC’s eigenvalues and those of the structural Laplacian’s is given by

In this paper we will fit for both the theoretical signal spectrum (Eq. 7) over all frequencies, and theoretical FC derived from the cross-spectral density (Eq. 9).

Data were collected as part of a multisite study aimed at better understanding the brain mechanisms underlying psychosis development and provided by our collaborators in the Brain Imaging and EEG Laboratory at the San Francisco VA Medical Center. Scanning was completed at the UCSF Neuroimaging Center using a Siemens 3T TIM TRIO with the following parameters.

High-resolution structural T1-weighted images: MPRAGE, repetition time (TR) = 2,300 ms, echo time (TE) = 2.95 ms, flip angle = 9 degrees, field of view (FOV) = 256 x 256, and slice thickness = 1.20 mm.

Resting fMRI: T2*-weighted AC-PC aligned echo planar imaging (EPI) sequence: TR = 2,000 ms, TE = 29 ms, flip angle = 75, FOV = 240 x 240, slice thickness = 3.5 mm, acquisition time = 6:22.

Diffusion-weighted MRI: b = 800 s/mm2, 30 diffusion sampling directions, TR = 9,000 ms, TE = 84 ms, FOV = 256 x 256, and slice thickness = 2.00 mm.

For the purpose of this study, we included only healthy subjects from the UCSF study who had both fMRI and DWI data, yielding 56 subjects (20 women; $23.8±8.3$ years).

The T1-weighted (T1w) images and T2*-weighted BOLD images were preprocessed using default procedures in fMRIPrep (Esteban et al., 2019), which is based on Nipype (Gorgolewski et al., 2011) and Nilearn (Abraham et al., 2014). For more details on this pipeline, see the section corresponding to workflows in fMRIPrep’s documentation. For completeness, we summarize the anatomical and functional preprocessing below.

T1-weighted (T1w) images were corrected for nonuniformity (Avants et al., 2008; Tustison et al., 2010) and were used as a reference image throughout the workflow. The T1w reference was skull stripped and segmented based on cerebrospinal fluid, white matter, and gray matter (Avants et al., 2008; Zhang et al., 2001). The T1w reference and template were spatially normalized to a standard space (MNI152NLin2009cAsym) with nonlinear registration (Avants et al., 2008).

FMRIPrep generated a BOLD reference image that was coregistered to the T1w reference using flirt and boundary-based registration (Greve & Fischl, 2009; Jenkinson & Smith, 2001). Head-motion parameters were estimated before spatiotemporal filtering using FSL’s mcflirt (Jenkinson et al., 2002). BOLD runs were slice-time corrected and resampled onto their original, native space by applying transforms to correct for head motion (Cox & Hyde, 1997). The BOLD time series were resampled into the same standard space as the T1w image (MNI152NLin2009cAsym). ICA-AROMA performed automatic detection of signal noise components, including motion artifacts, and saved for use during functional network generation (Pruim et al., 2015).

Average functional time series were extracted from 86 regions of interest (68 cortical, 18 subcortical) as defined by the Desikan–Killiany atlas (Desikan et al., 2006). Functional connectivity processing followed the ICA-AROMA with global signal regression (GSR) pipeline as described in a benchmarking paper (Ciric et al., 2017). This pipeline included smoothing with a 6 mm FWHM SUSAN kernel and regional time series bandpass filtering with a Butterworth filter from 0.01 to 0.25 Hz.

Entries of (zero-lag) functional connectivity matrices were defined as the Pearson’s correlation between regional time series. We also obtained the full cross-spectral density (CSD) of the time series, denoted by the frequency-resolved matrix $CS(ω)$⁠. Functional connectivity matrices were thresholded at the percolation threshold, which reduces the influence of noise in the network and maximizes the networks’ information relative to null models that preserve the strength distribution, degree distribution, and total weight (Bordier et al., 2017; Nicolini et al., 2020). Because the percolation threshold seeks to maximize the network’s information in each FC, every individual had a subject-specific percolation threshold applied to their FC, which was used for model fitting.

Raw diffusion MRI data were processed with T1-weighted (T1w) images using Nipype (Gorgolewski et al., 2011), which implemented functions from FSL, MRtrix3, and ANTS. The raw DWI and T1w images were reoriented and registered to MNI space. T1w image brain extraction was performed with FSL BET (Jenkinson et al., 2012). DWI were denoised with MRtrix3 ringing removal, DWI bias correction, and Rician noise correction (Tournier et al., 2019). Fractional anisotropy in each DWI voxel was estimated from a fit tensor. Streamlines were probabilistically generated with whole-brain seeding using DiPy and an NVIDIA GPU (Garyfallidis et al., 2014). Spherical deconvolution informed filtering of tractograms (SIFT-2) determined the cross-sectional area multiplier for each streamline such that the streamline densities in each voxel are close to the fiber density estimated using spherical deconvolution (Smith et al., 2015). The Desikan–Killiany atlas with 86 regions was linearly and nonlinearly registered to DWI space using ANTS with GenericLabel interpolation, which parcellated streamlines into region-to-region structural connectivity (Avants et al., 2008). Structural connectivity between regions was quantified as the probabilistic and SIFT-2 weighted number of streamlines between regions. Additionally, following Cummings et al. (2022), we added to this structural connectome small amounts of interhemispheric connections between laterally homologous structures, $≈5%$ of the typical connection weight, in order to recover some of the interhemispheric connections missing from DTI tractography outputs. Refer to Cummings et al. (2022) for detailed analysis and justification. Finally, structural connectome of each subject was normalized by the sum of all matrix entries, such that individual connections have arbitrary units that add up to 1 for all subjects.

We sought to evaluate the model’s robustness and capacity to analyze fMRI data with richer frequency content. To this end, we used the recently published publicly available dataset for Microstructure-Informed Connectomics (MICA-MICs) with 50 healthy human subjects (⁠$23$ women; $29.54±5.62$ years) (Royer et al., 2022). These data use similar pipelines to those described above (see Cruces et al., 2022 for processing details), but importantly, the fMRI data were collected with a 600-ms TR, meaning that the time series contains much richer frequency content. Since these data are without global signal regression (GSR), we performed GSR by extracting and removing the signal’s first principal component from all regions.

While the theoretical model is simple and straightforward to evaluate simultaneously on FC and spectra, several implementation details were found to improve the quality and reliability of the fits to individual subjects. As described below, not all eigenmodes in the theoretical model are equally important, and not all portions of the signal cross-spectral density are equally useful.

We wish to achieve match between the theoretical cross-spectral density matrix $CX(ω|τ,α)$ in Equation (9) and the corresponding empirical CSD denoted by $CS(ω)$⁠. However, a direct fitting to the 3D CSD array (⁠$regions×regions×ω$⁠) was found to be problematic: as depicted in Figure 2, both the empirical and theoretical CSD are extremely sparse matrices, with near-zero values at frequencies higher than around 0.10 Hz and especially in off-diagonal entries, but with a significant amount of noise coming from both the measured signal and highly challenging Fourier transform-related spectral estimation issues. Thus fitting to the entire CSD volume is extremely poorly posed. Therefore, we developed an alternative strategy whereby we identified the most signal-rich and reliable portions of the CSD volume—the zero-lag FC (given by the integration over all frequencies of CSD: $∫CX(ω)dω)$⁠; or a single slice of CSD $CX(ω0)$ at a fixed frequency $ω0$⁠, the frequency where the fMRI spectral response was assessed to be the maximum for a given subject. This approach allows the model fits to be informed by the most critical elements of interest to the modeler, without necessarily suffering from the challenges explained above.

In a similar vein, apart from the formal definition of PSD as the diagonal of the CSD, that is, $diag(CX(ω))$⁠, it is also useful to define simply the power of the signal under an explicit driving signal, that is, as $|X(ω)|$⁠, when $P(ω)=1$⁠.

Not all eigenmodes in the summation shown in Equation (7) are equally important or relevant for predicting FC. As others have noted, the SC matrix yields many eigenmodes that either do not participate in dynamic FC (Preti et al., 2017; Preti & Van De Ville, 2019) or do not contribute to models fitted to fMRI (Abdelnour et al., 2018; Cummings et al., 2022) or MEG data (Verma et al., 2022). Hence, we developed a simple technique for evaluating eigenmode weighting by projecting the FC into the Laplacian’s eigen space, thus performing a graph Fourier transform of the FC matrix: $Q=UTFC U$⁠. The diagonal entries of this matrix correspond to the strength with which each respective eigenmode participates in FC. Hence, we define for each eigenmode $k$ its “graph Fourier weight” (GFW) as $GFWk=|Qk,k|$⁠. Although this calculation requires the empirical FC matrix, it is not necessary to use individual FCs, which would be impractical and circular. Instead, we found that the weights are reliably stable between individuals, hence they were precomputed using average FC and SC matrices.

The model parameters $τ,α$ are optimized per subject to maximize two objectives simultaneously: the spectral correlation and the FC correlation, as defined below. In this manner, we can obtain a unified model that reproduces both the regional power spectra and the second-order covariance statistics given by the frequency-resolved FC matrix (i.e., cross-spectral density). The cost function we used is, therefore, defined as

Functional correlation was defined as the Pearson’s correlation between the theoretical and empirical FC (only the upper triangular portion thereof), while spectral correlation is defined analogously on the PSD. We implemented a constrained cost function minimization, available as the routine fmincon() in MATLAB version R2021b. The parameters $τ,α$ were given lower limits 0 (to ensure positive values). We used default options for fmincon(), with $50$ iterations and tolerance of $10−6$ as stopping criterion.

To compare our proposed model with random structural networks, we constructed a null distribution by randomly sampling our dataset $1,000$ times and generating a randomized SC network using the Brain Connectivity Toolbox function randmio_und_connected() (Rubinov & Sporns, 2010). This function preserves the network’s degree distribution while ensuring the graph remains fully connected, which is highly significant for the Laplacian eigenmodes. To compare our proposed model with random functional networks, we again constructed a $1,000$ sampled null distribution by randomly permuting the BOLD time series region ordering, hence randomizing the FC matrix. We use the randomized SC to predict empirical FC, and we use the empirical SC to predict the randomized FC and spectra ordering. We compare all the resulting distributions by two-tailed t-tests.

There are no current model that fits to and predicts the spectral content of fMRI, with the exception of spectral DCM (Friston et al., 2014; Razi et al., 2017). However, numerous models are available that predict FC, that is, the covariance structure of the signal at zero lag. We compared the FC portion of our results against archetypal model-based approaches for the prediction of FC that explicitly involve the structural connectome (SC). Therefore, we excluded spectral DCM method since they do not involve SC, but included both linear and nonlinear SC-based graph models. We also excluded recent matrix algebraic methods such as Becker et al. (2018), on the grounds that it uses, and needs to infer, a large number of model parameters, for example, a full $N×N$ rotation matrix (Becker et al., 2018; Deslauriers-Gauthier et al., 2020). Such models do not concord with the attribute of parsimony that is essential to our approach.

Due to their emerging popularity and conceptual similarity, algebraic graph models (i.e., those that involve a transformation of SC eigenvalues) are highly pertinent benchmarks for our current proposal. Although many such methods are now available (Deslauriers-Gauthier et al., 2020), we chose the two models which are both the earliest and the most parsimonious: the network diffusion model (Abdelnour et al., 2014) and Eigen-decomposition model (Abdelnour et al., 2018). Both models have similar parsimony to ours: one parameter and three parameters, respectively.

To compare the proposed SGM with the more commonly reported generative simulations involving nonlinear neural dynamics, we implemented the popular “Conductance based NMM” described by Breakspear et al. (2003; Honey et al., 2009) using the same model parameters previously used in Abdelnour et al. (2014, 2018). We refer to these original publications for modeling details; here we note the salient points from their paper. This NMM implements the nonlinear Wilson–Cowan system of voltage and conductance currents (Wilson & Cowan, 1972, 1973). The dynamical variables incorporated are the mean membrane potential of pyramidal cells and of inhibitory interneurons, and the average number of “open” potassium ion channels. The mean cell membrane potential of the pyramidal cells is governed by the conductance of sodium, potassium, and calcium ion through voltage- and ligand-gated membrane channels. The firing of excitatory and inhibitory cell populations feeds back onto the ensemble through synaptic coupling to open ligand-gated channels and raises or lowers the mean membrane potential accordingly. In the case of excitatory-to-inhibitory and inhibitory-to-excitatory connections, this is modeled as additional input to the flow of ions across the membrane channel, weighted by functional synaptic factors. Coupling between brain regions (i.e., graph nodes) is introduced as competitive agonist excitatory action at the same populations of NMDA and AMPA receptors.

Neuronal population dynamics were simulated at 0.2-ms resolution for 16 min using a system of 86 neural masses, each corresponding to a Desikan parcel. Each neural mass represents a population of densely interconnected excitatory and inhibitory neurons, in which the effects of both ligand- and voltage-gated membrane channels are accounted for. Node-level neural masses were coupled to each other via structural connectome. In our implementation, we chose not to resample raw fiber strengths of the connectome onto a Gaussian distribution, in order to keep the comparison fair. Resampling is a convenience, but is not biophysically justifiable. These connection strengths between region-pairs are multiplied by a single free global “coupling strength” parameter $c$ that modulates the inter-regional coupling, and we evaluate the model over a range of $c$ (⁠${0.05,0.1,0.15,0.2,0.25,0.3,0.35}$⁠) and compare this model with SGM for fMRI using the optimal $c$⁠. All other NMM parameters were kept at their optimal default values, as determined in the original publications.

To estimate BOLD signals for each neural mass, a Balloon–Windkessel hemodynamic model was convolved with the simulated 0.2 ms NMM time series. The resulting low-frequency time series data were detrended and a global mean signal was regressed out from each node’s time series. Finally we computed the implied FC, that is, the correlation between any two regions’ simulated-BOLD time series. This FC was then correlated against empirical FC of each subject to obtain our measures of model performance—Pearson R, geodesic distance, and mean square error.

Our results are organized as follows. First, we establish that there is indeed interesting frequency content in rs-fMRI data, even in conventional recordings obtained at 2 sec TR, confirming prior work on coherence FC (Curtis et al., 2005; Sun et al., 2004). We then demonstrate the power of the proposed model-based analysis of frequency-resolved fMRI by fitting the proposed SGM to empirical data. Our hypothesis is that even a global, spatially invariant model with only two biophysical parameters should be able to predict conventional zero-lag FC and power spectra. We show this first on slower, conventional acquisitions obtained at our institution, since the vast majority of fMRI data are acquired this way. Of course, the true power of the approach would be better leveraged when fitting to fast acquisitions with a wider frequency range—this is our final result, and it is demonstrated on a high-quality, public fMRI dataset consisting of 50 healthy subjects.

We first show that fMRI signal has frequency-dependent information content worth modeling. Figure 2A shows a set of cross-correlation signals $Rij(τ)$ over a range of $τ$⁠. The cross-correlation function is symmetric around $0$⁠, hence only the positive lags are shown. It is clear that while the maximal number of correlations peaked at near-zero lags, many region pairs gave maximal correlations at nonzero lags. Interestingly, region pairs whose zero-lag correlations are small or negative can sometimes produce high correlations at nonzero lags. Figure 2A also contains a histogram of the lags at which maximum correlation was achieved. The range of lags is widely distributed, and while $τ=0$ remains the single most frequent lag, more region pairs overall peak at nonzero lags. Figure 2B shows the $86×86×Nfreqs$ frequency-resolved cross-spectral density (CSD) array: this corresponds to the spectral power of the cross-correlation function. Due to its evident frequency-selective distribution, it contains information beyond the conventional zero-lag FC, which corresponds to the integral over all frequencies of the CSD array. Since the full CSD volume contains large components with little signal, we carefully chose the most information-rich features from the full CSD: the power spectral density (PSD, defined as the diagonal of the CSD volume) and CSD at a specific frequency $ω0$ at which the pairwise sum of CSD has the highest value for a given subject. These features are illustrated in the figure.

The vast majority of fMRI scans, like the previous results, have $TR≥2$ sec, hence we first demonstrate that the proposed analysis technique can achieve strong fits to current established fMRI acquisitions. For this purpose, we used local fMRI studies performed at our institution—see Section 3. Using SGM induced by the structural Laplacian matrix, we find that the model performs well, fitting simultaneously to zero-lag FC and BOLD fMRI spectra with mean (std) Pearson’s $R$ = 0.59 (0.08) for FC and Pearson’s $R$ = 0.70 (0.08) for BOLD spectra (Fig. 3A). We additionally report optimal parameter values mean (std) for $α$ = 0.80 (0.09) and $τ$ = 1.96 sec (0.56 sec) (Fig. 3B). In Figure 3C, we show three subjects’ example SC, empirical FC, and SGM-predicted FC, along with their individual Pearson’s $R$ similarity; in Figure 3D, we show the same three subjects’ empirical and predicted BOLD spectra also with its associated Pearson’s $R$ similarity. From these results, SGM clearly reproduces salient elements of both empirical FC and BOLD spectra, using a single model parametrization.

Although very few sites routinely acquire low-TR (i.e., fast) fMRI data currently, we were able to obtain the public MICA dataset with $TR≈0.6$ sec to demonstrate the feasibility of fitting the SGM to both FC and a richer frequency content than conventional acquisitions currently allow. As shown in Supplemental Figure 1, SGM for fMRI performed similarly on this new frequency-rich data—mean (std) $R=0.57(0.08)$ for FC; $R=0.87(0.04)$ for spectra—compared with our previous result. This not only supports the model but also serves as independent confirmatory evidence analogous to the main results in Figure 3.

Since our model was intended to work on frequency-rich data, it is useful to assess whether the choice of TR plays a role in the model fits. Hence, we performed a new experiment in which we repeated our SGM fitting to several downsampled versions of the MICA dataset. We systematically downsampled the TR of the above MICA dataset (Supplemental Fig. 1) from its the natural 600 ms TR to lower TRs in intervals of 200 ms up to 2000 ms, using MATLAB’s resample function. The fitted SGM parameters and their resulting correlations across each downsampled TR are shown in Supplemental Figure 2. It is noted that the quality of the fits did not change appreciably in this range of TRs; possibly due to the fact that our model fits involve both spectra, as well as FC at maximum amplitude, the latter of which typically occurs at lower frequencies. The main effect of downsampling TR was a slight reduction in the global time constant $τ$⁠.

Several algorithmic and practical considerations affected the final technique demonstrated above. In this section, we briefly explore the effect of these choices and provide empirical support for why they were chosen.

We weighted the model’s eigenvalues by GFW, which quantifies how much the eigenmode is expressed in the average FC’s projection matrix $Q$ (Fig. 4A). The group-averaged GFWs across the eigen spectrum are shown in Figure 4C. This reweighting of eigenvalues alters the eigenmode’s theoretical frequency response (Fig. 4B), and substantially improves the model’s overall performance (Fig. 4D). To our knowledge, the exact manner in which we have implemented eigenmode weighting has not appeared previously.

Of course, there are many alternative decisions we could have made regarding individual versus group-averaged eigenmode weights, as well as how to threshold the empirical FCs using the percolation approach. For completeness, Table 1 summarizes results for these different conditions. All results contained elsewhere in this paper comprised the best condition set from this table, highlighted in bold.

The covariance and Fourier formalism given in the Theory section accommodates multiple definitions of both what constitutes FC and what constitutes PSD, depending on the underlying assumptions. We have thoroughly evaluated the effect of these choices, summarized in Table 2. In these explorations, we used the best condition from Table 1, using individual thresholds for FC and group-mean eigenmode weighting.

The main results of this paper shown in Figure 3 and Supplemental Figure 1 employed the maximum signal slice of the CSD volume for FC, that is, $CX(ω0)$⁠. However, the formal definition of zero-lag FC, given by $∫CX(ω)dω$⁠, also gives a very good match to empirical data, with mean (std) Pearson’s $R=0.54(0.10)$ for FC and Pearson’s $R=0.69(0.08)$ for BOLD spectra. This is shown in the first two rows of Table 2. We did not have an a priori preference of zero-lag FC over specific FC evaluated at $ω0$ - —but from a practical point of view, the latter was found to be more reliable to fit and gave slightly better results. In a similar vein, both the formal PSD, as the diagonal of the CSD, that is, $diag(CX(ω))$⁠, and the direct signal power under an explicit driving signal, $|X(ω)|$⁠, when $P(ω)=1$⁠, produced good results. Table 2 shows that the direct PSD $|X(ω)|$ in fact does slightly but consistently better than $diag(CX(ω))$ at fitting to empirical BOLD power spectra.

As noted in the Section 2, our model specification requires only two parameters: $α$ and $τ$⁠, the latter being shared by both the Gamma-shaped cortical response function and by the network spread process. Our motivation was essentially parsimony; yet there was no special reason why the two processes could not be characterized by different time constants. To empirically support our choice, we performed an empirical comparison between the two versions. We split the singular time constant $τ$ of the main model in Equation (7) into two: $τ1$⁠, the network diffusion time, and $τ2$⁠, the cortical response time—contained in the alternative model Equation (8). The results are shown in Supplementary Figure 3 and compared against the original two-parameter model. For these comparisons, we initialized both time parameters to the same value and gave them an identical range. As it transpired, the two constants meaningfully diverged from the singular version, yet the end model was no better than the original, more parsimonious model. We deemed the net benefit of adding the new parameter to be negative; henceforth all results are presented for the two-parameter model.

We sought to assess how SGM for fMRI performed in comparison with previous linear SC-FC models implemented on the same dataset. Since existing models can only predict FC and not the PSD of fMRI signal, we used only the FC component of our predictions in these comparisons. We observed that SGM obtains better fits to individual FC compared with both previous linear graph models, the network diffusion model (Abdelnour et al., 2014) and the eigen-decomposition model (Abdelnour et al., 2018; Fig. 5). Then we assessed the performance of the connectome-coupled neural mass model (Breakspear et al., 2003; Honey et al., 2009), whose Pearson’s R was centered around 0.44—lower than the above linear models and far lower than the proposed SGM. Using a one-way analysis of variance (ANOVA), we found that there was a significant group-wise difference (⁠$p<<0.01$⁠) between the SGM and all other benchmark models. Importantly, SGM showed better performance than all alternative models.

To ensure that SGM for fMRI goodness of fits are specific to the structural network, we sought to test whether randomizing the SC matrix diminished SGM’s capacity to fit to empirical FC and BOLD spectra. Indeed, we found that similarity distributions using randomized samples of subject-specific SC were significantly different than the empirical SC fit to FC (⁠$p=0$⁠) and to BOLD spectra (⁠$p=1.17×10−13$⁠) (Fig. 6A, B).

We also sought to evaluate whether SGM for fMRI could predict FC derived from a time series obtained by randomly reordering regions (⁠$Trand$⁠) while maintaining the empirical SC. We found that randomly shuffing the ROI time series again destroyed SGM’s ability to predict FC (⁠$p=0$⁠) and worsened the prediction of BOLD spectra (⁠$p=2.84×10−7$⁠) (Fig. 6C, D).

The results from the null models require careful assessment. The predicted FC from a randomized SC is clearly separated from the predictions using empirical SC, yet the spectral predictions for randomized versus empirical SC are significantly different but not as well separated. It is likely that SC randomization destroys the spatial patterns and region order inherent in FC, while the spectral content is less affected since in theory the random graph can encompass eigenmodes and eigenvectors that may be nearly sufficient to produce a desired set of spectra. This is observed even more clearly in the second randomization, that of time series regional order. Here again, the implied FC is destroyed due to region reordering, but the spectra are preserved (albeit with a random reassignment to regions). This implies again that the predicted model’s spectra do not possess strong regional specificity.

There is a need for innovation in new-generation analysis tools for rs-fMRI focused on overcoming current gaps. In particular, we identified (a) the opportunity to use rich spectral content available in modern multiband fast acquisition techniques and (b) the opportunity to relate fMRI analysis products to biophysically relevant parameters, which might one day become useful as biomarkers of task, cognitive or disease states. The proposed spectral graph model uses the brain’s structural network to generate a synthetic fMRI signal; in this study it is implemented as a general-purpose yet biophysics-informed tool to analyze wideband fMRI and its induced FC simultaneously. It does so with the simplicity and parsimony of spectral graph theory, using only two global parameters that between them accommodate the variability inherent within healthy subjects’ FC and spectra: $α$⁠, which accounts for global coupling, and $τ$⁠, which accounts for signal propagation speed. The model naturally decomposes into a sum over the eigenmodes of the structural graph Laplacian. Remarkably, the model admits closed-form solutions of the regional signal’s frequency spectra as well as frequency-dependent FC—a useful departure from existing coupled neural systems that require lengthy numerical simulations. A key refinement of the approach involved a priori weighting of the Laplacian eigenmodes such that only those eigenmodes that contribute to the resting-state fMRI signal patterns in the entire cohort were retained in the summation (Fig. 4).

We reported excellent capacity of the model to predict empirical power spectra as well as FC, on both conventional acquisitions (Fig. 3) and fast multiband acquisitions (Supplemental Fig. 1). We found that frequency-resolved SGM produced significantly improved predictions of FC over previous linear and nonlinear SC-FC models, significantly exceeding them in both correlation and geodesic measures of performance (5). Most strikingly, SGM gave far better fits than a coupled system of nonlinear neural masses, supporting the idea that macroscopic neurophysiological data on a graph can be sufficiently modeled with linear metrics, and nonlinear methods may not be required for problems of such scale (Hartman et al., 2011; Hlinka et al., 2011). We demonstrated that SGM’s predictions are specific to the SC network topology by comparing with a null distribution of randomized SC matrices. Our results were far better than could be expected by chance (Fig. 6). In the following sections, we discuss the implications of this spectral graph model and contextualize these results against the current literature.

Our model was able to access most observable configurations of spatial and spectral patterns seen in real rs-fMRI data by tuning only two global and biophysically relevant parameters: coupling strength and neural time constant. That such a parsimonious and analytical model can produce a frequency-rich repertoire of activity may be a novel contribution to the field. As observed from Figure 3 and Supplemental Figure 1, parameter regimes typically recruit and “steer” multiple eigenmodes in such a manner that a small number of them can reproduce signal spectrum and FC at any relevant frequency. Possibly, this points to an essential characteristic of real brain activity, which is thought to accommodate a large repertoire of microstates and their concomitant spatial patterns. Available literature on wide-band frequency brain recordings demonstrates that different fMRI functional networks are preferentially encapsulated by different higher-frequency bands via phase and amplitude coupling (Deco et al., 2009; Ghosh et al., 2008a). Phase and amplitude coupling of oscillatory processes in the brain is evidently important for the formation of coherent wide-band frequency profiles of brain recordings and processing of information (Fries, 2005; Deco et al., 2009; Ghosh et al., 2008a; Schnitzler & Gross, 2005; Varela et al., 2001). Historically, such processes could be modeled via large scale nonlinear simulations of oscillatory neural activity (Honey et al., 2009; Jirsa & Haken, 1997; Jirsa et al., 2017). In contrast, our approach does not require high-dimensional model parameters nor extensive simulations. This raises the possibility that complex behavior may be achievable by simple and parsimonious mechanisms.

The present computational study is not intended to explore the neural mechanisms that might control the model’s biophysical parameters. Nevertheless, several possibilities may be considered. Coupling strength $α$ is a direct scaling of white-matter excitatory long range connections between neural populations in the brain. Parameterization of coupling strength between distant brain regions via the connectome is ubiquitous in connectivity-based models of BOLD fMRI (Abdelnour et al., 2014; Deco et al., 2009; Honey et al., 2007, 2009) and electroencephalography activity (David & Friston, 2003; Ghosh et al., 2008b; Spiegler & Jirsa, 2013). Furthermore, pathological FC patterns as a result of disconnections in the brain can be reproduced with decrease in coupling strength (Cabral et al., 2012; Jirsa et al., 2010). The other tunable parameter in our model, neural time constant $τ$⁠, is the aggregate lumped time constant that captures in a single number the overall effect of membrane capacitance, conductance speed, delays introduced by the dendritic arbors of the ensemble of neurons within a region, and the local processing delays of signals within the region. It also incorporates the convolution of the neuronal signal with the hemodynamic response pertaining to the BOLD signal. At this time, it is not possible to separately characterise, from a lumped aggregate measure such as $τ$⁠, the potential underlying neurobiological processes that may contribute to it. While the individual component processes are likely complex and challenging to infer in a model such as this, their overall lumping into a single time constant $τ$ enables us to capture emergent behavior without breaking the essential parsimony of the generative model. Parameter $τ$ may, therefore, still be considered a biologically relevant and interpretable quantity, capturing a characteristic time constant in the brain which is measurable and has actual units (in seconds).

We empirically investigated the possibility that an additional time constant would improve the model further, but this was not the case; see Supplemental Figure 3. Hence, it may be deemed that the two-parameter model is sufficient.

We review two types of prior models relevant to current work: graph models and generative models.

Recent graph models involving eigen spectra of the adjacency or Laplacian matrices of the structural connectome have greatly contributed to our understanding of how the brain’s structural wiring (SC) gives rise to its FC. Such models typically assume that SC and FC are not independent entities, yet their relationship is not direct (Honey et al., 2009). In addition to connection strength between regions, metrics such as anatomical distances (Alexander-Bloch et al., 2013), shortest path lengths (Goñi et al., 2014), diffusion properties (Abdelnour et al., 2018; Kuceyeski et al., 2016), and structural graph degree (Stam et al., 2016) were also found to contribute to the brain’s observed functional patterns. Eigenmodes and eigenvalues of a network’s adjacency or Laplacian matrix express an ordered, low-dimensional, representation of that network’s topology. Eigenmodes are thus conceived to represent spatial frequencies of the underlying network—the fundamental spatial patterns contained within the network’s topology (Abdelnour et al., 2018; Atasoy et al., 2016; Huang et al., 2018). Reassuringly, these SC eigenmodes are reproducible, with the first few lowest “spatial frequency” patterns being similar across different atlas parcellation resolutions, between different human subjects, and between two scans from the same subject (Wang et al., 2017). The success of eigen-mapping methods has been attributed to SC and FC sharing fundamental features in their eigen space (Abdelnour et al., 2018; Liégeois et al., 2020). A recent extension introducing transmission delays and a “complex Laplacian” showed that different sets of eigenmodes contributed in varying degrees to canonical RSNs (Xie et al., 2021). Higher-order walks on graphs, involving a series expansion of the graph matrices, have also been quite successful (Becker et al., 2018; Meier et al., 2016). The diffusion and series expansion methods are themselves closely related (Robinson et al., 2016), and almost all harmonic-based approaches may be interpreted as special cases of each other, as demonstrated elegantly in recent studies (Deslauriers-Gauthier et al., 2020; Tewarie et al., 2020). Yet another mapping between FC and SC eigenvalues via the Gamma function was demonstrated (Cummings et al., 2022).

In the Theory section, we were able to show mathematically that the eigen-mapping model central to these studies is essentially a special case of the current generative model, reducing to the former when the neural response function is removed and the driving signal is i.i.d. Gaussian. The key differentiator, however, is that prior graph models have not considered the spectral content of the signal, and in fact do not have the ability to give nontrivial spectra whatsoever. Our concept of a priori learning and weighting the importance of the eigenmodes that go into the SGM (Fig. 4) is also inspired from prior work. Indeed, taking projections of SC eigenmodes into a functional space is a simple, powerful means to identify eigenmode weightings. A recent study has used the signal’s eigenmode projection as a form of filtering (Preti & Van De Ville, 2019).

While keeping the parsimony of eigen-mapping graph models, our approach adds an explicit signal model with frequency content accommodating oscillatory frequency and phase shifts between brain regions. However, the SGM signal generation model is not the only example of generative modeling in fMRI analysis. In fact, many such generative models exist, under two broad forms: (a) dynamic causal modeling with an unknown effective connectivity matrix (Friston et al., 2014) to fit a multivariate autoregressive process to the time series (Vidaurre, Abeysuriya, et al., 2018; Vidaurre, Hunt, et al., 2018) or complex fMRI cross-spectra (“spectral DCM”) (Friston et al., 2014; Razi et al., 2017) or (b) neural mass model (NMM) simulations, which use coupled nonlinear NMMs, typically simulated first to much higher (up to 40 Hz) frequencies, then downsampled to BOLD frequency ranges via a Hilbert envelope (Cabral et al., 2014; Deco, Cabral, et al., 2017). Each method is briefly described below and compared against the current approach.

While the goal of VAR and DCMs in fMRI analysis as means to make model inferences about FC is similar to our proposal, the two frameworks are vastly different in terms of approach and dimensionality. DCMs examine the second order covariances of brain activity (Park et al., 2018; Razi et al., 2015), and it is only recent works with spectral and regression DCMs that have expanded the model coverage to the whole-brain scale (Frässle et al., 2018, 2017; Razi et al., 2017). Since they are not designed to be biologically grounded in the underlying connectome, DCMs have many more degrees of freedom than our work because they parameterize for different interactions within and between brain regions. Although we did not benchmark against spectral DCM for this reason, we note that SC has recently been incorporated into spectral DCM via empirical priors (Greaves et al., 2024), and that DCM and associated methods are often applied in a manner that is informed by SC (e.g., effective connections are inferred for structurally connected regions only (Upadhyay et al., 2008)). Hence benchmarking against structurally informed DCM would present an attractive opportunity for future work.

An effective way to address higher frequency information is to employ structural connectivity to couple anatomically connected neuronal assemblies (Wilson & Cowan, 1972, 1973). Numerical simulations of such neural mass models (NMMs) provide an approximation of the brain’s local and global activities, and are able to achieve moderate correlation between simulated and empirical FC (Nunez, 1974; Jirsa & Haken, 1997; Honey et al., 2009; Spiegler & Jirsa, 2013; Valdes et al., 1999). These coupled NMMs are perfectly capable of producing high-frequency content, but their existing use in fMRI modeling tends to focus on reproducing zero-lag FC exclusively and not on generating the fMRI signal PSD at nonzero frequencies. Moreover, these models rely on lengthy and expensive time series simulations of local neural masses to derive dynamical behavior, which are then used to generate effective connectivity matrices. Such approximations through stochastic simulations are unable to provide a closed form solution; their inference requires optimizing a large number of model parameters and exacerbates and inherits interpretational challenges. It is noteworthy that SGM for fMRI may be considered a highly simplified and linearized NMM, but with a focus on the coupling between remote neural masses rather than detailed modeling of those masses themselves. Despite only fitting to FC (and not spectra as SGM does), the coupled NMM we implemented here was unable to produce results on FC comparable with ours, as demonstrated in Figure 5E. By avoiding large-scale simulations of neuronal activity, our approach overcomes these practical challenges and aids biological interpretability.

The current theory is also quite applicable to faster acquisitions such as EEG and MEG, since it has the capacity to explicitly capture frequency response via graph spectra. A similar signal generation equation was recently proposed by us for capturing higher frequencies obtained from MEG (Raj et al., 2020), and has been successfully applied (Bernardo et al., 2024; Verma et al., 2022). The MEG context is different (high frequency regime) and that the model involved local excitatory and inhibitory neuronal subpopulations and phase lags arising from axonal conductance—elements which are not pertinent to the fMRI regime. It would be very interesting to explore whether such models can fit to both slow fMRI and fast M/EEG simultaneously. Despite several studies of associations between fMRI and MEG FC and phase amplitude coupling, there is not yet a single convincing model that can simultaneously capture both regimes.

Any analysis tool must eventually find purpose in practical and clinical domains. Since the proposed technique is highly parsimonious and admits closed-form solutions, it presents an unrivalled opportunity: easy and straightforward fitting to real data, over all frequencies of interest. As shown in Figures 3 and 5, and Supplemental Figure S1, the model achieves excellent fits to individual subjects, quantitatively superior to FC predicted by prior linear models, and additionally fits to regional power spectra which no current model can achieve. Our work may, therefore, find applicability in certain clinical and neuroscientific contexts where predicting functional patterns from structure is important (Fornito et al., 2015; Jiang, 2013), particularly in cases of epilepsy (Abdelnour et al., 2021; Coan et al., 2014), stroke (Kuceyeski et al., 2016; Rehme & Grefkes, 2013), schizophrenia (Cummings et al., 2022; Fryer et al., 2015), and neurodegeneration (Zimmermann et al., 2018). However, this aspect will require careful empirical assessment; both the model fitting strategy and the quality of empirical data would need to be further refined before a successful demonstration of this aspect.

Several assumptions and limitations are noteworthy. First, we note that the SGM uses a Hermitian Laplacian, which is invalid in cases where the underlying structural connectome is asymmetric (e.g., in tracer data from animals). In that case Equation (6) would not hold, and the eigenvectors of the predicted FC would not be the same as those of the structural Laplacian. Next, our theorized model relies on a uniform neural time constant throughout the brain. In reality, the amount of myelination and synaptic strength varies greatly in the brain, as do the number of neurons and their local electrophysiological properties. However, this strong assumption did not adversely affect the model’s predictive ability, which is similar or higher than currently reported in the literature. Certainly, predictive power could be improved by the use of regionally varying model parameters, but the current approach has the inestimable benefit of a low dimensional and biophysics-informed model.

While we have attempted to relate the model’s parameters at the aggregate level to various underlying neurobiological mechanisms, at the moment we do not have the ability to clearly distinguish between them. In particular, the parameter $τ$ encompasses multiple time-varying processes that contribute to the fMRI data, but given the frequency content of the time-series signals can only be assessed in aggregate. It will be a subject of future analyses to carefully unpack the aggregate-level model into its constituent elements and subprocesses, which will likely require novel experimental setups and much higher sampling rates than currently available.

A linear model such as SGM is limited in the repertoire of dynamics that it can produce. Interestingly, a recent comparison (Nozari et al., 2020) demonstrated that linear models predicted resting fMRI time series better than nonlinear models. They argued that this may be because of the linearizing affects of macroscopic neuroimaging due to spatial and temporal averaging, observation noise, and high dimensionality, indicating that a linear model may be sufficient to capture the observed fMRI dynamics.

The functional data used in the main study had $TR=2$ sec temporal resolution, which limits the frequencies resolvable to the Nyquist limit of 0.25 Hz. Yet, most fMRI data are collected with 2-sec TR, and conventional resting-state analysis does not explore frequencies above 0.1 Hz, hence the presented temporal resolution is appropriate given our aims. Interestingly, SGM performed similarly on both the conventional 2 sec TR acquisitions (Fig. 3) and the faster 0.6 sec TR data (Supplemental Fig. 1). The quality of the fit did not appear to depend much on the effective TR of the data (Supplemental Fig. 2). This is somewhat contrary to the expectation that a frequency-rich model would fit better to faster data. This warrants further investigation, and potentially opens the door to future improvements in the model itself. One possible explanation is that our model fitting involves both FC and the fMRI spectrum; the former was evaluated only at the peak amplitude frequency, which typically occurs below 0.08 Hz. Due to this aspect, the model fits would naturally favor lower frequencies than higher frequencies. If we assume that high frequencies add more variability to the underlying empirical FC, then we might expect that the performance would potentially be worse for FC when fitting to higher frequencies. As the field moves toward faster fMRI protocols, we expect that far more frequency-rich data would become available for a more thorough testing and validation than is currently possible.

It should be noted that the current framework for obtaining model fits to empirical data involves optimizing model parameters using gradient descent routines which are most effective for convex optimization problems, and may not easily generalize to nonconvex scenarios. Further, the point-optimization method used here does not constitute a full inference procedure, which requires estimation of error bounds (Bzdok & Ioannidis, 2019), and typically uses sampling methods to achieve inference of model posteriors, for example, Markov Chain Monte Carlo (MCMC) (Raftery & Lewis, 1996; Sengupta et al., 2016). Hence future effort will be required to achieve full inference. Due to the computational complexity of the model, it may not be realistic to perform MCMC. In that case, we might leverage recent successes in our laboratory in the use of novel neural network methods such as simulation-based inference (SBI; see e.g., Bernardo et al., 2024; Jin et al., 2023). Further, our cost function is based on Pearson’s correlation between predicted and empirical FC and spectra, which causes our fits to reproduce broad similarities in spectral and FC shape rather than absolute scaling, which currently remains outside the scope of the model.

Finally, while weighting eigenmodes constitutes a clear improvement and novelty, this requires prior access to a group-level SC-FC data. Therefore, our results could be strengthened in future studies by applying population averages from a separate large cohort.

Code for SGM for fMRI has been made available at https://github.com/Raj-Lab-UCSF/SGMforFMRI along with group-averaged SC. The full dataset can be made available upon request.

Ashish Raj: Conceptualization, Methodology, Writing—Original Draft, Writing—Review & Editing, Supervision, Project Administration, Funding Acquisition. Benjamin S. Sipes: Software, Formal Analysis, Investigation, Data Curation, Writing—Original Draft, Writing—Review & Editing, Visualization. Parul Verma: Formal Analysis, Investigation, Writing—Original Draft, Writing—Review & Editing. Daniel H. Mathalon: Resources, Writing—Review & Editing. Bharat Biswal: Writing—Review & Editing. Srikantan Nagarajan: Methodology, Supervision, Writing—Review & Editing.

The studies involving human participants were reviewed and approved by UCSF’s Institutional Review Board. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

This work was supported by NIH grants R01EB022717, RF1AG062196, R01AG072753.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00381.