Complex visual cognitive function based on a large-scale neurovascular and metabolic coupling mechanisms model in whole brain

1.   Introduction

Neurovascular and metabolic coupling (NVMC) refers to the intricate interplay among neurons, glial cells, and distal cerebral blood vessels, which collectively maintain the homeostasis of brain tissue function. Typically, an increase in cerebral blood flow (CBF) ensures adequate energy supply to active neurons and facilitates the clearance of metabolic byproducts during neuronal activity ^[1,2]. At its core, NVMC translates neural signaling into metabolic and energetic demands, subsequently generating chemical and mechanical signals to regulate blood flow. NVMC serves dual critical roles. First, it sustains neuronal activity by dynamically supplying energy substrates, thereby enabling prolonged neural function. Second, it forms the physiological basis of blood oxygen level-dependent (BOLD) signals in functional magnetic resonance imaging (fMRI) ^[3,4]. The cognitive process evoked by visual stimulations is a complex process in which the brain integrates visual information with related cognitive functions. In particular, under different categories of visual stimulus (e.g., famous faces, unfamous faces, and scrambled faces), neurons in different regions of the brain perform the complex tasks of perception, memory, cognition, and decision-making through precise temporal and spatial synergies. This process not only relies on neural activity in the brain but also involves NVMC mechanisms.

According to the NVMC mechanism, the relationship between blood flow and metabolic hemodynamic indices is highly temporally and spatially consistent. Under visual stimulation, blood flow in the visual cortex and related brain regions increases, and cerebral metabolic rate of oxygen (CMRO₂) and glucose metabolism are enhanced accordingly. Different visual stimulations evoke different degrees of neural activity, which in turn leads to different metabolic demands. For example, studies have shown that the recognition of famous faces typically involves brain regions associated with the memory system, such as the temporal lobe and hippocampus whereas unfamous faces rely more on processing in the visual cortex ^[5]. Scrambled faces or meaningless visual stimulation, in contrast, may provoke additional cognitive load, requiring more attention and perceptual resources ^[6]. Researchers also found that processing familiar faces (e.g., famous faces) involves higher levels of cognitive processing, such as memory extraction and affective responses, than processing unfamiliar faces or complex visual stimulation ^[7].

In order to fully understand the mechanisms of NVMC during visual cognition, it is difficult to reveal its complex phenomena with a single data type, and integrating multimodal data is essential. Advances in magnetic resonance imaging (MRI) have offered noninvasive methods to characterize large-scale dynamical study in the human brain at the structural and functional levels. fMRI measurements of BOLD signals have been the predominant method for extensively investigating the functional organization of human brain activity. Electroencephalogram (EEG) detects changes in neural electrical activity in the cerebral cortex to obtain effective connectivity between brain areas. Multi-modal techniques provide more comprehensive and rich information to reveal the hidden relationship between brain structure and function, while, the combination of multi-modal information can better understand the spatiotemporal properties of the brain. More importantly, data from different modalities can better expose the differences between individuals, which is especially important for the study of brain cognition, brain injury, and neurological diseases. Multi-modal techniques in depression ^[8] and schizophrenia ^[9,10,11] have provided a holistic characterization of these diseases ^[12].

Moreover, modeling of the brain also contributes to a comprehensive understanding of NVMC mechanisms and thus reveals the dynamic interactions between neural and vascular functions in visual tasks. The realization of brain functions, including cognition, memory, and perception, relies on the interactions among neurons spanning multiple temporal and spatial scales. Thus, models of the brain encompass microscale (individual neurons), mesoscale (groups of neurons), and macroscale (large-scale brain networks) ^[13]. Since NVMC not only reflects the local metabolic hemodynamic demands but also reflects the integration and distribution process of the brain through the functional connections in the whole brain network. The choice of model scale is a non-neglectable issue, which determines the scope of the field of view that we can analyze. Whole brain network models are widely utilized in current research, primarily leveraging multimodal data to construct networks and extract topological organizational features at both local and global scales. These models are widely applied, for instance, in studying network alterations and classification in disorders such as Alzheimer's disease (AD), attention-deficit/hyperactivity disorder (ADHD) and so on ^{[14,15,16,17]}. However, such frameworks lacked physiological model and remain limited in their capacity to support deeper theoretical exploration of underlying mechanisms. In summary, current whole-brain models exhibit limited biological interpretability, ignore heterogeneous physiological structures, and some models lack the ability to incorporate realistic multimodal data inputs. Furthermore, there are few studies on task states in the model.

To address these limitations, we integrated multimodal data into a large-scale brain network biophysical model, which can accurately analyze the dynamics of neural activity, metabolism, and hemodynamics, simulate brain function, and provide a novel understanding of how the brain processes visual information. This research not only explains existing observations but may also predict unseen phenomena. In this study, we first used a previously developed brain network model. Individualized physiological parameters were obtained by an optimized algorithm to simulate real NVMC activities. Second, we investigated the differences in the temporal and spatial dimensions of three different visual stimulations, as well as the different performances in terms of node and network structure. Finally, we combined the individualized physiological parameters with cognitive sub-networks, aiming to explore the heterogeneity of cognitive functions in different subnetworks.

2.   Materials and methods

2.1.   Ethics approval of research

The study of this publicly available dataset was approved by the Cambridge University Psychological Ethics Committee. Written informed consent was obtained from each participant prior to and following each phase of the experiment. Participants also gave separate written consent for their anonymized data to be freely available on the internet.

2.2.   Clinical trial registration

We used visual stimulation task state data from 16 healthy volunteers for model application ^[18]. These data included EEG, MRI, and fMRI, and this data was used to develop personalized multimodal brain network models, while the multimodal data was used to improve spatial and temporal resolution over any one modality. Visual stimulation data can be downloaded from https://openfmri.org/. The preprocessing of the data is described in detail in the supplementary material.

In visually stimulated EEG data: EEG were measured in a photomagnetically shielded room using a MediTech Neuromag Vectorview 306 system (Helsinki, FI). EEG data were recorded synchronously using a 70-lead Easycap EEG cap with an electrode arrangement conforming to the extended 10-10% system and a data acquisition frequency of 1100 Hz.

The MRI data were collected from a Siemens 3T TIM TRIO in the visually stimulated data (Siemens, Erlangen, Germany). A standard 1mm isotropic T1-weighted "structural" image was acquired using an MPRAGE sequence (TR 2,250 ms, TE 2.98 ms, TI 900 ms, 190 Hz/pixel; flip angle 9°) ^[18].

The fMRI data were acquired using an EPI sequence of 33, 3 mm-thick axial slices in the visually stimulated data (TR 2000 ms, TE 30 ms, flip angle 78°). Slices were acquired in an interleaved fashion, with odd then even numbered slices (where slice 1 was the most inferior slice) and a 25% distance spacing (increased where necessary to cover the whole of cortex), resulting in a range of voxel sizes of 3 × 3 × 3.75 mm³ to 3 × 3 × 4.05 mm³ across participants. A total of 210 volumes were acquired in each of 9 runs (note that the 3 initial TRs were discarded to allow saturation of T1 effects) ^[18].

2.3.   Dynamic modeling

To elucidate the mechanistic differences underlying brain activities under distinct visual stimuli, we established a whole-brain-scale dynamic model that translates neural mass firing into metabolic-hemodynamic responses. The model was developed using multimodal neuroimaging data, including EEG, MRI, and fMRI. The methodological framework employed in this study is illustrated in Figure 1.

Figure 1.  The methodological framework of whole-brain-scale computational model. A: The model inputs consist of multimodal data (e.g., EEG time series and MRI) as well as the source reconstruction process (source estimation of dipoles and central source calculation based on brain atlas). B: Whole-brain-scale computational framework and detailed node-level models. these node-level models include: Reconstruction of LFP, Reverse neural mass model, Metabolic Hemodynamic model and Simulated fMRI results. C: Empirical fMRI results. The simulated results are compared with empirically measured fMRI data for validation. An error cost function is subsequently applied to optimize the parameters of the node-level models.

DownLoad: Full-Size Img PowerPoint

First, based on the preprocessed and personalized MRI data, the averaged ERP data were used for source reconstruction, i.e., converting to dipoles in cortical space using the brainstorm3 toolbox (https://neuroimage.usc.edu/brainstorm/Introduction), and the method used for this process was the sLORETA inverse transformation. Second, the central dipoles of each brain region were obtained based on the Desikan-Killiany (DK) ^[19] partition in order to construct a simplified distributed dipole model. Last, we reconstruct local field potentials (LFP) based on dipole theory. Furthermore, the LFP was considered as the neural dynamic energy used to drive the model. See supplementary material for supplementary material.

The reverse neural mass model (RNMM) was improved by the neural mass model (NMM) ^[20], which decomposes the LFP into excitatory and inhibitory neural activities (excitatory postsynaptic potential (EPSP) and inhibitory postsynaptic potential (IPSP)) by optimizing key parameters of the model. Key parameters of the model included excitatory synaptic gain (A), inhibitory synaptic gain (B), and number of intercellular ground connections (C). Parameters A and B modulated the sensitivity of excitatory and inhibitory synapses and determined the excitatory and inhibitory properties of the model. Parameter C described the interaction between pyramidal cells and interneurons. The differential equations describing the above theory are given in the supplementary material.

The RNMM optimal process is described as follows: $LF{P}_{input}$ was normalized to ensure that its magnitude was similar to $LF{P}_{output}$. To ensure the accuracy of the PSP, the values of A, B, and C were adjusted to maximize the match between the $LF{P}_{input}$ and $LF{P}_{output}$. The optimization method was the simulated annealing method, and the error cost function was root mean square error (RMSE):

where n is the number of samples in the LFP.

The simulated neural activities (EPSP and IPSP) were fed to the metabolic hemodynamic model (MHM) ^[21] to simulate the hemodynamic indexes for each region of interest (ROI). Briefly, changes in PSP are linked to the corresponding changes in glucose consumption. Then, the total glucose consumption was calculated as a weighted average of the excitatory and inhibitory contributions. The glucose variables were then directly related to the metabolic rates of oxygen for excitatory and inhibitory activities, as well as to the total oxygen consumption. The equations relating these biological processes are given in the supplementary material.

As in the study of Sotero and Trujillo-Barreto ^[21], there were specific details for the derivation of the equations describing this relationship. The MHM modifies the previous model ^[22], with the additional assumption that CBF is not coupled to inhibitory activity.

Furthermore, the Balloon model ^[23] was employed for linking the output of the metabolic and vascular models to normalized cerebral blood volume (v) and deoxyhemoglobin content (q). Knowing q and v, the BOLD signal was calculated as in ^[23].

Given q and v, the BOLD signal was described as:

The empirical parameters of the above differential equation ground are described in the supplementary material.

2.4.   CBF calculation and functional connectivity analysis

In this study, the CBF was averaged from the model calculations to obtain the CBF values averaged for the task. Furthermore, we computed the Pearson correlation coefficient (N = 68) between the BOLD time series of the DK brain region pairs and thus obtained the functional connectivity matrix at the level of the whole grey matter for each subject. Since whole-brain negative correlations are controversial, we restricted the analysis to positive correlation values above a threshold of 0.2 to eliminate weak correlations with noise-induced correlations. This parametric restriction of the analysis was taken from ^[24], as lower thresholds may include false-positive connections, while higher thresholds may lose useful connections. Binarization was performed using the threshold limit, which was set to 1 if the functional connectivity was less than the threshold limit, and 0 otherwise. For brain region x, the functional connectivity strength (FCS) was calculated as the average of the functional connectivity between x and all other brain regions (including 0 entries).

2.5.   Brain network analysis

In this section, we calculated the global and local indicators separately. The indicators included average degree (AD), clustering coefficient (CC), characteristic path length (CPL), and betweenness centrality (BC). The AD describes the degree of edge denseness; a higher average degree indicates denser edges between nodes in the network. The CC is a measure of functional segregation, which is the ability for specialized processing to occur within densely interconnected groups of brain regions. The CPL is a measure of functional integration, which is the ability to rapidly combine pieces of specialized information from distributed brain regions. The BC quantifies the importance of a node within a graph by measuring its involvement in shortest paths ^[25]. The specific equations for the calculations are described in the supplementary material.

3.   Results

3.1.   Task-state CBF changes in various brain regions

Since the neural activity as the input to the model was normalized, the CBF calculated by the model was a variable quantity. In the resting state, the CBF calculated by the model was assumed to be constant, i.e., 1. The spatial distribution of CBF in the resting state is shown in Figure 2.

Figure 2.  Spatial distribution maps of CBF at rest.

DownLoad: Full-Size Img PowerPoint

The brain regions with higher spatial distribution of CBF were more similar under the three visual stimulations (Figure 3). The cortex with higher CBF was mostly located in the right inferior temporal, the left transverse temporal, and the right lateral occipital, with ratios to resting state as follows: famous faces (1.5707, 1.4991, 1.4265), scrambled faces (1.5408, 1.4692, 1.4248), and unfamiliar faces (1.5586, 1.4859, 1.4328).

Figure 3.  Spatial distribution maps of CBF for three visual stimulations: Famous faces (A), scrambled faces (B), and unfamiliar faces (C).

DownLoad: Full-Size Img PowerPoint

3.2.   The relationship between CBF and FCS in different visual stimulations

3.2.1.   Differences in CBF and FCS among individuals

In this study, spearman correlation was used to investigated individual difference variables (age, gender, stimuli number) with CBF and FCS. We found that under famous face stimulation, age was significantly correlated with FCS (r = 0.3016, P < 0.01), and there was no significance between any of the other individual variables and CBF and FCS. See Table 2 in the supplementary material for details of individual variables and CBF and FCS.

3.2.2.   Coupling analysis of CBF-FCS among groups

The changes in CBF-FCS coupling at the level of whole-brain gray matter under three visual stimulations are shown in Figure 4. Significant positive correlations in CBF-FCS coupling were found at the group level in subjects with visual stimulation of famous faces (r = 0.5083, P < 0.001), scrambled faces (r = 0.5957, P < 0.001), and unfamiliar faces (r = 0.5677, P < 0.001). There was little difference in correlation between the three groups from the coupling results.

Figure 4.  (A) to (C) are CBF-FCS coupling changes at the whole gray matter level. (A) for famous faces, (B) for scrambled faces, and (C) for unfamiliar faces. (D) The CBF-FCS correlation coefficients across brain regions for different visual stimulus groups.

DownLoad: Full-Size Img PowerPoint

3.2.3.   Analysis of the ratio of CBF/FCS among groups

The changes in CBF/FCS at the level of region under three visual stimulations are shown in Figure 5. The spatial distributions of the three metrics were more similar among subjects under different stimulations.

Figure 5.  CBF/FCS map, normalized to Z-score at level of whole gray matter. (A–C) results of CBF/FCS for three visual stimulations: Famous faces, scrambled faces, and unfamiliar faces, respectively.

DownLoad: Full-Size Img PowerPoint

Under the visual stimulation of famous faces, the larger regions of CBF/FCS were the left transverse temporal (0.6260) and the right inferior temporal gyrus (1); for scrambled faces, the larger regions of CBF/FCS were the left transverse temporal gyrus (1) and the right inferior temporal gyrus (0.9409); in terms of unfamiliar faces, the larger regions of CBF/FCS were the left transverse temporal gyrus (1), the right inferior temporal gyrus (0.86830), and lateral occipital (0.7497).

3.3.   Network topology analysis of whole-brain nodes in different visual stimulations

We also analyzed the local indicators and computed the betweenness centrality under the three task states for evaluating the importance of nodes in each brain region, and the results were normalized and are shown in Figure 6(A)–(C). It was found that the median centrality differed under the three visual stimulations: Under famous faces, the left precentral gyrus, the right supramarginal gyrus, and the right lateral occipital lobe had higher values; for scrambled faces, the right inferior parietal lobe, and the right lateral occipital lobe had higher values; and for unfamous faces, the left pars orbitalis, and the right internal olfactory region had higher values.

Figure 6.  Local and global indicators of the functional connectivity matrix. (A) Betweenness centrality under famous faces. (B) Betweenness centrality under scrambled faces. (C) Betweenness centrality under unfamiliar faces. (D) Global metrics under three different visual simulations, which are AD, CPL, and CC.

DownLoad: Full-Size Img PowerPoint

Moreover, three types of global metrics, average degree, feature path length, and global clustering coefficient, were calculated, as shown in Figure 6(D). It was found that the average degree distribution and feature path length were all distributed between the range of 0.9 to 1. The global clustering coefficient removed the anomalies, and the data were distributed between 0.92 and 1.01 when the subjects were subjected to the three different visual stimulations. After that, a t-test was performed between the three task groups, and P > 0.05 was obtained for all three global metrics, indicating that the three global metrics were not significantly different between the three task states.

3.4.   Optimized parameters for six functional sub-networks

To further investigate the link between personalized physiological parameters and cognitive subnetworks, the whole-brain ROIs were regrouped into 12 subnetworks, including the default mode network (DMN), attention network (ATT), visual recognition network (VRN), auditory network (ADN), sensory-motor areas (SMA), and subcortical network (SN) in both left and right hemispheres, according to the classification method of Tao et al (Figure 7(A)) ^[26].

Figure 7.  Classification of 12 sub-networks obtained in DK mapping and optimized results for three parameters in different sub-networks. (A) Subnetwork segmentation results. (B-D) Results of optimization of parameters A, B, and C in subnetwork scales. (*P < 0.05 vs. ATT. & P < 0.05 vs. ADN. #P < 0.05 vs. DMN. $P < 0.05 vs. SN.!P < 0.05 vs. VRN).

DownLoad: Full-Size Img PowerPoint

The RNMM uses a simulated annealing algorithm to optimize parameters A, B, and C. The results show significant differences in the different sub-networks, as shown in Figure 7(B)–(D). Parameter A represents excitatory synaptic gain, and the results showed that the largest values were in the sub-networks ATT and DMN, with values ranging from 4.25 to 5.25, followed by ADN, VRN, SN, and SMA, with values ranging from 1.25 to 2.75. The two-by-two significance between the sub-networks was also calculated, and it was found that there was a significant difference between the six sub-networks, but there was no significant difference between the task states.

Parameter B represents inhibitory synaptic gain, and the results showed that the largest values were found in the sub-networks ATT and DMN, with values ranging from 6.5 to 7.75; followed by ADN, VRN, SN, and SMA, with ranges from 2.0 to 4.0. The two-by-two significance between the sub-networks was also calculated, and it was found that there was a significant difference between the six sub-networks, but there was no significant difference between the task states.

Parameter C represents the average number of connections between pyramidal neurons and intermediate neurons, and the results showed that the largest values were found in the sub-networks ATT and DMN, with values ranging from 1.1 to 1.7, followed by ADN, VRN, SN, and SMA, with values ranging from 0.3 to 0.9. The two-by-two significance between the sub-networks was also calculated, and it was found that there was a significant difference between the six sub-networks, but there was no significant difference between the task states.

4.   Discussion

In this study, we integrated multimodal data into a macroscopic model of the brain that accurately analyzes the dynamics of neural activity and hemodynamics, simulates the brain's visual cognitive functions, and provides a new understanding of how the brain processes visual information. We computed and analyzed differences in CBF, parameters of NVMC (CBF-FCS, CBF/FCS), network topology parameters (AD, CPL, GCC), and individualized physiological parameters in response to three different visual stimulations (famous faces, unfamiliar faces, and confusing faces), The results are discussed in detail below.

4.1.   The network model predicts fMRI time series in personalized subjects

The brain network model used in this study is a dynamic model in which a single neural mass model used in individual brain regions simulates the overall neural activity, and a functional connectivity matrix obtained using LFP is used to connect the node models. First, in contrast to EEG, LFP recordings are less susceptible to artifacts because the electrodes are placed inside the brain. However, LFP recordings are mostly limited to animal studies ^[27]. LFPs primarily measure the synaptic potentials of neuronal populations in the vicinity of the recording electrode and are closely related to transmembrane currents ^[28]. Crucially, although the LFP signal correlates with EEG signals measured at the scalp, as both measure voltage from pooled synaptic potentials, the LFP signal can reveal a very different pattern of activity from the EEG signal ^[29]. We attempt to convert scalp EEG more accurately to the LFP activity in the brain source domain: Firstly, scalp EEG inverse is transformed to the brain source activity (dipole), while we use the dipole theory to calculate the average LFP of the dipole at the center of each brain region, which reduces errors in anatomical structures such as cerebrospinal fluid and skull. Second, in contrast to previous neural field models where the inputs are used to simulate neural activity with noise ^{[20,30,31,32,33]}, we use a network model of neural mass modeled by LFP driven by a functional connectivity matrix constituting the weights of the short- and long-duration connections between different nodes.

Key parameters of the model include A, B, and C, which were optimized. Given the variability in parameter values across literature ^{[20,31,32,34]}, we select intervals covering most reported ranges: ^[1,10] for parameter A, ^[5,20] for parameter B, and ^[1,2] for parameter C. Initial values are randomly sampled within these ranges. Key parameters are optimized via a simulated annealing algorithm, with convergence triggered when the RMSE between the calculated LFP and optimized LFP falls below 0.001. Parameter updates follow a normal distribution with mean 0 and standard deviation X (uniformly sampled from 5 to 10 per iteration), dynamically balancing global exploration and local refinement. For validation, the optimized algorithm is tested on resting-state data. Specifically, the RNMM exhibits parameter optimization errors below 0.001 across diverse datasets, demonstrating strong data adaptation. See Figure 1 in the supplementary material for details of optimization of resting-state LFP data. Moreover, artificial rectangular wave stimulations from the MHM reveal that stimulus duration critically modulates CBF, BOLD responses: long activations induce higher-amplitude and sustained CBF, BOLD activation (peak amplitude difference: 0.15, 0.23 respectively), whereas shorter activation result in rapidly decaying weak responses. See Figure 2 in the supplementary material for details of hemodynamic indices induced by neural activities of different durations. These findings validate the model's sensitivity to temporal dynamics of neural activity and its capability to capture stimulus-pattern-dependent hemodynamic regulation.

In addition, under famous face stimulation, age is significantly correlated with FCS, indicating that even in a normal population, the strength of brain functional connectivity increases from age 21 to 31 under famous face stimulation. This result is probably caused the increase in white matter volume peaking at age 43 and declining thereafter ^[35]. Moreover, during early adulthood, the brain may enhance functional connectivity via the hippocampal–neocortical loop in response to the accumulation of visual memory (e.g., famous faces). In summary, white matter structure and memory circuits jointly promote the global integration of functional connectivity.

4.2.   Whole-brain and regional-scale differences in NVMC across face stimuli

At the node level, we obtain similar brain activation regions, metrics of NVMC (CBF vs. FCS), and network topology metrics for different visual stimulus task states (famous, unfamiliar, and scrambled).

We consider the process of brain activation as the region where neural activity causes the FCS to increase and found that brain regions such as the caudal frontal middle gyrus, cuneate lobe, fusiform, inferior temporal gyrus, and parahippocampal gyrus undergo significant changes, suggesting that the execution of the different visual task states is associated with the corresponding brain regions of attentional control, visual stimulation, and face recognition. Neural networks associated with spontaneous activity related to primary visual areas (PVA) include visual association areas, involving precuneus, precentral/postcentral gyrus, middle frontal gyrus, fusiform, inferotemporal/medial gyrus, and parahippocampal gyrus ^[36]. Considering the functions of these regions, we obtain that the PVA-related spontaneous activity may be associated with memory-related mental imagery or the process of visual memory recognition consolidation. Furthermore, the activation of the corresponding caudal middle frontal in all three visual stimulus task states suggests that the frontal mid-session caudate is significantly activated for the execution of the functional task, indicating that during this process, the subject performs memory mobilization to retrieve whether the face stimulus that is seen is familiar or not. The performance of this activation is similar to the N-back task ^[37,38].

The hippocampus is primarily involved in the processing of memory ^[39], and the surrounding parahippocampal gyrus has been found to respond to visual stimulations. The study has shown the separation of the fusiform in characterizing the physical features of self-faces and identity determination, i.e., the left fusiform distinguishes self-faces with different features, and the right fusiform distinguishes the identity features of the individuals with the faces ^[40]. In this study, the activation of the fusiform occurs on both sides, and the activation of the right fusiform is greater than that of the left side, suggesting that the brain is more inclined to differentiate and determine the identity of faces under different visual stimulations of faces.

Our results are consistent with the literature ^[41,42] that the fusiform is more responsive to faces than to scrambled faces under different visual stimulations. In previous studies ^[40,43], NVMC mainly uses imaging data to extract the corresponding indexes for measurement. Here, for the first time, we use a noninvasive large-scale network model to analyze the differences in the degree of NVMC in the task state. CBF reflects changes in neural activity laterally based on the assumption of neural-vascular coupling. When neural activity is enhanced, this brain region consumes more oxygen and has a higher metabolic demand, leading to increased perfusion. FCS describes the whole-brain functional connectivity of each brain region from a global perspective, reflecting the information-transferring role of each brain region in the whole-brain network. Consistent with previous studies ^[44,45], healthy subjects show significant cross-brain region correlations between CBF and FCS, representing normal NVMC. Under different visual stimulations, it was found that the more familiar faces results in an approximately lower index (CBF-FCS correlation), indicating a lower degree of overall NVMC across the brain. The following reasons are probable reasons why this occurs: First, the neural response gradually weakens (i.e., habituation) when subjects are exposed to familiar faces ^[46,47,48], as Grill-Spector et al. (2006) found that the activation of brain regions responsible for face recognition (e.g., the fusiform) is significantly reduced during repeated presentations of familiar faces ^[42]. Second, among the three visual stimulations (scrambled, unfamiliar, and famous), as faces become more familiar, the neural representations (such as ERP) of the brain become more stable, requiring fewer neural resources for face recognition, less activation, and weaker hemodynamic responses ^[7,46].

For the other metric (CBF/FCS), which is analyzed at the brain region level, it was found that subjects have increased CBF/FCS ratios in the transverse temporal gyrus, inferior temporal gyrus, and lateral occipital lobe, where the lateral occipital lobe is responsible for the early perception of facial features, the inferior temporal gyrus is a core brain region for visual analysis, and the transverse temporal gyrus acts as a peripheral system realizing the processing for the auditory cortical regions ^[5]. Moreover, another explanation for the activation of the transverse temporal gyrus is that one sensory or cognitive stimulus leads to the involuntary activation of a second sensory or cognitive pathway called synaesthesia ^[49]. The above results demonstrate no significant differences in CBF over FCS metrics between the three sets of visual stimulations, suggesting that activation of brain regions is similar and less variable across face stimulations in normal human subjects, but the results of Wakeman and Henson's (2015) study demonstrated that compared to scrambled faces (scrambled), faces (famous + unfamiliar) had significant activation in the fusiform, and similar results were obtained in our study (i.e., Figure 8) ^[18]; however, this result does not indicate that other regions did not undergo activation during face stimulation, nor does it indicate that faces are more likely to be activated than scrambled faces.

Figure 8.  There is a significant difference between famous + unfamiliar and scrambled on the group level metrics of NVMC (CBF/FCS) (P < 0.001).

DownLoad: Full-Size Img PowerPoint

Similar to the above results, no significant differences were found in the global network metrics calculation results. Different key nodes under the three sets of visual stimulations are demonstrated in the local betweenness centrality. Betweenness centrality measures the importance of nodes, and for key nodes, it is important in terms of connectivity and functionality of the network. Similar to the above results, no significant differences were found in the global network metrics calculation results. The left precentral gyrus, the right supramarginal gyrus, and the right lateral occipital lobe are the key nodes under the famous face stimulations. Moreover, the high values of betweenness centrality in the precentral gyrus and supramarginal gyrus are used to prevent subjects from losing attention in the task state, and subjects are asked to press one of the two buttons (based on the degree of symmetry that they had to each image) with either their left hand or their right index finger ^[18]. The lateral occipital, as the visual cortex, is undoubtedly a key node under visual stimulation tasks. The left inferior frontal gyrus orbital area and the right internal olfactory area are critical brain regions under unfamiliar faces stimulation for the inferior frontal gyrus orbital area, which is activated in emotional tasks (a task in which faces make judgments) ^[50]. Making recognition decisions will also inevitably cause activation of this brain region. The medial olfactory region serves as an important input-output region of the hippocampus ^[51], and the internal olfactory area helps to match and integrate this visual information with information stored in memory under stimulation of familiar faces; under stimulation of scrambled faces, the key nodes are the right inferior parietal lobe, the right lateral occipital lobe, and the inferior parietal lobe, which under this stimulation, do not cause activation of emotional cognition, and only the visual cortex is important in the activation process.

4.3.   Heterogeneity of cognitive functions in different subnetworks

To accurately simulate the activity of each brain region, we optimize the parameters of the model. In studying the mechanism of NVMC under visual stimulation, we found that the optimized parameters A (average excitatory synaptic gain), B (average inhibitory synaptic gain), and C (synaptic gain between pyramidal cells and interneurons) differed significantly among cognitive subnetworks. Figure 7 shows the relationship of the three parameters among the sub-networks. It was found that the three parameters have a similar pattern under different visual stimulations, and all three sets of parameters are higher under the attention network and default mode network, followed by the auditory sub-network and visual network. Furthermore, we found no significant differences in the parameters between the three stimulation tasks, but significant differences between the six sub-networks, suggesting that the activity of these networks in processing visual information, integrating situational memory, and modulating attention is significantly enhanced ^[52]. This enhanced synaptic gain improves the efficiency of communication and the strength of connections between neurons, which enhances the brain's ability to process visual information and related cognitive functions. Understanding these mechanisms can help explain individual differences in performance on visual tasks and provide new ideas for the diagnosis and treatment of neurological disorders. For example, patients with Alzheimer's disease usually show a decline in facial recognition in the early stages. Early diagnosis of Alzheimer's disease can be aided by detecting patients' neural responses under different face tasks, especially by EEG and fMRI. In addition, patients with difficulty in face recognition can be helped to improve the areas of the brain responsible for visual cognition and social interaction by designing cognitive training tasks based on face stimulation. Through repetitive training in facial recognition, brain plasticity can be promoted, thus improving the cognitive function of patients.

The parameter optimization results revealed that the value of parameter A is smaller than that of parameter B. Furthermore, the values of both optimized parameters B and C are slightly smaller than those used previously, suggesting that the neural network maintains excitatory-inhibitory balance, processes information efficiently, avoids over-suppression, enhances task-specificity, and maintains network homeostasis by precise fine-tuning and moderate optimization when processing visual information ^[53,54]. These mechanisms ensure the stability and efficient functioning of neural networks in complex visual tasks, reflecting the flexibility and adaptability of the brain in response to different cognitive tasks.

4.4.   The clinical applicability of the model

We predominantly focus on normal individuals engaged in visual stimulus task conditions; however, the proposed modeling approach demonstrates potential cross-disease applicability.

From a methodological standpoint, the model's strengths are threefold. First, it employs a heterogeneity assumption across brain regions, utilizing a simulated annealing algorithm for parameter optimization, thereby aligning more closely with the brain's physiological traits. Second, the model integrates visual task states with brain cognitive subnetworks, enabling the exploration of subnetwork metric differences under task conditions. Finally, the integration of neurophysiological and imaging data lays the foundation for the interpretability of task-state NVMC mechanism research.

In clinical practice, similar methods have primarily been applied to normal individuals in resting states. For example, Wang et al. utilized a cortical large-scale dynamic circuit model to investigate recurrent connectivity between regions, correlating with microscopic cortical level metrics, thus offering insights into the macro-cortical hierarchy ^[55]. In the context of other brain disorders, Gustavo Patow et al. examined the effects of introducing Aβ and tau proteins on local node excitation-inhibition balance in Alzheimer's disease (AD), elucidating the dynamic and complex interplay between these proteins ^[56]. In epilepsy research, numerous researchers have simulated macroscopic whole-brain epilepsy activities, though without integrating hemodynamics ^[57,58,59].

In summary, the above applications predominantly target healthy or pathological resting states, with model inputs often resembling noise-like neural activity, which does not fully capture brain state dynamics. Our approach, applied to task states, incorporates real neuro-electrophysiological and imaging data and accounts for mechanisms of NVMC, thereby facilitating a more feasible interpretation of brain cognitive activities. Moving forward, this method can be extended to various task states in normal individuals and to patients with AD, epilepsy, etc., though it is crucial to consider the unique physiological mechanisms of different diseases (the neural mass types can be further enriched based on specific diseases in the RNMM, and quantifying pathogenic biomarkers (e.g., amyloid-beta) effects on vascular reactivity), the specificity of individual patients (diffusion tensor imaging (DTI) was introduced into the model for computing the matrix of personalized structural-functional fusion matrices), and their compatibility with the model (validation results using independent datasets).

4.5.   Limitation

Dynamical network models of the whole brain are effective tools for simulating brain dynamics and understanding brain disease mechanisms ^[13]. Up to now, there are no general rules for selecting specific choice types and mathematical equations when constructing multiscale models, and specific models are highly dependent on the specific problem and the type of data available ^[60]. Moreover, the DK brain template with 68 brain regions is used in this study. In addition to the AAL, Brodmann, and Hagmann templates, Proix et al. argued that, under certain assumptions, an atlas containing about 140 brain regions might produce good agreement with empirical data ^[61]. However, the correct selection of brain region templates depends on many factors, and the optimal method of dividing brain regions is not clear. Finally, the effects of structural terms are not considered when constructing brain networks. In our future work, on the one hand, we optimize the parameters in the model more extensively to improve the relevance of the model further; we explore whether the choice of brain partitioning will affect the whole-brain model; and we add structural matrices and consider the effect of the structural-functional fusion matrices to make the model more complete.

5.   Conclusions

In this study, we integrated multimodal data into a macroscopic whole-brain model to systematically explore and characterize heterogeneous neural mechanisms underlying visual stimulation in healthy individuals. This study furthers the understanding of NVMC processes in face-specific stimulus processing and provide theoretical foundations for personalized cognitive neuroscience research.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This study was supported by the National Key Research and Development Program of China (Grant Nos. 2021YFA1000200 and 2021YFA1000202) and the National Natural Science Foundation of China (Grant Nos. 12102014, 32271361 and 12202022).

Conflict of interest

The authors declare there is no conflict of interest.