A multimodal feature-constrained whole-brain dynamics modeling and parameter optimization method

Abstract

The present invention discloses a multimodal feature joint constraint whole-brain dynamics modeling and parameter optimization method, which is applied to the fields of computational neuroscience and brain science. The current whole-brain dynamics model has deficiencies in synchronously fitting multimodal scalp EEG signals and cortical and subcortical BOLD signals across spatiotemporal scales, which limits the model's ability to comprehensively and accurately characterize the true dynamic characteristics of brain activity at different scales. To overcome this shortcoming, the present invention, based on the whole-brain dynamics model and the parameter space rapid optimization method, incorporates the core element of multimodal feature joint constraint and designs a modeling method and parameter rapid optimization strategy for the multimodal feature joint constraint whole-brain dynamics model. By using a high-dimensional parameter rapid optimization algorithm, the whole-brain dynamics model's ability to synchronously fit the brain's multimodal spatiotemporal characteristic signals is further enhanced. This has important significance and reference value for developing a more biologically plausible and interpretable whole-brain dynamics modeling method and parameter optimization system.

Description

Technical Field

[0001] The present invention belongs to the field of computational neuroscience and brain science, and particularly relates to a whole-brain model construction technology. Background Art

[0002] Whole-brain modeling (WBM), a key branch of computational neuroscience, aims to approximate neural activity throughout the entire brain through numerical computation. WBM simulates neural activity at a macroscopic scale and can be directly compared with data observed from the human brain using non-invasive neuroimaging techniques (such as MEG, EEG, and fMRI). The brain is a complex dynamical system composed of hundreds of regions. The WBM paradigm views the brain as an interconnected network of regions, with nodes representing distinct structural / functional areas and edges defined by the presence and weights of anatomical connections ("connectomes"). Commonly used neural dynamics models include neural oscillation models, neural population models, and mean-field models. By coupling the anatomical connectivity of multiple brain regions, neural dynamics, and the interactions of external stimuli and perturbations, WBM recreates multi-scale neural activity and the collective activity of neural networks. Establishing whole-brain neural computational models that integrate multi-regional interactions for specific brain functions and states is a new informatics approach and a new approach for studying brain information processing mechanisms and intervening in brain diseases at a whole-brain level.

[0003] Traditional WBM model construction is generally categorized into two types: quasi-inversion and inversion models. Quasi-inversion models combine empirical parameters from biological experiments (physiologically meaningful parameters measured experimentally) with real connectome information coupled to the neural dynamics of various brain regions to numerically simulate whole-brain neural activity. In recent years, with the continuous advancement of neuroimaging technology, prior information on structure-function coupling, such as brain structure and functional connectivity matrices, has been obtained using MRI data. Models that improve model fitting accuracy by iteratively optimizing model parameters are known as inversion models. However, whether quasi-inversion or inversion, whole-brain models have limited ability to integrate multimodal brain information. This is limited not only by the quality stability of data acquisition and variability across individuals and populations, but also by the inability to simultaneously fit cross-modal data. Furthermore, the model parameter space after inversion iterative optimization is single and does not accurately reflect the true dynamic parameter space of the individual, meaning that the model parameters lack individual specificity. Therefore, further design of modeling strategies and parameter optimization methods that can effectively integrate the spatiotemporal characteristics of multimodal brain data and enhance the accuracy of cross-modal data fitting is needed.

[0004] The researchers found that the brain's spontaneous neural activity signals are not isolated, that is, there is a cross-temporal and spatial correlation between blood oxygen-dependent signals (BOLD signals in the fMRI modality) and electroencephalogram signals (EEG signals). The BOLD signal reflects the slower blood oxygen level fluctuation process in brain activity, which mainly corresponds to low-frequency neural activity; while the EEG signal can capture fast electrophysiological activity processes. The different frequency bands of its signals represent different rhythmic activities of the brain. The data of the two modalities jointly characterize the multi-scale dynamic characteristics of brain activity. Although these signals differ in spatiotemporal resolution and physiological processes, studies have shown that they do have synchronization and cross-scale coupling relationships. However, the current whole-brain dynamics model is not sufficient to complete the synchronous fitting of multimodal signals across spatiotemporal scales, making it difficult for the model to fully and accurately reproduce the true dynamic characteristics of brain activity at different scales. Summary of the Invention

[0005] To solve the above technical problems, the present invention proposes a multimodal feature-jointly constrained whole-brain dynamics modeling and parameter optimization method, which efficiently and synchronously fits signal features at different spatiotemporal scales while rapidly optimizing the key dynamic parameter space of the model.

[0006] The technical solution adopted by the present invention is: a multimodal feature-jointly constrained whole-brain dynamics modeling and parameter optimization method, comprising:

[0007] S1. Acquire brain image and electrophysiological multimodal data of a subject to be processed and perform preprocessing; the preprocessing result at least includes an empirical multimodal signal;

[0008] S2. Construct a whole-brain dynamics model. Specifically, different functional brain regions are considered nodes in a network. Each node is described by a set of neural dynamics equations describing the average level of collective activity of the excitatory-inhibitory neural population in that brain region.

[0009] S3, obtaining a simulated multimodal signal based on the whole-brain dynamics model constructed in step S2;

[0010] S4. Optimize the fitting performance of the whole-brain dynamics model to the empirical multimodal signal by setting a weighted loss function, and continuously iterate the spatial distribution of the dynamic parameters of the whole-brain dynamics model until the loss function converges to a stable value.

[0011] Beneficial effects of the present invention: Based on the whole-brain neural dynamics computational model, the present invention proposes a multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization method, which provides an important idea for the development of a more biologically reasonable and interpretable whole-brain model modeling method and parameter optimization system. By combining the cross-modal features of the BOLD signal in the source space (based on brain functional partitioning) and the EEG signal in the head surface space (based on scalp electrode partitioning), the model significantly improves the accuracy of simultaneously fitting multi-scale spatiotemporal feature brain signals. This fills the gap in the traditional whole-brain modeling system based on inversion / quasi-inversion in that it is insufficient to simultaneously fit multi-modal signal features and lacks individual model parameter specificity. By flexibly designing the weighted loss function of the joint constraint conditions, the optimization performance of the model in the high-dimensional parameter space significantly improves the singleness and limitations of the corresponding parameter space of the traditional model fitting. The high-dimensional parameter space obtained by the rapid optimization method is more individual-specific, which lays a good foundation for creating a more reasonable and personalized model with a specific parameter space. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 It is a whole-brain neural dynamics model structure based on rapid parameter optimization.

[0013] Figure 2 It is a process of optimizing specific high-dimensional parameter space and joint feature constraint model.

[0014] Figure 3 It is the result of multimodal feature joint constraint model construction and parameter optimization. DETAILED DESCRIPTION

[0015] To facilitate those skilled in the art to understand the technical content of the present invention, the present invention is further explained below with reference to the accompanying drawings.

[0016] Based on the theoretical system of whole-brain dynamics computational models, the present invention uses the preprocessing results of cross-modal feature data at different spatiotemporal scales as the input and key constraints of the whole-brain dynamics model. The classic neural group model Wilson-Cowan is used as the basic computational unit of neural dynamics for a single network node, and the whole-brain dynamics activity is coupled through a structural connection matrix. First, the initial parameter space search boundary of the rapid optimization algorithm is defined. Then, neural signals of different modalities are combined as feature constraints and added to the model optimization process. A weighted loss function is flexibly defined, and the model's high-dimensional dynamic parameter space is continuously iterated until the loss function converges to a stable fitting value. Therefore, the specific parameter space corresponding to the global optimal solution under the convergence condition is used as the optimization parameter of the model. The optimized parameters are further used to calculate the neural signals of whole-brain activity and evaluate the fitting accuracy of the model and real data for simultaneously fitting multimodal features. Ultimately, a whole-brain model construction method and parameter optimization system based on the joint constraints of multimodal spatiotemporal features are formed.

[0017] like Figure 1 As shown, the whole-brain neurodynamic model involved in the present invention is composed of an input module, a whole-brain neurodynamic model module, and an output module, which are specifically as follows:

[0018] 1. Input module: The multimodal data measured experimentally are used as input and necessary constraints for constructing a whole-brain dynamics model.

[0019] Multimodal data includes multimodal brain imaging and electrophysiological data; multimodal brain imaging includes dMRI, fMRI, and T1w; electrophysiological data is EEG signals.

[0020] dMRI is a diffusion magnetic resonance imaging technique used to assess brain network topology. Using fiber tract tracking techniques, tools such as FSL, MRTrix3, and ANTS are used to register dMRI images to a standard space (such as MNI space). Further denoising, artifact removal, and brain tissue segmentation are performed to obtain a standardized dMRI. The brain is then partitioned according to functional brain regions defined by brain atlases of varying precision (such as the AAL3 Atlas). Finally, a probabilistic tracking algorithm is used to derive a structural connectivity matrix corresponding to the density of white matter fiber tracts between brain regions.

[0021] fMRI is a functional magnetic resonance imaging technique used to assess the blood oxygenation level dependence (BOLD) of activity at the voxel level within functional brain regions, reflecting the blood oxygenation status of local brain activity. Using Dpabi and SPM12, the BOLD signals of all voxels in the fMRI image are preprocessed and registered to a standard space. The BOLD signals of all voxels within different functional regions are averaged to obtain the average BOLD signal across the brain regions. The Pearson correlation of these averaged BOLD signals across the brain regions is then calculated to obtain a functional connectivity matrix reflecting the correlations between brain region activity.

[0022] T1w is a T1-weighted image that clearly displays brain tissue structures, such as gray matter, white matter, and cerebrospinal fluid. By using this clear brain tissue structural information as the boundaries and tissue stratification basis for individual mesh head models, a mesh head model composed of millions of finite element triangles can be constructed. This is used to construct a conversion relationship between the contribution of source dipole activity to the superposition of cerebral cortical potentials. By using different regions of the brain atlas as the location of sources in the head model, the contribution of the source dipole activity information to the superposition of potentials in the corresponding region of the cerebral cortex measured by a single electrode channel of the EEG cap can be calculated, namely the transfer matrix Leadfield.

[0023] EEG signals are a neuroimaging technique used to record cortical electrical activity. Compared to magnetic resonance imaging (MRI), they offer higher temporal resolution and can record millisecond-level cortical activity. However, the measurement process is subject to significant artifacts and noise due to factors such as head and eye movements, which can affect the representation of true neural activity. Therefore, all channel signals must be filtered, re-referenced, normalized, downsampled, and independent component removed to obtain a pre-processed EEG signal that reflects the individual's true cortical electrical activity. Specifically, using the EEGLAB tool, the EEG signal is filtered to a bandpass frequency band (0.5-45 Hz, encompassing the delta segment (0.5-4 Hz), the theta segment (4-8 Hz), the alpha segment (8-13 Hz), the beta segment (13-30 Hz), and the low-frequency gamma segment (30-45 Hz). The signal is then re-referenced, Z-score normalized, downsampled, and independent component removed to obtain a standardized EEG signal. The coherence functional connectivity matrix (Coherence) is then calculated based on the cross-power spectrum of the standardized EEG signal. This coherence matrix represents the functional connectivity matrix of the EEG modality.

[0024] 2. Multimodal feature joint constraint whole-brain dynamics model module

[0025] This module includes two parts: the construction of a whole-brain neural dynamics model and the rapid optimization of high-dimensional parameter space based on the joint constraints of multimodal features.

[0026] 21. The process of constructing a whole-brain neurodynamic model is as follows: different functional areas of the whole brain are regarded as nodes in the network. Specifically, based on a brain map of a certain precision in a standard space (such as MNI space) (such as AAL3 Atlas), the whole brain is divided into brain network nodes according to different functional areas, and each brain functional area is regarded as a node in the network. Each node is described by a set of neurodynamic equations to describe the average level of collective activity of the excitatory-inhibitory neural population in that brain region. The average activity ratio of the neural population is controlled by the Sigmoid function. Sigmoid is a nonlinear activation function widely used in machine learning and neuroscience. It can control the probability of the neural population firing within the range of [0,1] to characterize the excitatory activity firing ratio of a neural population per unit time. The neural population activity process based on feedback inhibition can use the Sigmoid activation function to control the activity ratio of excitatory and inhibitory neural populations.

[0027] 22. Rapid optimization of high-dimensional parameter space based on multimodal features and joint constraints

[0028] The structural connectivity matrix characterizes the long-range projection coupling of excitatory neural groups across the entire brain. Gaussian white noise is added to simulate the noisy nature of information transmission in the nervous system. Numerical calculations of the model generate simulated neural activity signals for each brain region. These signals are then converted into blood oxygen level-dependent (BOLD) signals using a blood oxygen balloon model.

[0029] For the simulation of EEG signals, the simulated neural signals of each brain region are superimposed on the potential calculated by the transfer matrix Leadfield to obtain the EEG lead signals of the head meter.

[0030] By setting a weighted loss function, the degree of fitting of the signal at different scales and the optimization center of gravity are controlled, and the parameter space distribution is continuously iterated until the loss function converges and stabilizes.

[0031] 3. Output module:

[0032] The individualized specific parameter space obtained by the rapid parameter optimization algorithm is input into the neural dynamics model to generate neural activity signals that simultaneously fit the different spatiotemporal characteristics of BOLD and EEG.

[0033] like Figure 2 As shown, the present invention is a multi-modal feature joint constraint whole-brain dynamics modeling and optimization method, the specific implementation process is as follows:

[0034] S1, MRI, EEG multimodal data preprocessing modeling input

[0035] By preprocessing the real multimodal data measured in individual subjects' experiments, we obtained the model input structural connection matrix and the average BOLD signal of each brain region voxel, the functional connection matrix obtained by calculating the Pearson correlation of the average BOLD signal of each brain region voxel, the Leadfield transfer matrix for converting source space and head space signals, and the Coherence consistency matrix based on the cross power spectrum based on the standardized EEG signal.

[0036] S2. Constructing a whole-brain dynamics model and setting the boundaries of the optimized parameter space

[0037] Considering brain regions as network nodes, the neural oscillation dynamics of each brain node is characterized by the Wilson-Cowan neural oscillator. As shown in the following formula, where E j (t) and I j (t) is the ratio of excitatory and inhibitory neural group activities in brain region j. e_max and S i_max are the control functions for the maximum firing ratio of the excitatory and inhibitory neural groups, respectively. E and τI is the time constant, c1, c2, c3, and c4 are the parameters of the synaptic connections between excitatory and inhibitory neural groups in the brain region (c1 excitation-excitation, c3 excitation-inhibition, c2 inhibition-excitation, and c4 inhibition-inhibition), and c5 is the global coupling strength. The structural connectivity matrix A is a single-scan diffusion structural imaging dMRI data, which is segmented into several brain regions using the AAL3 brain atlas in the standard MNI space and normalized by region size. It is used to couple the excitatory long-range connection projection neural activity of other brain regions in the whole brain. In addition, noise and time delay are inherent characteristics of the nervous system, and the transmission delay of neural activity between brain regions The average level can be defined as about 10 m / s based on the transmission speed v of neural signals in white matter fiber bundles. In this embodiment, the physical distance between nodes in different brain regions can be converted to The transmission delay of neural signals is approximately 0.8-14.8 milliseconds. j (t) is modeled as an external stimulus, w j As the additive noise of the system, it follows a standard normal distribution with a noise intensity of σ. Therefore, considering the Wilson-Cowan whole-brain dynamics model with external stimulus input, noise, and delay characteristics, the dynamic evolution process of a single node j is described as follows:

[0038]

[0039] S e / i (x) is a nonlinear activation function based on Sigmoid. x represents the network activity input of the excitatory / inhibitory neural group of the brain node, which is used to control the ratio of the network node to convert the activity input into the activity of the neural group. For example, The x is used to refer to a e / i is the excitatory / inhibitory activity gain of the activation function, θ e / i It is the conversion threshold of excitatory / inhibitory neural group activity. By adjusting the synaptic parameters c1, c2, c3, and c4, the population is made to reach the excitation-inhibition (EI) balance, and then the model is controlled to achieve a specific neural oscillation mode. By adjusting the membrane time constant τ of the neural group E and τ I , global coupling strength c5, synaptic parameters c1, c2, and the dynamic parameter space composed of can drive the node neural group to produce different activity characteristics (such as criticality, limit cycle oscillation, bifurcation, fixed point, etc.).

[0040] Since the collected individual subjects are multimodal data in a resting state, that is, the EEG activity characteristics are expressed as α segment, and the frequency corresponding to the α segment is 8-13Hz; therefore, the whole-brain dynamics model should also generate brain activity rhythms in the α frequency band. Therefore, the parameters including the physiological significance of the α rhythm are set as follows: τ E =8ms, τ I =8ms, c1=16, c2=12, c3=15, c4=3, a e =1.3,a i =2,θ e =4,θ i =3.7,P j (t) = 1.25, v = 10, σ = 0.00001, the total simulation duration was aligned with the time of the empirically measured BOLD signal, and the simulation step size was set to 0.1 ms.

[0041] In order to make the whole-brain dynamics model produce neural signals that better fit multimodal characteristics, the present invention uses a fast parameter optimization algorithm (such as genetic algorithm, Bayesian optimization) to improve the fitting performance of the whole-brain dynamics model and sets the parameter space search range according to the α rhythm physiological parameters.

[0042] Since the whole-brain dynamics model works in a high-dimensional parameter space, that is, it is controlled by multiple parameters to produce specific activities; and the extremely high-dimensional parameter space is prone to produce local optimal solutions, therefore, according to the dynamic characteristics of Wilson-Cowan, the membrane time constant τ that affects the frequency of neural oscillation is mainly optimized E and τ I ([10,20]ms), the excitatory-inhibitory neural group synaptic connection strengths c1 and c2 ([5,25]) that control the synchronization characteristics of the neural group and feedback inhibition, and the global coupling strength c5 ([0.5,3.0]). The value spaces of these parameters constitute the dynamic parameter space for the model to synchronously fit the multimodal characteristics.

[0043] S3. Simulating multimodal neural signals using random initial model dynamic parameter space

[0044] The fast optimization algorithm is suitable for situations where the objective function is computationally expensive or the parameter space is high-dimensional. The algorithm first initializes the model dynamic parameter space based on the Gaussian process (GP) and constructs a prior distribution of the approximate objective function f((θ)) as a proxy model, where θ represents the model dynamic parameter space. The Gaussian process can provide a mean prediction μ((θ)) and uncertainty estimate σ for each parameter point. 2 (θ), that is:

[0045] f(θ)~GP((μ(θ),σ 2 (θ)

[0046] To effectively balance exploration and exploitation, we define an acquisition function ψ(θ). Common choices include expected improvement (EI) or upper confidence bound (UCB). The acquisition function determines the location of sampling points by combining the mean and uncertainty of the surrogate model. The goal is to maximize the probability of improving the objective function. For example, the formula for EI is:

[0047] ψ(θ)=E[max(f(θ)-f * ,0)]

[0048] where f * is the current optimal value,

[0049] The next sampling point is selected based on the acquisition function, and a true evaluation is performed on the objective function f(θ). The results are fed back to the proxy model to update its parameter distribution. By continuously iteratively updating the proxy model and acquisition function, the optimization process gradually approaches the global optimal solution. The parameter space corresponding to the global optimal solution can be further input into the neural dynamics model computational model to synchronously fit the neural signals under the joint constraints of multimodal features.

[0050] S4. Calculate the joint feature constraint loss function

[0051] The loss function based on multimodal feature constraints mainly includes the weighted loss functions of the two parts of the BOLD signal and the EEG signal. f1 is the loss function of the model's feature fitting of the BOLD signal in the source space, which mainly includes the Pearson correlation, mean square error and probabilistic metastable sub-state space distance of the functional connectivity matrix calculated by the BOLD signal. The model's fitting of the BOLD signal is better evaluated by calculating indicators in the dynamic brain state space and the static space brain network situation. f2 is the EEG signal feature fitting loss function simulated by the model in the head space, which mainly includes the calculation of the functional connectivity of the EEG channel signal based on the cross-power spectrum Coherence in the five EEG rhythm frequency bands (δ, θ, α, β, and low-frequency γ band), and further calculation of the Pearson correlation and mean square error of the functional connectivity.

[0052] In addition, in terms of signal power spectrum feature fitting, the power spectral density (PSD) of the model EEG signal and the empirical EEG signal is calculated based on the Pwelch algorithm, and the morphological characteristics and relative power distribution of the PSD are measured by calculating the Pearson similarity and KL divergence of the PSD.

[0053] The joint feature constraint loss function f = f1 + f2 is shown below. The loss function f represents the joint constraint of the model on the synchronous fitting of the source space and head space signals.

[0054] f1=argmin<a*(1-FC_Correlation)+b*FC_MSE+c*PMS_KLD>

[0055] f2=argmin<a*(1-Coherence_Correlation)+b*Coherence_MSE+c*Fitting_PSD>

[0056] Where FC_Correlation is the Pearson correlation of the functional connectivity matrices corresponding to the simulated and measured BOLD signals, FC_MSE is the mean squared error (PMS_KLD) of the functional connectivity matrices corresponding to the simulated and measured BOLD signals (based on the activity correlation matrix obtained by processing the measured data), and PMS_KLD is the probabilistic metastable substate space distance between the simulated and measured BOLD signals. f2 represents the multiscale feature fitting index for the EEG signal modality, where Coherence_Correlation is the Pearson correlation of the coherent functional connectivity matrices (Coherence) of the simulated and measured EEG signals, Coherence_Correlation is the mean squared error (PMS_KLD) of the coherent functional connectivity matrices (Coherence) of the simulated and measured EEG signals, and Fitting_PSD is the power spectrum fitting value of the simulated and measured EEG signals. a, b, and c are the weighted values of the multiscale features.

[0057] The calculation formula of the probabilistic metastable substate space distance is as follows:

[0058]

[0059] P emp and P sim The brain state space obtained by k-means clustering of the simulated BOLD signals and the measured BOLD signals.

[0060] The calculation formula of the consistency functional connectivity matrix is as follows:

[0061]

[0062] The consistency functional connectivity matrix is used to analyze the synchronization of signal pairs in the frequency domain, where S AA (f) and S BB (f) is the autopower spectral density of channels A and B, S AB (f) is the cross power spectral density of channel A and B signals.

[0063] The process of calculating the correlation coefficient and mean square error for a functional connectivity matrix (such as FC and Coherence) is as follows:

[0064] The Pearson correlation coefficient is a statistic that measures the degree of linear correlation between two variables. Its value ranges from -1 to 1, indicating a completely negative correlation to a completely positive correlation, and 0 indicates no linear correlation. Its calculation formula is as follows:

[0065]

[0066] Among them B i,t and B j,t are the observed values of two variables, and It is the sample mean of the two variables, corresponding to the calculation process of the similarity of the functional connection matrix. The variable represents the correlation between the functional activities of the nodes in the two brain regions.

[0067] The mean square error (MSE) is an indicator that measures the difference between the model's predicted value and the actual observed value. The calculation formula is as follows:

[0068]

[0069] Where Matrix sim,i and Matrix emp,i Represents the i-th edge of the empirical and simulated functional connectivity matrices (i.e., the functional activity correlation between two nodes), which is used to calculate the functional connectivity matrix of the BOLD modality signal and the consistency functional connectivity matrix of the EEG modality to calculate the difference between the model and experience.

[0070] S5. Iterative optimization of parameter space

[0071] For high-dimensional parameter space optimization, the total number of iterations and the number of seed points are first set. The number of iterations and seed points are determined based on computing resources, typically 30-50 iterations and 3-5 times the number of seed points. In each iteration, a new point is selected for evaluation based on the acquisition function. The newly acquired data point is then added to the dataset, the Gaussian process model is updated, and the surrogate function is then updated. This process is repeated until the global optimal solution for the target loss function f is found or the preset number of iterations is reached.

[0072] S6. Optimization parameter space output and model fitting accuracy evaluation

[0073] Through rapid parameter optimization, the loss function f is determined to converge to a stable value, which represents the model's closest fitting accuracy to the target. The model's dynamic parameter space corresponding to the global optimal solution during the iterative process is then used as the model input parameter to calculate a new optimal simulated neural signal that can simultaneously fit the characteristics of different spatiotemporal patterns of BOLD and EEG.

[0074] like Figure 3 As shown in the figure, the global optimal solution obtained by iteration under the multimodal joint feature-constrained whole-brain dynamics model framework and parameter optimization method is used to fit the model to the neural signal. The performance of simultaneously fitting the multidimensional features of BOLD and EEG at different spatiotemporal scales (Loss12 represents the joint feature loss function f1+f2) is significantly higher than the accuracy of model fitting using only BOLD features (Loss1 represents the single-modal feature loss function f1) or EEG features (Loss2 represents the single-modal feature loss function f2).

[0075] Those skilled in the art will appreciate that the embodiments described herein are intended to aid the reader in understanding the principles of the present invention, and it should be understood that the scope of the present invention is not limited to such specific descriptions and embodiments. Various modifications and variations are readily apparent to those skilled in the art. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention are intended to be included within the scope of the claims.

Claims

1. A multimodal feature-constrained whole-brain dynamics modeling and parameter optimization method, characterized in that: include: S1. Acquire brain imaging and electrophysiological multimodal data of the subject to be processed and perform preprocessing; The preprocessing result at least includes an empirical multimodal signal; the empirical multimodal signal is a standardized BOLD signal and a standardized EEG signal; S2. Construct a whole-brain dynamics model. Specifically, different functional brain regions are considered nodes in a network. Each node is described by a set of neural dynamics equations describing the average level of collective activity of the excitatory-inhibitory neural population in that brain region. S3. Obtain simulated multimodal signals based on the whole-brain dynamics model constructed in step S2; specifically: The whole-brain dynamics model generates brain activity rhythms in the alpha band, including parameters of the physiological significance of the alpha rhythm, which are set as follows: the total simulation duration is aligned with the time of the empirically measured BOLD signal, and the simulation step size is set; A fast parameter optimization algorithm is used for calculation, and the parameter space search range is set according to the physiological parameters of the α rhythm; Optimize the membrane time constant τ that affects the frequency of neural oscillations based on the Wilson-Cowan dynamics E and τ I , the excitatory-inhibitory synaptic connection strengths c1 and c2 that control the synchronization characteristics of the neural group and feedback inhibition, as well as the global coupling strength c5. The value spaces of these parameters constitute the dynamic parameter space of the model's synchronous fitting of multimodal characteristics; The process of fast optimization algorithm is: Initialize the model dynamic parameter space based on the Gaussian process and construct a prior distribution of the approximate target function f(θ) as a proxy model, where θ represents the dynamic parameter space obtained by the model; the Gaussian process can provide a mean prediction μ((θ)) and uncertainty estimate σ for each parameter point 2 (θ), that is: f(θ)~GP((μ(θ),σ 2 (i)) In order to effectively balance exploration and utilization, the acquisition function ψ(θ) is defined. The acquisition function determines the location of the sampling point through the mean and uncertainty of the surrogate model, and the goal is to maximize the possibility of improving the objective function. The next sampling point is selected based on the acquisition function, and a true evaluation is performed on the objective function f(θ). The results are fed back to the proxy model to update its parameter distribution. By continuously iteratively updating the proxy model and acquisition function, the optimization process gradually approaches the global optimal solution. The parameter space corresponding to the global optimal solution is further input into the neural dynamics model calculation model to synchronously fit the neural signal under the joint constraints of multimodal features. S4. Optimize the fitting performance of the whole-brain dynamics model to the empirical multimodal signal by setting a weighted loss function, and continuously iterate the spatial distribution of the dynamic parameters of the whole-brain dynamics model until the loss function converges to a stable value.

2. The method for multimodal feature-constrained whole-brain dynamics modeling and parameter optimization according to claim 1, characterized in that: Brain imaging and electrophysiological multimodal data include dMRI, fMRI, T1w, and EEG signals; The preprocessing process is: Based on dMRI registered to the MNI standard space, the whole brain is divided into different functional areas to obtain different brain functional regions. Then, a probabilistic tracking algorithm is used to obtain the density of white matter fiber bundles between brain functional regions. The streamlines of fiber bundles from one brain functional region to another are regarded as connections between nodes. The connection strength between each brain functional region is calculated based on the number and length of streamlines to obtain the whole-brain structural connection matrix. The BOLD signals of different brain functional areas in fMRI images were subjected to head motion correction, spatial smoothing, denoising, registration, and normalization to obtain the standardized BOLD signals of all voxels in each brain functional area. The standardized BOLD signals of all voxels in the brain functional area were then averaged as the average activity level of the brain functional area. The Pearson correlation of the BOLD signals between brain functional areas was calculated to obtain the functional connectivity matrix (FC) reflecting the activity correlation between brain functional areas. Based on the brain tissue structure displayed by T1w, a grid head model is constructed. Different brain functional areas are used as the locations for the source dipoles in the head model. The contribution of the activity information of the source dipole to the superposition of the potential of the corresponding cerebral cortex area measured by all electrode channels of the EEG cap is calculated. In other words, the transfer matrix Leadfield of the conversion relationship between brain functional areas and EEG channel potentials is obtained. The EEG signal is filtered, re-referenced, normalized, downsampled, and processed with independent component removal to obtain a standardized EEG signal. Based on the standardized EEG signal, the signal power spectrum and the coherence functional connectivity matrix based on the cross power spectrum are calculated.

3. The method for multimodal feature-constrained whole-brain dynamics modeling and parameter optimization according to claim 2, characterized in that: In step S2, each node is described by a set of neural dynamics equations describing the average level of collective activity of the excitatory-inhibitory neural population in the brain region, which is expressed as: Among them, t is the time variable, E j (t) and I j (t) is the ratio of the activity of the excitatory and inhibitory neural groups in brain region j at time t, S e / i (x) is the activation function based on Sigmoid, which is used to control the firing ratio of the excitation-inhibition neural groups of the network nodes; E j (t) is the ratio of the excitability of the jth brain functional area; I j (t) is the proportion of inhibitory neural group activity in the jth brain functional area; a e / i is the excitatory / inhibitory activity gain of the activation function, θ e / i is the threshold of excitatory / inhibitory neural group activity; c1, c2, c3, c4 are synaptic connection strength parameters; τ E and τ I is the membrane time constant; c5 is the global coupling strength; P j is the external stimulus; A is the structural connection coupling matrix, w j (t) and v j (t) is Gaussian white noise, and σ is the noise intensity.

4. The method for multimodal feature-constrained whole-brain dynamics modeling and parameter optimization according to claim 3, characterized in that: The simulated multimodal signal includes a simulated BOLD signal and a simulated EEG signal.

5. The method for multimodal feature-constrained whole-brain dynamics modeling and parameter optimization according to claim 4, characterized in that: The process of outputting the simulated BOLD signal from the whole-brain dynamics model is as follows: The Euler algorithm was used to numerically solve the firing ratio of excitatory neural groups in the neural dynamics model. The total duration was calculated as the duration of the empirical fMRI BOLD signal, and the simulation time step was 0.1ms. The firing ratio of the neural group corresponding to the brain functional area at each time t was obtained, which was represented as a neural activity signal. The neural activity signals of all brain functional areas are simulated by the whole-brain dynamics model, and the blood oxygen balloon model is used to convert the neural activity signals into blood oxygen level-dependent BOLD signals.

6. The method for multimodal feature-constrained whole-brain dynamics modeling and parameter optimization according to claim 5, characterized in that: The process of outputting simulated EEG signals from the whole-brain dynamics model is as follows: By transferring the Leadfield matrix, the neural signals of the brain areas simulated in the source space are converted into EEG signals based on the head surface electrode channels through the consistent functional connection matrix Coherence.

7. The method for multimodal feature-constrained whole-brain dynamics modeling and parameter optimization according to claim 6, characterized in that: The dynamic parameter space optimized for the whole-brain dynamics model includes: τ E , τ I , c1, c2, c5.

8. The method for multimodal feature-constrained whole-brain dynamics modeling and parameter optimization according to claim 7, characterized in that: The weighted loss function is expressed as: f=f1+f2 Where f is the weighted loss function, f1 represents the loss function of the whole-brain dynamics model fitting the unimodal characteristics of the BOLD signal, and f1 = argmin<a*(1-FC_Correlation)+b*FC_MSE+c*PMS_KLD> , FC_Correlation is the Pearson correlation of the functional connectivity matrices corresponding to the simulated BOLD signal and the standardized BOLD signal, FC_mSE is the mean square error of the functional connectivity matrices corresponding to the simulated BOLD signal and the standardized BOLD signal, PMS_KLD is the probabilistic metastable sub-state space distance between the simulated BOLD signal and the standardized BOLD signal, f2 represents the loss function of the whole-brain dynamics model fitting the unimodal features of the EEG signal simulated in the head space, f2 = argmin<a*(1-Coherence_Correlation)+b*Coherence_MSE+c*Fitting_PSD> , Coherence_Correlation is the Pearson correlation of the consistent functional connectivity matrix Coherence of the simulated EEG signal and the standardized EEG signal, Coherence_Correlation is the mean square error of the consistent functional connectivity matrix Coherence of the simulated EEG signal and the standardized EEG signal, Fitting_PSD is the power spectrum fitting value of the simulated EEG signal and the table-transformed EEG signal, a, b, c are the weighted values of the multi-scale features.

Images