Structure-function models of temporal, spatial, and spectral characteristics of non-invasive whole brain functional imaging.

We review recent advances in using mathematical models of the relationship between the brain structure and function that capture features of brain dynamics. We argue the need for models that can jointly capture temporal, spatial, and spectral features of brain functional activity. We present recent work on spectral graph theory based models that can accurately capture spectral as well as spatial patterns across multiple frequencies in MEG reconstructions.

1. Introduction

One of the outstanding questions in the field of neuroscience is to understand how dynamic brain functional activity arises despite being constrained by the static anatomical structure (Bassett and Bullmore, 2009; Cao et al., 2014; Fornito et al., 2015; Suarez et al., 2020). To explore this relationship, both structure and functional activity are estimated from non-invasive neuroimaging techniques such as magnetic resonance imaging (MRI), functional MRI (fMRI), and magneto-/electroencephalography (M/EEG). Given that such techniques cannot capture the underlying mechanistic biophysics of the brain activity, computational approaches are widely used to fill this gap. Many of these approaches are based on interpreting the brain structure as a graph with nodes as different brain regions and graph edges as the white matter fibers connecting the brain regions. Subsequently, various graph-based statistical as well as mathematical modeling techniques have focused on capturing the temporal, spectral, and spatial patterns of the brain in different brain states as well as diseases.

A key question in computational neuroscience is the elucidation of the exact relationship between the regional neural signals and their functional connectivity (FC) organization on one hand, and the underlying structural or anatomic connectivity (SC) on the other. This is what is referred to as the structure-function relationship in network neuroscience. Although these relationships may be addressed using purely statistical tools e.g., by evaluating the correlation between SC and FC (Hagmann et al., 2008), it is far more elegant to explore these relationships using biophysical mathematical models to serve as a formal interface between the physiology of the brain and its anatomic and connectivity properties (Honey et al., 2010). Moreover, mathematical models can be used for exploring biophysics underlying the whole brain imaging activity, and identify biophysically interpretable markers of disease (Honey and Sporns, 2008; Alstott et al., 2009; de Haan et al., 2012; Yang et al., 2016; Zimmermann et al., 2018; Singh et al., 2020). In this brief review we discuss different types of computational models at whole-brain and network level. We will explore their ability to recapitulate important aspects of empirical neural recordings, and end with new advances aimed at disease contexts and in the study of dynamic brain activity and behavior.

The contribution of this paper is threefold. First, it summarizes computational approaches that seek to understand how rich spatial and temporal structures can arise in brain activity using a set of underlying biophysical processes and their parameters. Second, it highlights several ways in which computational models are able to recapitulate empirical data across relevant scales ranging from electrophysiological scales corresponding to M/EEG signals to hemodynamic scales corresponding to fMRI signals. Our focus here is on whole brain network models of brain activity, rather than neuron or local circuit models. For readers interested in fine-scale or multi-scale modeling, we refer to a recent review (Glomb et al., 2021). Third, it discusses recent advances that attempt to make these complex computational models relevant to clinical and practical applications, often by deliberately simplifying these models and only keeping the most essential elements. The scope of the review is more focused on biophysical models rather than statistical graph theoretic metrics of SC-FC relationships. The field is greatly expanding and contains promising advances both from research and clinical standpoints. We believe that this article will serve to inform a wide readership about recent modeling work that are important from a practical point of view, yet are at risk of being inaccessible to a wide audience due to their mathematical or computational nature.

1.1. Magneto- and electroencephalography

We begin our report with a quick introduction to MEG and its older cousin EEG. These are non-invasive neuroimaging techniques that measure the postsynaptic potentials of pyramidal neurons in the brain (Biasiucci et al., 2019) by electrodes placed on the scalp (da Silva, 2013; Tivadar and Murray, 2019) or sensitive magnetic field detectors placed on or near the scalp (Hämäläinen et al., 1993; Gross et al., 2021). Both MEG and EEG possess high temporal resolution in the millisecond timescale, with spatial resolution determined by algorithms that solve the electromagnetic inverse problem (Biasiucci et al., 2019; Tivadar and Murray, 2019). The propagation of the signal from deep brain tissue to the scalp takes place due to the conductive properties of brain and head tissues, a phenomenon known as volume conduction which is also referred to as the electromagnetic forward problem (Kajikawa and Schroeder, 2011). To solve the inverse problem, source reconstruction algorithms (Michel et al., 2004) have been developed to estimate neural activity originating in deeper tissue, from the sensor signals and then superposing estimated functional activity on structural brain MRI, giving voxel-level signals. These algorithms have turned M/EEG into brain imaging technologies. However, it should be noted that source reconstruction is a mathematical ill-posed inverse problem; the reliability and spatial resolution of reconstructed signals are often functions of properties of the forward model as well as the signal and measurement noise. Once brain-wide neural sources and their activity have been reconstructed, further analysis might proceed analogously to imaging data like fMRI. The relationship between the functional activity from MEG/EEG/fMRI and the SC have been explored using various mathematical models, which we summarize in the next section.

2. Mathematical models

2.1. Coupled or networked non-linear neural mass models (NMMs)

At the mesoscopic scale, a lumped version of neural field was first described by Wilson and Cowan (1973); Cowan et al. (2016), consisting of an inhibitory and an excitatory population. Variants of this model were proposed by Da Silva et al. (1974) and Jansen and Rit (1995). Neural mass and neural field models are able to reproduce a range of dynamical behaviors that are observed in M/EEG, like oscillations in typical frequency bands (David and Friston, 2003), phase-amplitude-coupling (Onslow et al., 2014; Sotero, 2016), and evoked responses (Jansen et al., 1993; Jansen and Rit, 1995; David et al., 2005). Frequency-domain NMMs have also been proposed in order to model the M/EEG power spectrum (David and Friston, 2003; Bojak and Liley, 2005; Moran et al., 2007; Razi and Friston, 2016; Pereira et al., 2021).

By coupling together more than one lumped neural mass, one can start investigating the effect that inter-regional or remote connectivity and the resultant delays have on neural activity (Jirsa and Haken, 1996). Almost three decades ago, Jansen and Rit (1995) had already coupled together two neural mass models in order to simulate the effect of interactions between cortical columns on their activity. The Wilson-Cowan single oscillator model (Wilson and Cowan, 1972) too has evolved into a family of macroscopic coupled NMMs, with extensions to neocortical dynamics (Cowan et al., 2016; Deco et al., 2019b), controllability of brain networks (Muldoon et al., 2016), biomarkers in disease (Zimmermann et al., 2018), and second order statistics of observed brain signals (Deco et al., 2009; Abeysuriya et al., 2018; Wang et al., 2019; Byrne et al., 2020; Singh et al., 2020).

It is also common to use chaotic oscillators with little biophysical relevance (e.g., Yeung and Strogatz, 1999; Kuramoto, 2003), as the local neural mass oscillating around a given natural frequency (e.g., gamma band). Although there is no direct biological counterpart to these artificial models, they are sufficient to produce realistic synchrony and correlation properties (Deco et al., 2009; Cabral et al., 2014), e.g., the Kuramoto order parameter was linked to how dispersed firing is within a population (Byrne et al., 2017, 2020). Node-level Kuramoto oscillators were coupled using realistic brain connectivity and time delays determined by fiber lengths (Cabral et al., 2014; Finger et al., 2016).

Coupled NMMs are rapidly moving from mathematical curiosities to real applications in healthy and diseased brains (Honey and Sporns, 2008; Alstott et al., 2009; de Haan et al., 2012; Yang et al., 2016; Zimmermann et al., 2018; Singh et al., 2020) and can potentially be used to investigate biological processes such as neuromodulation (Shine et al., 2018, 2021; John et al., 2022). Such NMMs have also been simulated via the Virtual Brain platform—an open source neuroinformatics tool that provides access to brain network simulation tools (Ritter et al., 2013; Sanz-Leon et al., 2015). This tool has been used to study various diseases such as tumor (Aerts et al., 2018), Alzheimer's disease (Zimmermann et al., 2018), stroke (Falcon et al., 2015), and epilepsy (Jirsa et al., 2017).

Over the past few years, a vast literature has been developed on coupled neural models in a framework called dynamic causal modeling (DCM) (Kiebel et al., 2008; Pinotsis et al., 2012). DCM is a Bayesian framework for estimating directional interactions between coupled brain regions, i.e., the effective connectivity. A key concept underlying DCM is to treat the brain as a non-linear dynamic system that accepts multiple inputs and produces multiple outputs (i.e., MIMO model). This neuronal MIMO model is augmented with a regionally specific forward model that describes the mapping from neuronal activity to observations (MEG, EEG, local field potential; LFP, or fMRI). The neuronal and observation model together comprise a generative model that is subject to Bayesian inference. Stochastic DCMs are the most common generative models that have been used to test competing hypotheses about networks comprising only a few regions. Stochastic DCMs have usually proved to be highly successful at inferences of small networks (<10 nodes) than at large networks spanning the entire brain. However, increasing the number of regions or nodes in a stochastic DCM results in models with an enormous number of free parameters and profound conditional dependencies among the parameters, with exponential increases in computational time. Approaches such as spectral domain DCMs have been proposed to address this prohibitive estimation problem (Pereira et al., 2021), but there have been few attempts at whole-brain structural connectome integration with DCMs for M/EEG data.

2.2. Statistical models and graph metrics analysis

If the goal of structure-function modeling is merely to reproduce the second order covariances of neural activity, it is not necessary to employ a detailed biophysical signal generation model. Many recent studies have therefore focused on low-dimensional processes involving diffusion or random walks on the structural graph instead of mean-field models, under the assumption that such walks or network-diffusion processes can roughly simulate functional connectivity. These class of models were able to match or exceed the predictive power of complex coupled NMMs in predicting empirical FC (Abdelnour et al., 2014). Higher-order walks on graphs have also been quite successful; typically these methods involve a series expansion of the graph adjacency or Laplacian matrices (Meier et al., 2016; Becker et al., 2018). The diffusion and series expansion methods are themselves closely related (Robinson et al., 2016; Deslauriers-Gauthier et al., 2020; Tewarie et al., 2020), and naturally employ the so-called eigenmodes, or harmonics, of graph adjacency or Laplacian matrix. In these models, the eigenvectors of SC and FC are shared, and the functional eigenvalues are related by an exponentially-decaying function of the structural Laplacian eigenvalues. Hence these methods were generalized to yield spectral graph models whereby e.g., Laplacian harmonics were sufficient to reproduce empirical FC, using only a few eigenmodes (Atasoy et al., 2016; Abdelnour et al., 2018; Preti and Van De Ville, 2019). Although these models were initially demonstrated on fMRI-based FC, they have been readily demonstrated on M/EEG (Meier et al., 2016; Tewarie et al., 2019). Such spectral graph models are computationally attractive due to low-dimensionality and more interpretable analytical solutions. The SC's Laplacian eigenmodes may be thought of as the substrate on which functional patterns of the brain are established via a process of network transmission (Atasoy et al., 2016; Robinson et al., 2016; Abdelnour et al., 2018; Preti and Van De Ville, 2019; Glomb et al., 2020). Another class of models, the linear time-invariant system models based on the SC, have also been widely used to estimate brain controllability and energy required to transition to different states (Gu et al., 2015, 2017; Tang and Bassett, 2018; Stiso et al., 2019; Srivastava et al., 2022).

3. Challenges in inference of non-linear simulation models

The impressive richness and dynamic repertoire manifested by the above coupled neural mass and field models comes at a very high cost: that of model inference. These models are often mathematically complex (non-linearities), computationally complex (too many parameters) and statistically complex (ill-posed inference, non-unique solutions, broadly distributed posteriors). In addition, since they require long-time simulations involving numerical integration of highly coupled ODEs, they tend to be computationally costly, taking inordinately long time and computer power to create simulations that cover the entire parameter space.

The critical factor underlying these challenges is the coupling between a large number of non-linear NMMs. The complexity of coupled non-linear systems accommodates chaotic behavior which certainly imparts these models with a rich dynamic repertoire. Typically, coupled non-linear systems are characterized by exhibit discontinuous and abrupt shifts to oscillating, unstable, or chaotic behavior when network coupling or external driving force push the system over a bifurcation point. This produces a non-convex cost function with many local minima with discontinuous switching between regimes and model behavior—all adding up to highly challenging parameter inference. Recently, these challenges have been addressed using a two-step optimization procedure. First, all biophysical parameters of the local NMMs are set manually, without the use of empirical data, to be near an appropriate Hopf bifurcation point that gives the correct frequency (e.g., alpha band) (Sanz-Leon et al., 2015; Breakspear, 2017). In the second stage, optimization of the remaining few (global) parameters are performed by fitting to empirical FC metrics such as pairwise correlation or synchrony (Deco et al., 2009; Honey et al., 2009; Abeysuriya et al., 2018; Zimmermann et al., 2018; Demirtas et al., 2019; Wang et al., 2019). Even this optimization problem is usually no, zimmermann2018differentiationn-convex, hence grid search is employed instead of more sophisticated sampling methods. It is noteworthy that this two-step fitting procedure typically does not attempt to fit directly to actual time series of individual brain regions, only their second order FC, or in some cases the switching in time of granular “brain states” (Hansen et al., 2015; Deco et al., 2019a).

Two recent studies provide a thoughtful analysis of these inference problems. Hartoyo et al. (2019) highlighted the problem of unidentifiabiliy in whole brain models, where different parameter combinations can generate similar model predictions. Out of the 22 unknown parameters of a linearized network model they implemented, only one parameter was found to be identifiable when fitted to EEG data. Xie et al. (2019) show, using the Wilson-Cowan system as an archetype, that a large system of coupled NMMs may be pre-specified for firing rates pertaining to a specified frequency, but extending such a model to the network scale does not allow accurate, reliable, or efficient optimization. In particular, the manual/grid search fitting described above does not then also predict regional signals' power spectra. It was rightly noted in a recent review that coupled NMMs do not necessarily aim at maximizing the fit to the empirical brain signals but to reproduce specific temporal, spatial, or spectral features of the empirical data emerging at the macroscopic scale whose underlying mechanisms remain unclear (Glomb et al., 2021).

Due to these challenges, to date, coupled NMMs have not found widespread or routine usage in clinical or translational settings. One way to make the leap to practice is to carefully assess which modeling aspects may be considered useful, and which may be discarded. For instance, if the goal is to reproduce covariances rather than the signal spectra, then is a full-strength coupled NMM even necessary?

4. Recent advances in structure-function models

Below we highlight a few recent attempts that seek to enhance the practical applicability of complex computational models. Often these attempts have taken a step back from the complexity of current models and have instead explored deliberate simplifications rather than enhancements to the theoretical models. For example, many of the following coupled NMMs attempts involve fixing the local neural mass model parameters and then tuning only a small set of global parameters, simplify them to consist of fewer parameters, or linearize local masses to allow better inference.

4.1. Advances in coupled NMMs

In one study, the non-linear dynamics of coupled NMMs was transformed by applying a Lotka-Volterra type transform (Galadi et al., 2021). This enabled tracking of the temporal evolution of the global attractor and the non-stationary attractor landscape. Using this approach, the authors were able to compare awake state with different stages of sleep with BOLD fMRI data (Galadi et al., 2021). A recent study modeled noisy Stuart-Landau oscillators to simulate BOLD fMRI and found that low-level states of consciousness were associated with structural hub regions losing stability and an increase in homogeneity in the brain regional dynamics (López-González et al., 2021). Another study on the levels of consciousness demonstrated that in addition to the descriptive distance metrics classifying the different levels of consciousness, the extent of perturbation required to switch between the levels provide complementary information about the levels. Here, the extent of perturbation was quantified by introducing a perturbation term in the coupled NMMs (Sanz Perl et al., 2021). A DCM approach was used to investigate modulation of gain by the cognitive control systems in early psychosis (Burgher et al., 2021).

In addition to exploring brain states, a recent study demonstrated how NMMs can be enriched by incorporating brain-wide transcriptional data that captures the regional inhibitory and excitatory receptor gene expression. This approach ensures model parsimony while incorporating regional heterogeneity to simulate rich BOLD fMRI dynamics (Deco et al., 2021). While these studies explored biological and clinical applications of the NMMs, a recent study also showed how NMMs can be leveraged to infer virtual FC or SC. The inferred virtual connectomes's classification of healthy vs. AD as well as different age classes was comparable to the classification done by the empirical connectomes (Arbabyazd et al., 2021). Lastly, a next generation NMM was employed to investigate the role of structural connectivity in regulating the emergent functional dynamics (Gerster et al., 2021).

Recent advances were also made using the Virtual Brain platform. One of the studies demonstrated selectively stimulating different fiber tracts in a personalized virtual brain model to simulate network effects in deep brain stimulation for patients with treatment-resistant depression (An et al., 2022). Another application of the virtual brain platform was to demonstrate similarity between neuronal cascades of firing rate fluctuations and resting state networks dynamics (Rabuffo et al., 2021).

Coupled NMMs have also been leveraged to fit to M/EEG. Recently, coupled phase oscillator models was used to reconstruct coupling between frequency bands (interlayer connectivity) obtained from MEG (Tewarie et al., 2021). Another study focused on estimating coupled NMMs from EEG recordings of epilepsy patients to develop a closed-loop deep brain stimulator (Chang et al., 2021).

4.2. Advances in statistical models and graph metrics analysis

To improve the structure-function model fits, a new proposal relates the eigenvalues of the structural connectivity and functional networks using the Gamma function instead of the exponential decay, producing a reliable prediction of functional connectivity with a single model parameter (Cummings et al., 2022). They also found that adding small levels of long range interhemispheric connectivity greatly improves model performance, since it is well-known that DTI-based fiber connectivity methods significantly under-estimate inter-hemispheric connections. Another recent study extended the notion of graph Laplacian by introducing global connectome coupling and a conductance speed parameter, which gave rise to a “tunable” Laplacian (Xie et al., 2021) (Figure 1A). Since delays becomes phases in Fourier space, this study proposed a so-called “complex Laplacian”. They showed that with the right selection of the two model parameters, it is possible to “steer” the (complex) eigenmodes of this Laplacian in such a manner that a small number of them can reproduce any given canonical functional network (like default mode, salience, etc.) at relevant frequencies. Examples of these complex eigenmodes were shown in Figure 1B, copied from the original publication (Xie et al., 2021). This is exciting, since it can accommodate a large repertoire of microstates and their concomitant spatial patterns—an essential characteristic of real brain activity. This rich repertoire could be engaged using only two steerable model parameters, suggesting the possibility that complex behavior may be achievable by simple and parsimonious mechanisms. A latest study compared the performance of 40 structure-function coupling models in capturing the variance in FC (Zamani Esfahlani et al., 2022). They found that having regional heterogeneous coupling is more accurate than having a global coupling between structure and function. Moreover, specific models were more accurate in capturing the coupling between specific brain regions. Lastly, this coupling decreased with age. Another study showed a stronger association of structural eigenmodes (in comparison to structural connections) with phase and amplitude connectivity of MEG at short time-scales (Tewarie et al., 2022).

Figure 1

(A) Overview of the complex Laplacian eigenmode decomposition. The complex Laplacian matrix is computed from the structural connectivity matrix (C) and distance adjacency matrix (D) which were extracted from diffusion MRI. Then, eigendecomposition of the Laplacian provided the complex structural eigenmodes. The spatial similarities between these eigenmodes and the canonical functional networks in fMRI was done for model parameter tuning. (B) Complex Laplacian eigenmodes for different coupling strength and wave number. Top row: eigenmodes from the real Laplacian matrix with coupling strength = 1. Bottom three rows: Complex Laplacian eigenmodes for small wave number (top), high wave number (middle), and high wave number and coupling strength (bottom). Figure is extracted from Xie et al. (2021).

Network control models have also been recently extended and applied to investigate controllability in drug-resistant focal epilepsy in children (Chari et al., 2022), functional dynamics mediated by controllability of the structural connectivity (Gu et al., 2022), prediction of aphasia recovery using controllability measures (Wilmskoetter et al., 2022), association of controllability measures with normative negative affect variability (McGowan et al., 2021), prediction of positive psychosis spectrum symptoms with controllability measures (Parkes et al., 2021), role of controllability measures in language tasks (Medaglia et al., 2021), altered energy landscape and stability of working memory representations in schizophrenia (Braun et al., 2021), structural control energy in psychosis vulnerability (Zoller et al., 2021), and dynamics of controllability measures during seizure progression (Scheid et al., 2021).

4.3. Spectral graph model for neural oscillations

Above-noted spectral graph models were successful in capturing steady-state, stationary second order covariances in brain activity, but did not possess oscillatory behavior. Hence their extension to accommodate a larger repertoire of dynamically-varying microstates or rich power spectra at higher frequencies typically observed on EEG or MEG would require a full accounting of axonal propagation delays as well as local neural population dynamics within graph models, as previously advocated (Cabral et al., 2011). Band-specific MEG resting-state networks were successfully modeled with a combination of delayed neural mass models and eigenmodes of the structural network (Tewarie et al., 2019), suggesting delayed interactions in a brain's network give rise to functional patterns constrained by structural eigenmodes.

Recently another effort was undertaken to characterize wide-band brain activity using graph harmonics in closed form (i.e., requiring no time-domain simulations), a rarity in the field of computational neuroscience (Raj et al., 2020). This spectral graph model (SGM) of wideband brain activity produced realistic power spectra that could successfully predict both the spatial as well as temporal properties of MEG empirical recordings (Raj et al., 2020). Remarkably, the model has very few (seven) parameters, all of which are global and not dependent on local oscillations. This method therefore exemplifies the power of graph methods in reproducing more complex and rich repertoire of brain activity, while keeping to a parsimonious approach that skirts the need for high-dimensional and non-linear oscillatory NMMs. SGM can be interpreted as a linearized coupled NMM.

The SGM consists of two levels. At the mesoscopic (local or node) level the excitatory and inhibitory neuronal populations are characterized by time constants τe, τ and gains gee, gei, and gii. At the macroscopic level long-range pyramidal excitatory neuronal populations are characterized by time constant τG, and a coupling constant α. A pictorial depiction of this model is shown in Figure 2A. Remarkably, the SGM is available in closed-form in Fourier space as a summation over the eigenmodes of the complex Laplacian :

Figure 2

(A) Flow chart showing the structure of SGM and inference of SGM parameters by comparing modeled frequency spectra with the empirical frequency spectra captured by MEG. Input p(t) is white Gaussian noise and is the same for both excitatory as well as inhibitory signals. (B) Spectral and spatial correlations obtained after optimization for all the subjects with three different cost functions: (1) spectral + spatial correlation, 2) spatial correlation, 3) spectral correlation. (C) Comparison of the empirical and modeled spatial distribution of the alpha frequency band, averaged over all the subjects. (D) Comparison of empirical and modeled frequency spectra, averaged over all brain regions and subjects.

where ω is the angular frequency, u(ω) are the eigenmodes and λ(ω) are the eigenvalues. The term Hlocal(ω) refer to the mesoscopic model's transfer function, P(ω) to the input noise spectrum, and FG(ω) is the neural impulse response transfer function (Gamma shaped). Details and derivation of this equation from first principles are contained in Raj et al. (2020) and Verma et al. (2022a). The structure of SGM is similar to that of Jansen-Rit model (Jansen and Rit, 1995; Sanz-Leon et al., 2015). Further methodological comparisons can be found elsewhere (Verma et al., 2022b).

Since the development of the SGM, there have been various extensions as well. We recently improved SGM by enhancing the biophysical interpretability of the mesoscopic model (Verma et al., 2022a), and subsequently demonstrated various stability properties of this improved model (Verma et al., 2022b). In the latter, we also outlined a strategy to capture temporal fluctuations in MEG by fitting the SGM parameters to MEG frequency spectra at various time points. The mesoscopic SGM model parameters were also estimated to differentiate Alzheimer's disease with healthy controls (Ranasinghe et al., 2022), as well as different stages of sleep (Fan et al., 2022).

4.4. Parameter inference of SGM from empirical MEG Data

The key benefit of linear graph models like the SGM over coupled NMMs is that parameter inference is far simpler and easily accomplished in a data-driven fashion without the need for manual selection. While prior works showed that simple simulated annealing-type optimization can accomplish this task, there is still an open question regarding the exact form of the cost function and which data features to use for inference. Here we highlight very recent developments in this area. Broadly, prior NMM methods have focused on FC as the target of optimization, whereas prior SGM studies used wideband MEG spectra directly. Here we show that SGM is capable of fitting simultaneously to both spectra and spatial distributions of the alpha frequency band.

We used the same MEG dataset reported previously (Xie et al., 2020; Verma et al., 2022a), comprising of 36 healthy adult subjects (23 males, 13 females; 26 left-handed, 10 right-handed; mean age 21.75 years, age range 7–51 years), having 1 min MEG recordings at rest and eyes closed. SGM was computed on the structural connectomes of these subjects, and both modeled and empirical MEG were converted into power spectral densities (PSD) in dB scale. Pearson's r between modeled and MEG PSD was computed for all 68 brain regions. Its average r across all regions is referred to as the spectral correlation. Next we calculated the spatial correlation by obtaining the regional distribution of alpha band (8-12 Hz) raw power of both model x and MEG y. Then, the spatial correlation was defined as xǁ(C+wI)ǁy, where C is the row degree normalized structural connectivity matrix, I is the identity matrix, w is an empirical weight, and ǁ(C+wI)ǁ is the row normalized version of C+wI. The objective function for optimization and estimation of model parameters was the sum of spectral and spatial correlations. The details of the optimization procedure are in the Supplementary material.

4.4.1. Model parameter estimation results

We firstly tested different values of w in the spatial correlation weight of C+w*I. The values being tested were 0.01, 0.1, 0, 1, and 10. For each of the values of w, no statistically different spectral correlations were observed based on a Tuckey's HSD test on Fisher's r-to-z transformed correlations. However, the spatial correlations were highest for w = 10. Therefore, we finally chose a value of w = 10 out of 0.01, 0.1, 0, 1, and 10 for subsequent analysis. After performing model parameter optimization for all the subjects with w fixed at 10, we see that the spectral as well as spatial correlations are high for all of them, as shown in the left-most subplot in Figure 2B. In addition, the spectral correlation reduces significantly when spatial correlation alone is used as a cost function for optimization (middle subplot in Figure 2B), and the spatial correlation reduces significantly when spectral correlation alone is used as a cost function (right-most subplot in Figure 2B). Therefore, we used an unweighted sum of spectral and spatial correlations for final analysis. Visually, we see an accurate match for the spatial distribution of the alpha band (Figure 2C) as well as the spectral shape (Figure 2D).

This represents a noteworthy advance over prior publications, since the inferred model is informed equally by spectral as well as spatial features of empirical MEG—unlike most comparable methods described in this paper. We believe future computational approaches may similarly benefit from the use of multiple feature classes of empirical functional data.

In addition to the approach above, NMMs have also been linearized previously. Such linearized NMMs, both at a regional level and for the entire brain network (Moran et al., 2007; Deco et al., 2014; Abeysuriya and Robinson, 2016; Gabay et al., 2018; Raj et al., 2020; Verma et al., 2022a) provide for tractable parameter inference in comparison to their non-linear counterparts. However, this simplification limits the repertoire of dynamical solutions that these models can exhibit, similar to SGM.

5. Conclusions

We briefly surveyed the burgeoning field of computational neural mass models on the whole brain connectivity network. Such models have been widely used to infer the underlying biophysics of the whole brain activity from non-invasive neuroimaging data. They have been used to study various neurological disease as well as brain states.

We noted that most current methods were designed to achieve resemblance to functional connectivity—a second-order feature based on measuring co-activity between brain regions—rather than to the data itself. Although this is justified in specific use-cases, we argued that it would be necessary to develop more complete models that can achieve both goals. We noted potential scientific and clinical impact of current approaches, but enumerated several challenges associated with model inference. In summary, the challenges are: (i) parameter identifiability issues when model parameters are heterogeneous across brain regions, (ii) parameter inference difficulty due to abrupt shift in model behavior with small shifts in parameters, and (iii) fitting directly to actual time series instead of second order FC.

We identified many recent and exciting advances, in both non-linear coupled NMM approaches, as well as the recent emergence of linearized spectral graph models. The latter in particular present an attractive trade-off between model complexity, accuracy and practicality. We conclude that while the field has seen various advancements with applications, further advances are needed to overcome identified issues. Therefore in this review, we suggest the following directions for further advancement of this field:

Develop computational models that can capture temporal, spatial, as well as spectral features of whole brain imaging.

Develop parsimonious biophysical computational models that can be used in clinical or translational settings.

We are excited by the recent highlighted advances, that have together taken meaningful steps in this direction.

Author contributions

AR, PV, and SN conceptualized and designed the SGM results. AR and SN supervised the SGM results. PV generated and compiled the SGM results. AR created the first draft. All authors contributed to the writing and editing of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was partially supported by NIH grants R01AG072753, R01NS092802, R01NS018341, R01EB022717, R01DC013979, R01NS100440, R01DC176960, R01DC017091, R01AG062196, UCOP-MRP-17-454755, DOD CDMRP Grant W81XWH1810741, Alzheimer's Association Grant AARFD-22-923931, and an industry research contract from Ricoh MEG Inc.

Conflict of interest

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.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.