Digital twin brain reduced model fast modeling method for brain image data integration

Abstract

The application provides a fast digital twin brain modeling method based on a linear dynamics model, and the method comprises the following steps: S1, constructing an initial DTB model, using a linear dynamics oscillator to construct an oscillator network model with separable coupling terms as a linear system to be optimized; S2, calculating an empirical covariance matrix, using a covariance matrix of BOLD signals between brain regions as the empirical covariance matrix; S3, constructing a target covariance matrix, converting the empirical covariance matrix into a target covariance matrix of the current linear system based on a dynamics model prior; S4, constructing a dynamics parameter evaluation model for evaluating the rationality of the dynamics parameters; and S5, synchronously optimizing the dynamics parameters and the coupling matrix by using an optimization algorithm. The method can significantly improve the performance of the DTB model and improve the modeling efficiency.

Description

Technical Field

[0001] This invention relates to the interdisciplinary field of biomedical engineering and brain science, and in particular to a rapid modeling method for digital twin brain reduction models used for brain imaging data integration. Background Technology

[0002] The brain is a complex dynamic system composed of dynamic interactions between different functional subregions. To study the brain's cognitive functions, information processing, and disease mechanisms, it is urgent to integrate brain imaging data at multiple spatiotemporal scales and in multiple modalities. The data-driven, biologically constrained digital twin brain (DTB) model can effectively integrate multi-scale, multi-modal neuroscience data and play an important role in revealing brain information processing and disease mechanisms and advancing personalized precision medicine.

[0003] The brain connectome plays a central role in elucidating the information processing mechanisms of the nervous system and constructing digital twin brain models. Structural connectivity (SC), a commonly used technique based on diffusion magnetic resonance imaging (dMRI), provides a basic framework for understanding the physical connections between brain regions by tracing the anatomical paths of white matter fiber tracts. However, on the one hand, SC can only characterize the general features of axonal pathways and lacks information on bilateral brain region connections, thus failing to represent the true anatomical connections of axonal fiber tracts. On the other hand, SC cannot directly reflect complex coupling effects such as synaptic transmission efficiency and neurotransmitter modulation effects between neuronal clusters. This lack of information leads to structure-function decoupling: more than 40% of the spatiotemporal dynamic variability of experimentally observed blood-oxygen-level dependent (BOLD) signal time series cannot be explained by static SC.

[0004] The limitations of SC (Simultaneous Coupling) have a significant impact on digital twin brain modeling. The calibration of dynamic parameters in existing brain dynamics models (such as the Wilson-Cowan model and Hopf oscillator networks) relies on the precise setting of coupling weights. When SC is directly used as the whole-brain coupling matrix (hereinafter referred to as the "coupling matrix"), the errors between the simulation-generated oscillation patterns, synchronization behaviors, and other key dynamic features and the real EEG / MEG signals generally exceed acceptable thresholds (e.g., mean squared error > 10%). DTB models lacking realistic coupling relationships and dynamic parameters struggle to reproduce neural processes dependent on coupling dynamics in the brain, such as consciousness state transitions and seizure propagation. Therefore, it is necessary to introduce Effective Connectivity (EC) to replace SC for quantitative characterization of inter-brain coupling effects and to use heterogeneous dynamic parameters to simulate the dynamic differences between different brain regions in the real brain, thereby achieving biologically plausible model simulations.

[0005] Functional connectivity (FC) provides important clues for uncovering potential brain-related coupling (EC) by calculating statistical indices of brain region bounded-down (BOLD) signals (such as the Pearson correlation coefficient). Studies have shown a significant but highly nonlinear positive correlation between the strength of functional connectivity and coupling effects between brain regions. Therefore, decoding potential coupling effects from functional connectivity holds promise for improving the biological plausibility of brain-related brain-tightening (DTB) models. However, existing EC-based DTB modeling often uses time-complexity EC prediction algorithms such as machine learning or causal analysis, failing to simultaneously calibrate kinetic parameters during EC prediction and struggling to simultaneously guarantee the interpretability and simulation effectiveness of DTB models. These limitations severely restrict the practical application of DTB models in personalized medicine. Summary of the Invention

[0006] This invention aims to solve the problems of slow EC prediction, poor simulation results, and difficulty in synchronously calibrating dynamic parameters in the process of DTB modeling. It provides a rapid modeling method for digital twin brain reduction models for brain imaging data integration, and realizes rapid synchronous optimization of coupling matrix and dynamic parameters in DTB model.

[0007] This invention proposes a rapid modeling method for digital twin brain reduction models used for brain imaging data integration. Based on the covariance matrix of BOLD signals in brain regions and an analytically solvable whole-brain linear dynamics model, it achieves rapid calculation of the dynamic coupling matrix through approximate analytical solution and combines this with an optimization algorithm to search for dynamic parameters. The method includes the following steps:

[0008] Step S1: Construct an initial DTB model. Use linear dynamic oscillators to construct an oscillator network model with separable coupling terms as the linear system to be optimized.

[0009] Step S2: Calculate the empirical covariance matrix, using the covariance matrix of the BOLD signal in the brain region as the empirical covariance matrix;

[0010] Step S3: Construct the target covariance matrix. Based on the prior knowledge of the dynamic model, convert the empirical covariance matrix into the target covariance matrix of the current linear system.

[0011] Step S4: Construct a dynamic parameter evaluation model to assess the rationality of the dynamic parameters. The dynamic parameter evaluation model includes a dynamic matrix construction module, a coupling matrix solving module, and a loss function calculation module. The dynamic matrix construction module receives a set of dynamic parameters to be evaluated and constructs an uncoupled system dynamic matrix based on the dynamic model. The coupling matrix solving module uses the system dynamic matrix, the target covariance matrix, and the dynamic parameters to be evaluated to construct a Sylvester equation with the coupling matrix as the variable. It approximates the unique solution of the Sylvester equation by minimizing the residuals and converts it into a dynamic coupling matrix. The loss function calculation module adds the obtained coupling matrix to the DTB model, calculates the system covariance matrix and the functional connectivity matrix, and constructs a loss function using one or more statistical indicators to evaluate the simulation effect of the parameter combination and its corresponding coupling matrix.

[0012] Step S5: Use an optimization algorithm to simultaneously optimize the dynamic parameters and coupling matrix, set the optimization range of the dynamic parameters, use the loss function obtained from the above dynamic parameter evaluation model as the standard, use the optimization algorithm to quickly explore the parameter space, find the optimal dynamic parameters and their corresponding optimal dynamic coupling matrix, and construct the optimized DTB model.

[0013] Furthermore, in step S1, a linearized Stuart-Landau oscillator model is used as the nodal dynamics model, and its expression is as follows:

[0014]

[0015]

[0016] Among them, u=[u1,…,u N ] T Let A represent the model state vector and its node components, where N is the number of nodes in the DTB model, i.e., the number of oscillators in the oscillator network, δu is the oscillation of the model's state vector near the fixed point, and A is the Jacobian matrix of the 2N×2N system at the fixed point. jk Let a be the element in the j-th row and k-th column of the Jacobian matrix A. j W represents the bifurcation parameter of node j. jk Let η be the coupling weight between node j and node k in the coupling matrix W, where j = 1, ..., N, k = 1, ..., N, and η = [η1, ..., η]. 2N ] represents the vector representation of Gaussian white noise at each node, ω j Let be the eigenangular frequency of node j; specifically, the oscillator model in the oscillator network, i.e. the node dynamics model, is not limited to the linearized Stuart-Landau model, but is applicable to any oscillator network model that satisfies the condition that the coupling terms can be separated.

[0017] Matrix A can be represented in the following block form:

[0018]

[0019] Among them, A xx A xy A yx A yy All are N×N matrices and have the following relationship:

[0020] A xx =A yy =diag(aS)+W,

[0021] A yx =-A xy =diag(ω),

[0022] Where diag(v) represents a diagonal matrix with vector v as its main diagonal, and a = [a1, ..., a2]. N ] represents the vector representation of the bifurcation coefficients of each node. Let ω be the degree of each node in the coupling matrix W, where ω = [ω1, ..., ω]. N [ ] represents the vector representation of the intrinsic angular frequencies of each node;

[0023] Further, let A be represented as:

[0024] The two diagonal blocks of L are the graphical Laplace matrices of the coupling matrix W. Thus, the system Jacobian matrix A has been decomposed into the uncoupled system dynamics matrix A0 and the matrix L containing coupling information.

[0025] Furthermore, in step S2, the covariance matrix C v elements in The calculation method is as follows:

[0026]

[0027] Where T represents the total duration of the BOLD signal, B i (t) represents the BOLD signal value of brain region i at time t. This represents the average BOLD signal in brain region i.

[0028] Furthermore, in step S3, the variance of the signal at each node of the linear system is set to be proportional to the variance of the empirical BOLD signal. Then, the system's target covariance matrix is ​​expressed as: Where m is the scaling factor.

[0029] Furthermore, the specific function of the dynamic matrix construction module in step S4 is to receive a set of dynamic parameters to be evaluated, input the dynamic parameters into the dynamic model constructed in step S1, and thus construct the uncoupled system dynamic matrix corresponding to the dynamic parameters.

[0030] Furthermore, the specific function of the coupling matrix solving module in step S4 is as follows:

[0031] The Jacobian matrix A and the covariance matrix C of the system under steady state are obtained from the Lyapunov equations of the linear system. vs The following relationship must be satisfied:

[0032] AC vs +C vs A T =-Q,

[0033] Where Q = σ 2 I is the covariance matrix of Gaussian white noise, σ 2 Let I be the variance of the Gaussian white noise, and I be the identity matrix. After separating the coupling terms in the Jacobian matrix, we obtain:

[0034]

[0035] Since L is a symmetric matrix, then L T =L = diag(S) - W, let k = -L = W - diag(S), then we get:

[0036]

[0037] This equation is the Sylvester equation DK+KE=-Γ in D=E=C. vs In a special form, the Sylvester equation has a unique solution when D and -E do not have the same eigenvalues.

[0038] The system is configured to satisfy the following condition when it has the optimal coupling matrix: according to Solve for K, and from this we obtain the optimal coupling matrix W, which is the off-diagonal part of the unique solution K.

[0039] Furthermore, the loss function in step S4 is calculated as follows:

[0040] The system covariance matrix is ​​obtained by solving the Lyapunov equations of the linear system. Dividing each element of the system covariance matrix by the product of the standard deviations of the corresponding two nodes yields the simulated FC matrix FC, defined based on the Pearson correlation coefficient. sim By calculating FC sim The empirical FC matrix FC calculated based on measured BOLD signals. empUsing the KS distance (ksdist) and the Pearson correlation coefficient (corr), the joint loss function is obtained as: Loss = ksdist + (1 - corr), where the KS distance is the absolute value of the maximum difference between the empirical and simulated cumulative distribution functions of the FC matrix, and the Pearson correlation coefficient is defined as... in, Representing the empirical FC matrix FC emp The mean, To simulate the FC matrix FC sim The mean.

[0041] Finally, the optimized DTB model was applied to neuroimaging data analysis, prediction of neuromodulation efficacy, and quantitative diagnosis and subtype classification of neurological diseases. The analysis results provide doctors with auxiliary diagnostic references.

[0042] This invention proposes a rapid modeling method for digital twin brain reduction (DTB) models used for brain imaging data integration. It constructs a linear dynamic oscillator network model with separable coupling terms as the linear system to be optimized. Then, it uses the covariance matrix of the BOLD signals in brain regions to construct the target covariance matrix of the model and builds a dynamic parameter evaluation model. An approximate solution algorithm is used to calculate the unique optimal dynamic coupling matrix corresponding to the parameter combination, and the dynamic parameters and dynamic coupling matrix are optimized through a loss function. This enables simultaneous optimization of the dynamic parameters and coupling matrix of the DTB model, significantly improving the simulation effect of the DTB model while greatly reducing modeling time. The optimized DTB model can be used in personalized medicine fields such as personalized neuroimaging data analysis, prediction of neuromodulation efficacy, and quantitative diagnosis and subtype classification of neurological diseases, significantly improving the accuracy of medical information analysis and effectively promoting the auxiliary diagnosis of diseases. Attached Figure Description

[0043] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0044] Figure 1 This is a flowchart of a rapid modeling method for digital twin brain reduction models used for brain imaging data integration, provided by an embodiment of the present invention.

[0045] Figure 2 This is a schematic diagram of the eigenvalues ​​of the system target covariance matrix provided in an embodiment of the present invention;

[0046] Figure 3This is a schematic diagram illustrating the modeling effect of the method of the present invention in a group of healthy individuals' average data, provided by an embodiment of the present invention. In this diagram, (a) is a histogram of the optimal kinetic parameters, (b) is the optimal coupling matrix of this embodiment, and (c) is the empirical FC matrix of this embodiment. emp With analog FC matrix FC sim The histogram is shown in (d), which is the simulated FC matrix obtained in this embodiment.

[0047] Figure 4 This is a schematic diagram comparing the time consumption of the method of the present invention and the traditional modeling method under the same effect, provided by the embodiments of the present invention. Detailed Implementation

[0048] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0049] like Figure 1 As shown, the present invention proposes a rapid modeling method for digital twin brain reduction models for brain imaging data integration, comprising the following steps:

[0050] Step S1: Construct an initial DTB model. Use linear dynamic oscillators to construct an oscillator network model with separable coupling terms as the linear system to be optimized.

[0051] In this embodiment of the invention, a linearized Stuart-Landau oscillator model is used as the nodal dynamics model, and the expression is as follows:

[0052]

[0053]

[0054] Among them, u=[u1,…,u N ] T Let N represent the model state vector and its node components, where N is the number of nodes in the DTB model (i.e., the number of oscillators in the oscillator network), δu represents the oscillation of the model's state vector near the fixed point, and δx is a component of N. i ,δy i Let i = 1, ..., N be the fluctuations of the two state variables of node i in the model near the fixed point, and A be the Jacobian matrix of the 2N×2N system at the fixed point. jk Let a be the element in the j-th row and k-th column of the Jacobian matrix A. j W represents the bifurcation parameter of node j. jkLet η be the coupling weight between node j and node k in the coupling matrix W, where j = 1, ..., N, k = 1, ..., N, and η = [η1, ..., η]. 2N ] represents the vector representation of Gaussian white noise at each node, ω j Let be the eigenfrequency of node j. A coupled linear Stuart-Landau oscillator network model of this form can be called a linearized Hopf whole-brain model.

[0055] Matrix A can be represented in the following block form:

[0056]

[0057] Among them, A xx A xy A yx A yy All are N×N matrices, with the following form:

[0058] A xx =A yy =diag(aS)+W,

[0059] A yx =-A xy =diag(ω),

[0060] Where diag(v) represents a diagonal matrix with vector v as its main diagonal, and a = [a1, ..., a2]. N ] represents the vector representation of the bifurcation coefficients of each node. Let ω be the degree of each node in the coupling matrix W, where ω = [ω1, ..., ω]. N ] is the vector representation of the intrinsic angular frequencies of each node.

[0061] A can be further represented as:

[0062] In this embodiment, the two main diagonal blocks of L are the graph Laplacian matrices of the coupling matrix W. Thus, the system Jacobian matrix A has been decomposed into an uncoupled system dynamics matrix A0 and a matrix L containing coupling information. An oscillator network model that allows for the above-mentioned separation operation on the Jacobian matrix is ​​called an oscillator network model with separable coupling terms. Separable coupling terms are a necessary condition for using this invention. Specifically, the oscillator model in the oscillator network, i.e., the node dynamics model, is not limited to the linearized Stuart-Landau model in this embodiment, but is applicable to any oscillator network model that satisfies the condition of separable coupling terms.

[0063] Step S2, calculate the empirical covariance matrix using the covariance matrix C of the brain region BOLD signal. v As an empirical covariance matrix;

[0064] Covariance matrix C v elements in The calculation method is as follows:

[0065]

[0066] Where T represents the total duration of the BOLD signal, B i (t) represents the BOLD signal value of brain region i at time t. This represents the average BOLD signal in brain region i.

[0067] Step S3: Construct the target covariance matrix. Based on the prior knowledge of the dynamic model, convert the empirical covariance matrix into the target covariance matrix of the current linear system.

[0068] For the linearized Hopf whole-brain model used in this embodiment, the intrinsic angular frequencies at each node are uniform, i.e., ω i When ω, i = 1, ..., N, assuming that the variance of the signal at each node of the linear system is proportional to the variance of the empirical BOLD signal, the system objective covariance matrix has the following form: Where m is the scaling factor.

[0069] Step S4: Construct a dynamic parameter evaluation model to assess the rationality of the dynamic parameters. The dynamic parameter evaluation model includes a dynamic matrix construction module, a coupling matrix solving module, and a loss function calculation module, which are specifically defined as follows:

[0070] Dynamics Matrix Construction Module: This module receives a set of dynamic parameters to be evaluated, inputs these parameters into the dynamic model constructed in step S1, and thus constructs an uncoupled system dynamics matrix corresponding to the dynamic parameters. The dynamic parameters depend on the dynamic model used. In this embodiment, the form of the uncoupled system dynamics matrix is ​​as shown in A0 in step S1, and the dynamic parameters to be evaluated are a = [a1, ..., a2]. N ] represents the bifurcation coefficient of each node.

[0071] Coupled Matrix Solving Module: This module constructs a Sylvester equation with the coupling matrix as the variable using the system dynamics matrix (i.e., the system Jacobian matrix A), the target covariance matrix, and the dynamic parameters to be evaluated. It then uses the minimum residual approximation to obtain the unique solution to the Sylvester equation and converts it into a dynamic coupling matrix. The specific process is as follows:

[0072] From the Lyapunov equations for linear systems, we can obtain the Jacobian matrix A and the covariance matrix C of the system under steady state. vs Satisfying the following relationship: AC vs +C vs A T=-Q, where Q = σ 2 I is the covariance matrix of Gaussian white noise, σ 2 Let I be the variance of the Gaussian white noise, and I be the identity matrix. Separating the coupling terms in the Jacobian matrix yields:

[0073]

[0074] Since L is a symmetric matrix, L T =L = diag(S) - W, let K = -L = W - diag(S), and simplifying, we get:

[0075]

[0076] This equation is the Sylvester equation DK+KE=-Γ in D=E=C. vs In a special form, the Sylvester equation has a unique solution when D and -E do not have the same eigenvalues.

[0077] The expected system satisfies the following when it has the optimal coupling matrix: Based on the positive semidefiniteness of the covariance matrix, Since the eigenvalues ​​are non-negative, and because there is no strictly linear relationship between the nodes of the system, Eigenvalues ​​are strictly positive, therefore, as Figure 2 As shown, and Since the eigenvalues ​​are opposites of each other, it is impossible for two eigenvalues ​​to be identical. Therefore, the above special form of the Sylvester equation has a unique solution. Thus, it can be solved according to... Solving for K reveals that for any combination of dynamic parameters that makes the system stable, there exists a unique optimal coupling matrix W, where W is the off-diagonal part of the unique solution K.

[0078] Using an approximate solver (such as the solve_sylvester function in Python's Scipy.linalg computation library), the approximate solution K of the Sylvester equation above can be quickly obtained. Taking its off-diagonal part yields the optimal coupling matrix of the DTB model, i.e., the dynamic coupling matrix.

[0079] Loss function calculation module: This module adds the obtained dynamic coupling matrix into the DTB model, calculates the simulation covariance matrix and FC matrix, and uses one or more statistical indicators to construct a loss function to evaluate the simulation effect of parameter combinations and their corresponding coupling matrices.

[0080] In this embodiment, the system covariance matrix is ​​obtained by solving the Lyapunov equation of the aforementioned linear system. Each element of the system covariance matrix (off-diagonal elements represent the covariance between two nodes, and diagonal elements represent the variance of that node) is divided by the product of the standard deviations of the corresponding two nodes to obtain the simulated FC matrix FC based on the Pearson correlation coefficient. sim By calculating FC sim The empirical FC matrix FC calculated based on measured BOLD signals. emp The KS distance (ksdist) and Pearson correlation coefficient (corr) are used to define the joint loss function as: Loss = ksdist + (1 - corr), where the KS distance is the absolute value of the maximum difference between the cumulative distribution functions of the empirical and simulated FC matrices, and the Pearson correlation coefficient is defined as... in, Representing the empirical FC matrix FC emp The mean, To simulate the FC matrix FC sim The mean.

[0081] Step S5: The dynamic parameters and coupling matrix are simultaneously optimized using an optimization algorithm. An optimization range for the dynamic parameters is set. Using the loss function obtained from the dynamic parameter evaluation model as the standard, the optimization algorithm quickly explores the parameter space to find the optimal dynamic parameters and their corresponding optimal dynamic coupling matrix, thus constructing the optimized DTB model. In this embodiment, using group average data from a set of healthy individuals, the Success-History based Adaptive Differential Evolution with Linear Population Size Reduction (L-SHADE) algorithm was used to search for the optimal combination of dynamic parameters and its unique corresponding optimal dynamic coupling matrix within a reasonable range of the parameter space.

[0082] The specific implementation method of step S5 varies depending on the optimization algorithm used. The optimization algorithm that can be used is not limited to the specific machine learning algorithm described in this embodiment. For those skilled in the art, the optimization algorithm used in this step can have various modifications and variations.

[0083] The principle of the L-SHADE algorithm used in this embodiment is as follows:

[0084] Step S51, Initialization: Randomly generate a series of dynamic parameter combinations as the initial population, and initialize historical memory and external archive;

[0085] Step S52, Mutation: For each individual in the population, perform random mutation based on the current population individuals and the archived individuals (eliminated individuals) to generate mutated individuals;

[0086] Step S53, crossover: randomly swap the values ​​of certain dimensions in the parameter combinations of the mutant individual and the original individual to generate the experimental individual, that is, randomly select a position in the two parameter vectors and swap the value of that position in the two vectors.

[0087] Step S54: Select and evaluate whether the experimental individual has a lower loss function value than the corresponding original individual. If so, eliminate the original individual and put it into the external archive; otherwise, retain the original individual, which is the "survival of the fittest" principle.

[0088] Step S55: Update the successful historical memory and save the evolutionary algorithm coefficients of the successfully eliminated original individuals;

[0089] Step S56: Linear population size reduction, eliminating the number of inferior individuals with the worst loss function that exceed the new population size upper limit;

[0090] Step S57: Iterate, repeating steps S52 to S56 until the optimal loss function of the entire population converges or the maximum number of iterations is reached, then output the optimal combination of parameters for the individual.

[0091] The final optimization result and modeling effect of this embodiment are as follows: Figure 3 As shown, where, Figure 3 In the diagram, (a) is a histogram of the optimal dynamic parameters (i.e., the optimal values ​​of the bifurcation parameters in this embodiment). Figure 3 In this embodiment, (b) represents the optimal coupling matrix. Figure 3 In this embodiment, (c) represents the empirical FC matrix FC. emp With analog FC matrix FC sim Histogram, Figure 3 In the figure, (d) is the simulated FC matrix obtained in this embodiment. The above labels indicate the similarity between the simulated FC matrix and the empirical FC matrix (represented by the Pearson correlation coefficient, with a similarity of up to 98.7%) and the difference in the cumulative distribution function (represented by the KS distance, with a difference of no more than 0.5%). It can be seen that the DTB model modeled by this invention has better performance and its simulation effect is greatly improved.

[0092] Comparison of time consumption between the embodiments of the present invention and traditional modeling methods to achieve the same effect: Figure 4As shown, in terms of modeling speed, based on the average whole-brain BOLD signal of a single individual / group using 200 brain region maps, this method predicts the corresponding EC in 0.5 seconds for each group of dynamic parameters, while the traditional iterative optimization method takes 18.0 seconds to estimate the EC with the same effect. This method can achieve a speedup of 32 times compared to traditional modeling and optimization methods.

[0093] In summary, this method leverages the analytical solvability of linear models to provide a rapid and high-quality DTB modeling and optimization approach. Compared to traditional modeling and optimization methods, it significantly reduces the time required for modeling, optimizing, and analyzing individualized brain imaging data from a large number of subjects, greatly improves the fitting effect of DTB models, enhances model interpretability and result reliability, and promotes the implementation of personalized medical methods based on whole-brain dynamics models, such as individualized neuroimaging data analysis, prediction of neuromodulation efficacy, and quantitative diagnosis and subtype classification of neurological diseases.

[0094] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A rapid modeling method for digital twin brain reduction models used for brain imaging data integration, characterized in that, The method includes: Step S1: Construct an initial DTB model. Use linear dynamic oscillators to construct an oscillator network model with separable coupling terms as the linear system to be optimized. Step S2: Calculate the empirical covariance matrix, using the covariance matrix of the BOLD signal in the brain region as the empirical covariance matrix; Step S3: Construct the target covariance matrix. Based on the prior knowledge of the dynamic model, convert the empirical covariance matrix into the target covariance matrix of the current linear system. Step S4: Construct a dynamic parameter evaluation model to assess the rationality of the dynamic parameters. The dynamic parameter evaluation model includes a dynamic matrix construction module, a coupling matrix solving module, and a loss function calculation module. The dynamic matrix construction module receives a set of dynamic parameters to be evaluated and constructs an uncoupled system dynamic matrix based on the dynamic model. The coupling matrix solving module uses the system dynamic matrix, the target covariance matrix, and the dynamic parameters to be evaluated to construct a Sylvester equation with the coupling matrix as the variable. It approximates the unique solution of the Sylvester equation by minimizing the residuals and converts it into a dynamic coupling matrix. The loss function calculation module adds the obtained coupling matrix to the DTB model, calculates the system covariance matrix and the functional connectivity matrix, and constructs a loss function using one or more statistical indicators to evaluate the simulation effect of the parameter combination and its corresponding coupling matrix. Step S5: Use optimization algorithms to simultaneously optimize dynamic parameters and coupling matrices, set the optimization range of dynamic parameters, use the loss function obtained from the above dynamic parameter evaluation model as the standard, use optimization algorithms to quickly explore the parameter space, find the optimal dynamic parameters and their corresponding optimal dynamic coupling matrices, and construct the optimized DTB model. In step S1, a linearized Stuart-Landau oscillator model is used as the nodal dynamics model, and its expression is as follows: Among them, u=[u1,…,u N ] T Let A represent the model state vector and its node components, where N is the number of nodes in the DTB model, i.e., the number of oscillators in the oscillator network, δu is the oscillation of the model's state vector near the fixed point, and A is the Jacobian matrix of the 2N×2N system at the fixed point. jk Let a be the element in the j-th row and k-th column of the Jacobian matrix A. j W represents the bifurcation parameter of node j. jk Let η be the coupling weight between node j and node k in the coupling matrix W, where j = 1, ..., N, k = 1, ..., N, and η = [η1, ..., η]. 2N ] represents the vector representation of Gaussian white noise at each node, ω j Let be the eigenangular frequency of node j; Matrix A can be represented in the following block form: Among them, A xx A xy A yx A yy All are N×N matrices and have the following relationship: A xx =A yy =diag(aS)+W, A yx =-A xy =diag(ω), Where diag(v) represents a diagonal matrix with vector v as its main diagonal, and a = [a1, ..., a2]. N ] represents the vector representation of the bifurcation coefficients of each node. Let ω be the degree of each node in the coupling matrix W, where ω = [ω1, ..., ω]. N [ ] represents the vector representation of the intrinsic angular frequencies of each node; Further, let A be represented as: The two diagonal blocks of L are the graphical Laplace matrices of the coupling matrix W. Thus, the system Jacobian matrix A has been decomposed into the uncoupled system dynamic matrix A0 and the matrix L containing coupling information.

2. The method according to claim 1, characterized in that, In step S2, the covariance matrix C v elements in The calculation method is as follows: Where T represents the total duration of the BOLD signal, B i (t) represents the BOLD signal value of brain region i at time t. This represents the average BOLD signal in brain region i.

3. The method according to claim 2, characterized in that, In step S3, the variance of the signal at each node of the linear system is set to be proportional to the variance of the empirical BOLD signal. Then, the target covariance matrix of the system is expressed as: Where m is the scaling factor.

4. The method according to claim 1, characterized in that, The specific function of the dynamic matrix construction module in step S4 is to receive a set of dynamic parameters to be evaluated, input the dynamic parameters into the dynamic model constructed in step S1, and thus construct the uncoupled system dynamic matrix corresponding to the dynamic parameters.

5. The method according to claim 3, characterized in that, The specific function of the coupling matrix solving module in step S4 is as follows: The Jacobian matrix A and the covariance matrix C of the system under steady state are obtained from the Lyapunov equations of the linear system. vs The following relationship must be satisfied: AC vs +C vs A T =-Q, Where Q = σ 2 I is the covariance matrix of Gaussian white noise, σ 2 Let I be the variance of the Gaussian white noise, and I be the identity matrix. After separating the coupling terms in the Jacobian matrix, we obtain: Since L is a symmetric matrix, then L T =L = diag(S) - W, let K = -L = W - diag(S), then we get: This equation is the Sylvester equation DK+KE=-Γ in D=E=C. vs In a special form, the Sylvester equation has a unique solution when D and -E do not have the same eigenvalues. The system is configured to satisfy the following condition when it has the optimal coupling matrix: according to Solve for K, and from this we obtain the optimal coupling matrix W, which is the off-diagonal part of the unique solution K.

6. The method according to claim 5, characterized in that, The loss function in step S4 is calculated as follows: The system covariance matrix is ​​obtained by solving the Lyapunov equations of the linear system. Dividing each element of the system covariance matrix by the product of the standard deviations of the corresponding two nodes yields the simulated FC matrix FC, defined based on the Pearson correlation coefficient. sim By calculating FC sim The empirical FC matrix FC calculated based on measured BOLD signals. emp Using the KS distance (ksdist) and the Pearson correlation coefficient (corr), the joint loss function is obtained as: Loss = ksdist + (1 - corr), where the KS distance is the absolute value of the maximum difference between the empirical and simulated cumulative distribution functions of the FC matrix, and the Pearson correlation coefficient is defined as... in, Representing the empirical FC matrix FC emp The mean, To simulate the FC matrix FC sim The mean.

Images