Multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization method

Abstract

The invention discloses a multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization method which is applied to the field of computational neuroscience and brain science. Aiming at the defects of a current whole brain dynamic model in the aspect of synchronously fitting multi-modal scalp EEG signals and cortical and subcortical BOLD signals across spatial and temporal scales, the capability of the model for comprehensively and accurately representing real dynamic characteristics of brain activities on different scales is limited; in order to overcome the short plate, the core element of multi-modal feature combined constraint is added on the basis of a whole-brain dynamic model and a parameter space rapid optimization method, and a modeling method and a parameter rapid optimization strategy for the multi-modal feature combined constraint whole-brain dynamic model are designed; by applying a high-dimensional parameter rapid optimization algorithm, the capability of the whole-brain dynamic model in synchronous fitting of brain multi-modal spatial-temporal characteristic signals is further improved; the method has important significance and reference value for developing a whole-brain dynamic model modeling method and a parameter optimization system which are more biological rationality and interpretability.

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) is an important branch of computational neuroscience. Its goal is to approximate the neural activity of the entire brain through numerical calculations. WBM realizes the simulation of neural activity on a macroscopic scale and can be directly compared with the data observed from the human brain by non-invasive neuroimaging techniques (such as MEG, EEG, fMRI, etc.). The brain is a complex dynamic system composed of hundreds of regions. The research paradigm of WBM regards the brain as an interconnected regional network. The nodes in the network represent different structural / functional regions of the brain, and the edges of the network are defined by the existence and weight of anatomical connections ("connectomes"). Commonly used neural dynamic models include neural oscillation models, neural group models, and mean field models. By coupling the anatomical connection structure of multiple brain regions, neural dynamics, and the interaction of external stimulus perturbations, WBM reproduces multi-scale neural activities, collective activities of neural networks, etc. For specific brain functions and brain states, establishing a whole-brain neural computational model of multi-brain region interaction is a new informatics method and new approach to study brain information processing mechanisms and brain disease intervention treatment from the whole brain level.

[0003] Traditional WBM model construction is usually divided into two types of whole-brain dynamic models: quasi-inversion and inversion models. Quasi-inversion models combine empirical parameters from biological experiments (experimentally measured parameters with physiological significance) and use real connection group information to couple the neural dynamics of each brain region to numerically simulate whole-brain neural activity. In recent years, with the continuous development of neuroimaging technology, the prior information of structure-function coupling such as brain structure and functional connection matrix is ​​obtained with the help of MRI data. The model that improves the model fitting accuracy by iteratively optimizing model parameters is called inversion model. However, whether it is a quasi-inversion or inversion model, the whole-brain model has very limited ability to integrate multimodal information of the brain. It is not only limited by the quality stability of data acquisition and the variability across individuals / populations, but also most of them do not have the ability to synchronously fit cross-modal data. In addition, the model parameter space after inversion iterative optimization is single and cannot accurately reflect the individualized real dynamic parameter space, that is, the model parameters are not individual-specific. Therefore, it is necessary to further design modeling strategies and parameter optimization methods that can effectively integrate the spatiotemporal characteristics of multimodal brain data and enhance the fitting accuracy of cross-modal data.

[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 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, and 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] In order to solve the above technical problems, the present invention proposes a multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization method, which can efficiently and synchronously fit signal characteristics of 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 multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization method, comprising:

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

[0008] S2. Construct a whole-brain dynamics model. Specifically, different brain functional areas of the whole brain are regarded as nodes in the network. Each node is described by a set of neural dynamics equations to describe the average level of collective activity of the excitatory-inhibitory neural population in the brain functional area.

[0009] S3, obtaining a simulated multimodal signal according to 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 signals 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 developing 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 partitions) and the EEG signal in the head surface space (based on scalp electrode partitions), the model significantly improves the accuracy of simultaneously fitting multi-scale spatiotemporal characteristic brain signals. Thereby filling the gap in the traditional whole-brain modeling system based on inversion / quasi-inversion in that it is not enough 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 individualized 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] The present invention is based on the theoretical system of whole-brain dynamics computational model, takes the preprocessing results of cross-modal feature data of different spatiotemporal scales as the input and key constraints of the whole-brain dynamics model, uses the classic neural group model Wilson-Cowan as the basic neural dynamics computational unit of a single network node, and couples the whole-brain dynamics activity through the structural connection matrix. First, the initial parameter space search boundary of the fast optimization algorithm is defined, and then the neural signals of different modalities are combined as feature constraints to join the model optimization process, the weighted loss function is flexibly defined, and the high-dimensional dynamics parameter space of the model 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, and the optimized parameters are further used to calculate the neural signals of the whole-brain activity, and the fitting accuracy of the model and the real data for simultaneously fitting the multi-modal features is evaluated. Finally, a whole-brain model construction method and parameter optimization system based on the joint constraints of multi-modal 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 in the experiment are used as the input and necessary constraints for building the whole-brain dynamics model.

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

[0020] dMRI is a diffusion magnetic resonance imaging technology used to evaluate the topological structure of brain networks. Through fiber tract tracking technology, dMRI images are registered to standard space (such as MNI space) using tools such as FSL, MRTrix3, and ANTS, and further denoising, artifact removal, and brain tissue division are performed to obtain standardized dMRI. Then, according to the brain functional partitions divided by brain maps of different precision (such as AAL3 Atlas), the structural connection matrix corresponding to the density of white matter fiber bundles between brain regions is finally obtained using a probabilistic tracking algorithm.

[0021] fMRI is a functional magnetic resonance imaging technique used to evaluate the blood oxygen level dependence (BOLD) of the activity of voxel-level sources in the functional divisions of the brain, reflecting the blood oxygen situation of local brain activity. Using Dpabi and SPM12, the BOLD signals of all voxels in the fMRI image are preprocessed and registered to the standard space. By averaging the BOLD signals of all voxels in different functional divisions, the average BOLD signal of the brain region voxels is obtained. The Pearson correlation of the average BOLD signal of the brain region voxels is further calculated to obtain the functional connection matrix reflecting the association of brain region activity.

[0022] T1w is T1-weighted imaging, which can clearly display brain tissue structures, such as gray matter, white matter, cerebrospinal fluid, etc. By using clear brain tissue structure information as the boundary and tissue stratification basis of the individual grid head model, a grid head model composed of millions of finite element triangles can be constructed to construct the conversion relationship between the activity of the source dipole and the contribution of the superposition of the cerebral cortex potential. By using different areas of the brain map as the location of the source in the head model, the contribution of the activity information of the source dipole to the superposition of the potential on the cerebral cortex of the corresponding area measured by a single electrode channel of the EEG cap can be calculated, that is, the transfer matrix Leadfield.

[0023] EEG signal is a neuroimaging technology used to record the electrical activity of the cerebral cortex. Compared with magnetic resonance imaging, it has a higher time resolution and can record millisecond-level cerebral cortical electrical activity. However, during the measurement process, a large amount of artifacts and noise will be generated due to factors such as head movement and eye movement, which will affect the representation of real neural activity. Therefore, it is necessary to filter, re-reference, standardize, downsample, and remove independent components of all channel signals to obtain pre-processed EEG signals to reflect the individual's real cortical electrical activity. Specifically: Using the EEGLAB tool, the EEG signal is filtered to a bandpass frequency band (0.5-45Hz, including δ segment 0.5-4Hz, θ segment 4-8Hz, α segment 8-13Hz, β segment 13-30Hz, low-frequency γ segment 30-45Hz), re-referenced, Z-score standardized, downsampled, and independent components removed to obtain a standardized EEG signal. Then, based on the standardized EEG signal, the consistent functional connectivity matrix Coherence based on the cross power spectrum is calculated. The consistent functional connectivity matrix Coherence 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 the 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 the brain area. 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 to be released in the range of [0,1] to characterize the excitatory activity release 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 the excitatory and inhibitory neural populations.

[0027] 22. Rapid optimization of high-dimensional parameter space based on multi-modal feature joint constraints

[0028] The structural connection matrix is ​​used to characterize the long-range connection projection coupling of excitatory neural groups in other areas of the whole brain, and Gaussian white noise is added to simulate the noisy characteristics of information transmission in the nervous system. Through the numerical calculation of the model, the simulated neural activity signals of each brain region are obtained, and the blood oxygen balloon model is used to convert the neural activity signals into blood oxygen level dependent (BOLD) signals.

[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 table.

[0030] The weighted loss function is set to control the degree of fitting of signals at different scales and the optimization center of gravity, and the parameter space distribution is continuously iterated until the loss function converges to a stable value.

[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 a neural activity signal that simultaneously fits the different spatiotemporal characteristics of BOLD and EEG.

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

[0034] S1. MRI and EEG multimodal data preprocessing modeling input

[0035] By preprocessing the real multimodal data measured by individual subjects, 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 source space and head space signal conversion, and the Coherence consistency matrix based on the standardized EEG signal, and calculated the cross-power spectrum.

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

[0037] The brain regions are regarded as network nodes, and the dynamic process of neural oscillation of each brain node is characterized by Wilson-Cowan neural oscillator. As shown in formulas (1) and (2), 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 of the maximum firing ratio of the excitatory and inhibitory neural groups, respectively. Eand τ I is the time constant, c1, c2, c3, c4 are the parameters of the synaptic connections between the excitatory and inhibitory neural groups in the brain region (c1 excitation-excitation, c3 excitation-inhibition, c2 inhibition-excitation, 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 by 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 activities 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 according to the transmission speed v of the neural signal in the white matter fiber bundle. In this embodiment, according to the physical distance of nodes in different brain regions, it can be converted into The transmission delay of neural signals is approximated in 0.8-14.8 milliseconds. j (t) is modeled as an external stimulus, w j As the additive noise of the system, it obeys the standard normal distribution and the noise intensity is σ. 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. It is used to control the network node to convert the activity input into the activity ratio 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 can reach an excitation-inhibition (EI) balance, and then control the model 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 composed of the dynamic parameter space 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 the α 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, so 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 generate neural signals that better fit multimodal characteristics, the present invention uses a fast parameter optimization algorithm (such as a genetic algorithm and 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 extremely high-dimensional parameter spaces are 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 synaptic connection strengths c1 and c2 ([5,25]) of the excitatory-inhibitory neural groups that control the synchronization characteristics of the neural groups and feedback inhibition, and the global coupling strength c5 ([0.5,3.0]). The value spaces of these parameters constitute the dynamic parameter space of the model for synchronously fitting the multimodal characteristics.

[0043] S3, random initial model dynamic parameter space simulation of multimodal neural signals

[0044] The fast optimization algorithm is suitable for situations where the objective function is expensive to calculate 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 that approximates the objective function f(θ) as a proxy model, where θ represents the model's 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] In order to effectively balance exploration and exploitation, the acquisition function ψ(θ) is defined. The general choices include expected improvement (EI) or upper confidence bound (UCB). The acquisition function determines the location of the sampling point through the mean and uncertainty of the surrogate model. The goal is to maximize the possibility of improving the objective function. For example, the formula of EI is:

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

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

[0049] The next sampling point is selected according to the acquisition function, and a true evaluation is performed on the objective function f(θ), and the result is fed back to the proxy model to update its parameter distribution. By continuously iteratively updating the proxy model and the acquisition function, the optimization process gradually approaches the global optimal solution. Using the parameter space corresponding to the global optimal solution, the neural dynamics model can be further input into the 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 features of BOLD signals and EEG signals. 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 substate space distance of the functional connection 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 brain network situation. f2 is the EEG signal feature fitting loss function simulated by the model in the head-table space, which mainly includes the calculation of the functional connection of the EEG channel signal based on the cross-power spectrum Coherence in five EEG rhythm frequency bands (δ, θ, α, β, and low-frequency γ bands), and further calculating the Pearson correlation and mean square error of the functional connection.

[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 psd are measured by calculating the Pearson similarity and KL divergence of psd.

[0053] The joint feature constraint loss function f=f1+f2 is shown as follows. The loss function f represents the joint constraints 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] Among them, FC_Correlation is the Pearson correlation of the functional connectivity matrix corresponding to the simulated BOLD signal and the measured BOLD signal, FC_MSE is the mean square error of the functional connectivity matrix corresponding to the simulated BOLD signal and the measured BOLD signal (the activity correlation matrix obtained based on the measured data processing), PMS_KLD is the probabilistic metastable substate space distance between the simulated BOLD signal and the measured BOLD signal. f2 represents the multi-scale feature fitting index for the EEG signal mode, among which Coherence_Correlation is the Pearson correlation of the coherent functional connectivity matrix Coherence of the simulated EEG signal and the measured EEG signal, Coherence_Correlation is the mean square error of the coherent functional connectivity matrix Coherence of the simulated EEG signal and the measured EEG signal, and Fitting_PSD is the power spectrum fitting value of the simulated EEG signal and the measured EEG signal. a, b, c are the weighted values ​​of the multi-scale 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 simulated BOLD signals and measured BOLD signals.

[0060] The calculation formula of the consistent functional connection matrix is ​​as follows:

[0061]

[0062] The consistent 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 auto-power 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 complete negative correlation to a complete positive correlation, and 0 indicates no linear correlation. The 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 of the functional activities between the nodes of the two brain regions.

[0067] The mean square error (MSE) is a measure of 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), used for calculating the functional connectivity matrix for BOLD modality signals and the consistent functional connectivity matrix for 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, first set the total number of iterations and the number of seed points. The number of iterations and the number of seed points are determined according to the computing resources. Generally, the number of iterations is set to 30-50 times, and the number of seed points is 3-5 times the number of iterations. In each iteration, a new point is selected for evaluation according to the acquisition function, and then the newly obtained data point is added to the data set, the Gaussian process model is updated, and then the substitution function is updated. This process is repeated until the global optimal solution of 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 converges to a stable value, which is the fitting accuracy of the model closest to the target. The model dynamic parameter space corresponding to the global optimal solution in the iterative process is further used as the model input parameter to calculate the new optimal simulated neural signal, which 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, and its corresponding parameter space is used for model fitting of neural signals. The performance in simultaneously fitting 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 (Loss1 represents the single-modal feature loss function f1) features 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 help readers understand the principles of the present invention, and should be understood that the scope of protection of the present invention is not limited to such specific statements and embodiments. For those skilled in the art, the present invention may have various changes and variations. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of the claims of the present invention.

Claims

1. A multi-modal feature-jointly constrained whole-brain dynamics modeling and parameter optimization method, characterized in that: include: S1, obtaining brain images and electrophysiological multimodal data of the subject to be processed and performing preprocessing; The preprocessing result includes at least an empirical multimodal signal; S2. Construct a whole-brain dynamics model. Specifically, different brain functional areas of the whole brain are regarded as nodes in the network. Each node is described by a set of neural dynamics equations to describe the average level of collective activity of the excitatory-inhibitory neural population in the brain functional area. S3, obtaining a simulated multimodal signal according to the whole-brain dynamics model constructed in step S2; S4. Optimize the fitting performance of the whole-brain dynamics model to the empirical multimodal signals 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. A method for whole-brain dynamics modeling and parameter optimization with multi-modal feature joint constraints 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 the dMRI registered to the MNI standard space, the whole brain is divided into different functional areas to obtain different brain functional areas; then the probability tracking algorithm is used to obtain the density of white matter fiber bundles between brain functional areas. The streamlines from one brain functional area to another are regarded as the connection between nodes. The connection strength between each brain functional area is calculated according to the number / 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 standardization 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 association of brain functional areas. According to the brain tissue structure displayed by T1w, a grid head model is constructed, and different brain functional areas are used as the locations where the source dipoles in the head model are placed. 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, that is, 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, standardized, 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 consistent functional connectivity matrix Coherence based on the cross power spectrum are calculated.

3. The method for multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization according to claim 2, characterized in that: The empirical multimodal signals are the normalized BOLD signal and the normalized EEG signal.

4. The method for multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization according to claim 3 is characterized in that: In step S2, each node is described by a set of neural dynamics equations, which are expressed as follows: 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 an activation function based on Sigmoid, which is used to control the firing ratio of the excitation-inhibition neural groups in 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.

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

6. A method for multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization according to claim 5, characterized in that: The process of outputting the simulated BOLD signal from the whole-brain dynamics model is: The Euler algorithm is used to solve the numerical solution of the firing ratio of the excitatory neural group in the neural dynamics model. The total duration is calculated as the duration of the empirical fMRI BOLD signal, and the simulation time step is 0.1ms. The firing ratio of the neural group corresponding to the brain functional area at each time t is obtained, which is 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.

7. The method for multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization according to claim 6, characterized in that: The process of outputting simulated EEG signals from the whole-brain dynamics model is as follows: Through the transfer matrix Leadfield, the simulated brain area neural signals in the source space are converted into EEG signals based on the head surface electrode channels through the consistent functional connection matrix Coherence.

8. The method for multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization according to claim 7, characterized in that: The dynamic parameter space optimized by the whole-brain dynamic model includes: τ E , τ I , c1, c2, c5.

9. The method for multi-modal feature joint constraint whole-brain dynamics modeling and parameter optimization according to claim 8, 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 substate 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-surface 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_Correnlation 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 transformed EEG signal, a, b, c are the weighted values ​​of multi-scale features.

Images