Whole brain kinetic parameter fitting method based on Landau-Stewart oscillator network and application

Abstract

The invention discloses a whole brain kinetic parameter fitting method based on a Landau-Stewart oscillator network, which mainly comprises the following steps of: generating time sequence data for simulating a blood oxygen level dependent BOLD signal in real brain imaging data as simulation data; extracting data characteristics of the simulation data, and based on the data characteristics, selecting an initial value of a global coupling parameter # imgabs0 # and an initial value of a bifurcation parameter # imgabs1 #; implementing a gradient descent process and gradient correction based on data features; and fitting of a global coupling parameter # imgabs2 # and a bifurcation parameter # imgabs3 # is completed. According to the method, individual fitting can be carried out on the BOLD signal time sequence of each subject, the fitting performance of the model can be remarkably improved, and the method is used for evaluating and enhancing the applicability and accuracy of the model.

Description

Technical Field

[0001] The present invention belongs to the field of brain nerve technology, and specifically relates to a whole-brain dynamics parameter fitting method applicable to a Landau-Stewart oscillator network and its application. Background Art

[0002] For the time series data corresponding to the blood oxygen level dependent (BOLD) signals of a group of samples of test patients extracted from functional magnetic resonance imaging (fMRI) image data, the traditional method usually combines or averages the data features corresponding to the BOLD signal time series of the entire sample set, and then uses these integrated data to fit a whole-brain dynamics model applicable to the entire sample group, including but not limited to the Landau-Stewart oscillator network.

[0003] However, those skilled in the art have found that the fitting methods in the prior art have significant deficiencies in robustness and accuracy, and the fitting methods in the prior art are insufficient in capturing data features and cannot accurately reconstruct the dynamic characteristics of the BOLD signal time series.

[0004] The above information disclosed in this Background section is only for enhancement of understanding of the background of the invention and therefore it may contain information that does not form the prior art that is already known to a person of ordinary skill in the art. Summary of the invention

[0005] In order to overcome the above problems of the prior art, the present invention discloses a whole-brain dynamic parameter fitting method based on a Landau-Stewart oscillator network, comprising the following steps: S100, generating time series data of blood oxygen level-dependent BOLD signals in simulated real brain imaging data as simulation data, wherein the simulation data is used for global coupling parameters of the Landau-Stewart oscillator network and bifurcation parameters The fitting of S200, extracting data features of simulation data; S300, select global coupling parameters based on the data characteristics of the simulation data The initial value of S400, select bifurcation parameters based on the data characteristics of the simulation data The initial value of S500, calculating a gradient descent index of the simulation data according to the data features of the simulation data; S600, calculating an adaptive learning rate according to data features of the simulation data; S700, for the data features of the simulation data, introducing a data feature correction value for gradient correction of the data features; wherein, Introduced during gradient descent The data feature correction value is used as a correction measure for data feature fitting to accurately match the data features: ; in, The subscript t represents the bifurcation parameter Iterate the tth step. The superscript sim indicates the simulation value, and the superscript emp indicates the empirical value. Represents the bifurcation parameter The data feature correction value at the t-th iteration step, Represents the bifurcation parameter At the t-th iteration, the average value of the empirical FCD is Represents the bifurcation parameter At the t-th iteration, the average value of the simulated FCD is calculated; S800, construct an approximate loss function; S900, implement the gradient descent process, and use the data feature correction value And approximate loss function, perform gradient correction of data features during gradient descent, and find parameters and to complete the parameters of the Landau-Stewart oscillator network and The fitting of

[0006] Preferably, in step S100, The simulation data is generated using a Landau-Stewart oscillator network and randomly generated connection structures and parameter values ​​of various simulated brain networks.

[0007] Preferably, the method further comprises the following steps: S1000: Further use various simulated brain networks and random parameters to conduct a series of tests to evaluate the universality and accuracy of the fitting; when the universality and accuracy need to be improved, repeat S100 to S900 until the universality and accuracy meet the requirements.

[0008] Preferably, step S100 includes: (1) Setting up the simulated brain network: The connection density of the simulated brain network is set according to the connection density of the real brain network; the connection weight is uniformly distributed The simulated brain network is randomly selected from the network and is an undirected network.

[0009] Preferably, When the simulated brain network is a symmetric undirected network with M nodes, the SC matrix is ​​an M×M symmetric matrix, in which each row and column represents a node in the network; each element SC in the SC matrix ij represents the connection weight from node i to node j.

[0010] Preferably, In order to ensure the diversity and unpredictability of the simulated brain network, two network models were selected during the testing process: a connection structure based on the Watts-Strogatz (WS) network model, and a connection structure based on the Barabási-Albert (BA) network model; The fitting parameters were confirmed by comparing the connection structures of the two network models with the connection structures based on the real brain network. and effectiveness; among them, For the Watts-Strogatz (WS) network model, to change the reconnection probability That is, the probability of cutting off and randomly reconnecting the edges of the network is introduced, and additional random elements are introduced to experiment with multiple probability values; For the Barabási-Albert (BA) network model, scale-free networks with different connection densities are constructed, in which the number of connections between nodes is unevenly distributed: some nodes have a large number of connections, while other nodes have fewer.

[0011] Preferably, step S100 further includes: (2) Generate simulated data: A bifurcation parameter is randomly generated for each node in the simulated brain network As the benchmark bifurcation parameter value for evaluating the goodness of fit, the parameter set As a verification benchmark, it is input into the Landau-Stewart oscillator network, and the number of generated nodes is , the time length is The simulated data size is ; The generated simulation data is subjected to feature fitting and parameter estimation, wherein: the accuracy and robustness of the adopted parameter estimation method are evaluated by comparing the fitting results with the original parameters used to generate the simulation data, and the method is continuously improved; if the fitted parameters are highly consistent with the parameters used to generate the simulation data, the simulation data is trusted, wherein the trusted simulation data is not only used for the parameters of the Landau-Stewart oscillator network and It is suitable for clinical and other scenarios.

[0012] Preferably, step S200 specifically includes: The following key data features are determined based on the simulation data, and the parameters of the Landau-Stewart oscillator network are adjusted based on these key data features: and The fitting is performed, and these multiple key data features are also used as quantitative measures of the goodness of data fitting; The multiple key data features are determined in the following manner: (1) Determine functional connectivity (FC): The Pearson correlation coefficient between the simulated data of each pair of nodes in the simulated brain network is calculated to form a symmetric matrix, in which, Each element in the matrix is ​​the correlation coefficient r between the node pair xy , calculated as follows:

[0013] in, and Is a node and nodes in The value of the simulation data at time, and The nodes are and nodes The mean of the simulated data; is the time length of the simulated data; Using the correlation coefficient r xy represents the size of functional connectivity (FC); (2) Determine the functional connectivity matrix (dFC) at each moment: Using the Hilbert Transform Compute Node , p respectively at time The phase of the simulated data and , and calculate the phase difference between the node pair of node n and p at each time point, and then use the cosine function to calculate the phase coherence matrix Elements : ; in, The matrix size is , is the number of nodes, is the time length of the simulated data; (3) Determine dynamic functional connectivity (FCD): To compare the empirical and simulated FCDs, the Kolmogorov-Smirnov (KS) distance, which is used to quantify the maximum difference between the cumulative distribution functions of two samples, is used to measure the difference between the elements of the upper triangular part of the empirical FCD matrix calculated from the data and the simulated FCD matrix, where The range of Kolmogorov-Smirnov or KS distance is , the smaller the KS distance is, the better the fitting degree of the data features is. The empirical FCD matrix is ​​a matrix based on phase consistency; wherein, The empirical FCD matrix is ​​a A symmetric matrix, where is the length of the sliding window time window; The elements of the empirical FCD matrix is the bifurcation parameter Iteration Step and The functional connection matrix corresponding to the step and The upper triangle part and The cosine similarity between is expressed as follows: ; In this way, the empirical FCD matrix is ​​calculated from the empirical data; wherein the empirical data is obtained by generating simulation data in S100; (4) Calculate the probabilistic metastable substate (PMS): Step a. Feature vector extraction: From the bifurcation parameters The first feature vector ( ); The steps to extract the first eigenvector of the matrix are as follows: 1. Definition of matrix: Suppose there is a matrix , this is a The matrix represents samples and Features 2. Calculate the covariance matrix: First, calculate the matrix The covariance matrix of ;if Each column represents a feature and the row represents a sample, then the covariance matrix Calculated by: ,in yes The transposed matrix of 3. Calculate eigenvalues ​​and eigenvectors: Solve the covariance matrix The eigenvalues ​​and eigenvectors of satisfy the following equations: ,in is the eigenvalue, is the corresponding eigenvector; 4. Extract the first eigenvector: Sort the eigenvalues ​​from large to small. The eigenvector corresponding to the largest eigenvalue is the first eigenvector. ,in, Represents the direction with the largest directional variance in the data set; 5. Standardization: The feature vector Normalize it so that its modulus is 1, that is, its length is 1; Step b. Cluster analysis: Performing cluster analysis using the feature vectors extracted in step a to discover common and recurring connection patterns between different time points and different samples; wherein the cluster analysis uses a k-means clustering algorithm to cluster data points, i.e., dFC matrices at different time points, into multiple clusters based on the similarity of the feature vectors; Step c. Define the metastable substate: Considering that each cluster in the clustering result represents a unique metastable substate, the centroid of each cluster, i.e., the cluster center, is defined as the metastable substate; Step d. Probability and transition probability calculation: For each metastable substate, calculate its probability of appearing at all time points, and analyze the transition probability between different substates, that is, the probability of converting from one substate to another; Step e. Measure goodness of fit: The KL distance is used to compare the empirical data and the fitted data to obtain the probability of the metastable state, where Metastable substate The empirical probability of is the simulation probability on the same metastable substate; The KL distance between the empirical data and the fitted data is calculated to evaluate the fit of the model to the actually observed brain state probability distribution: .

[0014] Preferably, step S500 specifically includes: The frequency spectrum characteristics of each node extracted from the simulation data are selected as key indicators to characterize the numerical value of the bifurcation parameters of each node; wherein, This indicator is used to formulate the basis and evaluation criteria for parameter updating during the gradient descent process, so as to effectively guide the adjustment of bifurcation parameters during the gradient descent optimization process; This indicator is the node heterogeneity and is expressed as , where each element f of the indicator vector i is the spectrum feature of the i-th brain region; The process of calculating this indicator using simulated data includes: (1) Demean: The simulated data of each node is de-meaned so that the mean of the processed simulated data at all time points on each node is zero, which helps to remove the fixed deviation in the data and better analyze its fluctuations and changes; (2) Filtering processing: Filter the simulated data after removing the mean value of each node; (3) Detrend: Perform de-linear trend processing on the simulated data after de-averaging and filtering to eliminate possible linear changes in the data and ensure the accuracy of the analysis; (4) Fast Fourier Transform (FFT): Perform fast Fourier transform on the simulated data after de-linearization to convert it from the time domain to the frequency domain, thereby obtaining the spectrum characteristics of each frequency component; (5) Calculate the proportion of low-frequency spectrum features of each node in the simulated data: A low-frequency frequency band range where the blood oxygen level dependent (BOLD) signal time series is located is selected, and the value of the spectral characteristics that can characterize the simulated data within this range is calculated. The proportion of low-frequency spectral characteristics of each node of the simulated data is calculated by dividing the sum of the spectral characteristics of this frequency band by the sum of the spectral characteristics of the full frequency band range; Among them, the low-frequency frequency band range is generally concentrated between 0.01Hz and 0.1Hz. This frequency range is considered to be most relevant to neural activity in a resting state, mainly because the changes in blood oxygen caused by neural activity are relatively slow in time, and there is a large amount of physiological noise in the high-frequency band. Therefore, the BOLD signal in fMRI analysis mainly focuses on low-frequency fluctuations to effectively reflect the functional state of the brain and its network activity; Among them, the proportion f of any spectrum feature is calculated by the following formula: , Among them, in digital signal processing, the full frequency range is half of the sampling frequency; according to the Nyquist theorem, in order to unambiguously capture and reconstruct all frequency components in the signal, the sampling must be performed at least at twice the rate of the highest frequency component in the signal; when TR=2s, the full frequency range is defined as 0.01-0.25Hz; Among them, TR stands for repetition time, which refers to the interval between two adjacent radio frequency pulses during magnetic resonance imaging; specifically, TR is the time interval from the start of one imaging signal acquisition to the start of the next imaging signal acquisition, usually expressed in seconds or milliseconds; (6) Calculate the proportion of spectrum features in different frequency bands for each node in the simulated data: In addition to the low-frequency band where the BOLD signal time series is located, an arbitrary frequency band is selected, and the ratio of its spectral characteristics to the spectral characteristics of the full frequency band is calculated as a way to characterize the node dynamics; wherein the ratio is also calculated according to the proportion f of any spectral feature described above; Furthermore, this ratio is processed by standardization, exponential or multiple scaling to obtain an index more suitable for specific analysis or application, so as to better characterize the size of each node bifurcation parameter; (7) Calculate the spectrum characteristics of the full frequency band of each node of the simulated data: The frequency center, bandwidth, spectrum flatness, spectrum roll-off, and spectrum entropy of the entire frequency band are calculated to characterize the size of the bifurcation parameters of each node.

[0015] Preferably, The low-frequency amplitude fraction (fALFF) is preferentially selected as the index for calculating the gradient in the gradient descent process to achieve parameter update and improve the goodness of fitting of data spectrum characteristics, so as to effectively guide the adjustment of bifurcation parameters in the gradient descent optimization process. fALFF focuses on the amplitude of the low-frequency signal, that is, the ratio of the amplitude of the signal in the low-frequency signal band where the BOLD signal time series is located to the amplitude of the full-band signal, which is the node heterogeneity and is expressed as .

[0016] Preferably, step S600 specifically includes: According to the calculated standard deviation, mean, and amplitude characteristics of the indicators, an individualized adaptive learning rate strategy is selected for gradient descent, and the learning rate is adjusted according to the characteristics of the specific empirical data of the sample to improve the efficiency and universality of model learning; Since the parameter range is proportional to the standard deviation of the spectral feature index f, the standard deviation is selected to adjust the iterative learning rate of gradient descent:

[0017] in, is the coefficient of the learning rate, which is selected according to the characteristics of the simulated data; for example, , The spectral feature index f is derived from the spectral features of each node extracted from the simulation data, which is used to characterize the numerical value of the bifurcation parameters of each node; the spectral feature index f is used as a gradient descent index to formulate the basis and evaluation criteria for parameter updating during the gradient descent process, so as to effectively guide the adjustment of the bifurcation parameters during the gradient descent optimization process.

[0018] Preferably, step S800 specifically includes: (1) Introducing the true loss function: In order to more accurately determine whether the iteration of gradient descent has reached the target value, a loss function is introduced as an evaluation criterion; Bifurcation parameter value is known and serves as the benchmark bifurcation parameter value for evaluating the goodness of fit. Each node randomly generates a bifurcation parameter , get the parameter set As a validation benchmark, it is used to construct a true loss function to verify the performance of parameter fitting in terms of quantitative differences and qualitative correlations; Based on this, the true loss function consists of two parts: a. Differences between fitting parameters and benchmark parameters: ; b. Similarity between the fitted parameters and the benchmark parameters: ; in, It is step The fitting parameter value when ; The actual loss function calculation formula is defined as follows: ; in, refers to the true loss function, is the bifurcation parameter The difference between the fitted parameters and the true parameters at the t-th iteration step, is the bifurcation parameter The similarity between the fitting parameters and the benchmark parameters at the t-th iteration step, std means the standard deviation, and mean means the mean; (2) Construct an approximate loss function: An approximate loss function consisting of the goodness of fit of data features, the convergence degree of spectral features and the goodness of fit of parameters is constructed to evaluate the degree of fit between parameters and data features, thereby obtaining the parameter combination with the minimum loss function value: ; in, refers to the approximate loss function, is the bifurcation parameter The difference in data characteristics between the fitted data and the empirical data at the t-th iteration step, is the bifurcation parameter The difference between the spectral characteristic index of the fitting data at the t-th iteration step and the spectral characteristic index of the empirical data, yes The correlation between the moment fitting parameters and the spectral characteristic indicators of the empirical data; a. Goodness of fit of data characteristics: According to the fitting requirements of data characteristics, the corresponding indicators that can quantify the goodness of fit are selected, including: the correlation between the simulated and benchmark FC matrices, the KS distance between FCD matrices, and the KL distance between PMS vectors; Preferably, for the dynamic functional connectivity characteristics of the data, the KS distance used to quantify the maximum difference between the cumulative distributions of the empirical and simulated FCD matrices is used to measure the difference between the elements of the upper triangular part of the empirical matrix and the simulated FCD matrix: , represents the data feature distance between the empirical data and the fitted data at step t, Represents the dynamic functional connectivity features of empirical data, represents the dynamic functional connectivity features of the fitted data at step t; b. Convergence degree of spectrum characteristics: According to the fitting parameter values ​​G and Calculate the corresponding fitted data value; then, calculate At each node, the empirical data The maximum absolute value of the difference is used as a quantitative measure of the goodness of fit of the spectral feature. : , in, For parameters At the t-th iteration, the spectral feature index extracted from the fitting data is: is the spectral characteristic index of the empirical data; max means taking the maximum value; c. Parameter fit goodness: In the absence of known actual parameter values, qualitative fitting performance is evaluated using the correlation between the fitted parameters and an empirical metric as part of the loss function: ; Based on the true loss function of the benchmark bifurcation parameter and the exploration of ways to quantify the goodness of fit from various angles, an approximate loss function of the gradient descent process is constructed to evaluate the degree of fit between the parameters and the data features, thereby obtaining the parameter combination with the minimum loss function value to ensure that the model exhibits good fitting ability and wide applicability in different application scenarios:

[0019] in:

[0020]

[0021]

[0022] d. Other standardization methods: When calculating the approximate loss function, different standardization methods can be used to optimize the performance of the model by adjusting the mean and standard deviation of the data features to have a uniform scale and distribution, including: i. Max-min normalization: scaling all features to a specific range (usually 0 to 1):

[0023] in, is the maximum value operation, is the minimum operation; ii. Mean standardization: Divide by the mean value of the feature to change the scale of the data, making the data scale uniform and easier for algorithm processing:

[0024] in, is the mean operation; iii. Robust scaling: Use the median and interquartile range (IQR) of the data to scale the data, which is more robust to outliers:

[0025] in, is the median operation; iv. Mean normalization: .

[0026] Preferably, step S900 specifically includes: During the gradient descent process, the parameters of each generation are used Generate new simulated data and update the tuning parameters by calculating the spectral characteristics , to better match the benchmark bifurcation parameter values, including: (1) Initial parameter setting: The optimal parameter combination obtained from the initial value selection in the gradient descent process and Start to proceed; (2) Substitute the dynamic equations of the brain model (such as the Watts-Strogatz (WS) network model or the Barabási-Albert (BA) network model) to generate simulation data: The parameters for the iteration Substitute the dynamic equations and solve them through numerical integration (such as Euler method, Runge-Kutta method) to obtain the time evolution of the system under given parameters and generate simulation data corresponding to each node. ;in , ; (3) Calculate the spectral characteristics of the simulated data: The generated simulated data Perform spectrum analysis and calculate ; The process includes performing a fast Fourier transform (FFT) on the analog data and then extracting the relevant spectral features fALFF (it should be noted that frequency center, bandwidth, spectral peak, spectral flatness, spectral roll-off and spectral entropy can also be used); (4) Calculate the spectral feature difference: calculate and The difference between them is used as the updated gradient of the parameters: ; (5) Calculate the data feature difference: ; (6) Parameter update: Using the difference in spectral features To update the parameters , so that it gradually approaches the value that can better fit the real data; Among them, the parameter update rule is based on a simple gradient descent method, and the updated parameters are expressed as:

[0027]

[0028] in, is the learning rate, which controls the step size of parameter updates and is calculated by the standard deviation of the indicator; (7) Iteration process: Repeat the above steps and calculate the value of the approximate loss function at each iteration; (8) Select the parameter value at the minimum loss function value: The bifurcation parameter value at the position where the approximate loss function value is the smallest is the value closest to the benchmark bifurcation parameter, that is, the better fitting parameter obtained.

[0029] In addition, the present invention also discloses a method for extracting blood oxygen level dependent BOLD signal time series data, comprising the following steps: S1, data preprocessing; Wherein, the data preprocessing is used to correct various potential influencing factors in fMRI data to make the data suitable for further analysis; S2, spatial standardization; Wherein, the spatial normalization is used to convert fMRI data into standard spatial coordinates to ensure that signals of different samples can be compared; S3, define template; S4, extract BOLD time series; S5. Remove confounding factors; It is used to remove possible confounding factors from the extracted time series.

[0030] Preferably, in step S1, the data preprocessing includes: S11. Timing correction: Correct the timing of different slices collected at different times during the fMRI scan; S12, motion correction: correcting or compensating for the signal deviation caused by the head movement; S13, spatial registration: aligning functional images with anatomical images and mapping different samples to the same standard template; S14, spatial smoothing: Improve the signal-to-noise ratio by smoothing the data.

[0031] Preferably, in step S2, Standard templates used include the Montreal Neurological Institute (MNI) template or the Talairach template to map sample images to the standard template.

[0032] Preferably, step S3 specifically includes: According to the research objectives of interest, a brain region template is defined or selected, and each brain region in the template corresponds to each node in the subsequent model; among them, the brain region segmentation templates used include AAL template, Harvard-Oxford template, and Dosenbach template.

[0033] Preferably, step S4 specifically includes: Use the brain region mask in the template to perform segmentation in the preprocessed and normalized fMRI images; The average BOLD signal time series of all voxels within each brain region defined in the template was calculated to obtain a time series representing the region, corresponding to the time series of each node in the subsequent model.

[0034] Preferably, Confounding factors include head movement, whole-brain average signal, CSF signal, etc. Regression analysis was used to remove confounding factors to ensure that the extracted time series more accurately reflected neural activity.

[0035] Preferably, Fix all bifurcation parameters to , a parameter scanning method is used to determine the optimal global coupling constant that minimizes the Kolmogorov-Smirnov (KS) distance of the FCD matrix fitting of the time series .

[0036] Preferably, in step S400, Fixed optimal Value, yes exist The parameter scan with a step size of 0.01 was performed in the range to find the best fit for the BOLD signal time series FCD. The initial value of .

[0037] Preferably, In order to overcome the possible loss of individual information caused by group averaging, a set of individualized parameters is estimated for each subject.

[0038] Preferably, Select the initial value suitable for each subject , ensuring the personalization and accuracy of BOLD signal time series processing.

[0039] Preferably, in step S500, A convergence evaluation criterion was developed based on the low-frequency amplitude fraction (fALFF) indicator extracted from the BOLD signal time series to guide the gradient descent optimization of the bifurcation parameters.

[0040] Preferably, in step S600, According to the standard deviation of the fALFF index, an individualized adaptive learning rate strategy is selected for gradient descent, and the learning rate is adjusted according to the specific BOLD signal time series characteristics of the subjects to improve the efficiency and universality of model learning.

[0041] Preferably, in step S800, During the gradient descent process, an approximate loss function is constructed based on the true loss function of the benchmark bifurcation parameters to evaluate the degree of fit between the parameters and the time series characteristics of the BOLD signal, thereby obtaining the parameter combination with the minimum loss function value and ensuring that the model exhibits good fitting ability and wide applicability in different application scenarios.

[0042] Preferably, in step S900, Select the parameter combination that minimizes the approximate loss function, and this set of parameters is the optimal fitting parameters.

[0043] Compared with the prior art, the present invention has the following advantages: The present invention can perform individualized fitting of the BOLD signal time series of each subject, significantly improve the fitting performance of the model, and be used to evaluate and enhance the applicability and accuracy of the model. The present invention can also better capture and retain the unique imaging data features of each sample, enhance the credibility of brain data, and broaden the application scenarios of brain data. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] By reading the detailed description of the preferred specific embodiments below, various other advantages and benefits of the present invention will become clear to those of ordinary skill in the art. The drawings in the specification are only for the purpose of illustrating the preferred embodiments and are not considered to be limitations of the present invention. Obviously, the drawings described below are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can also be obtained based on these drawings without creative work. Moreover, the same reference numerals are used to represent the same components throughout the drawings.

[0045] Figure 1 is a schematic diagram of using a generative model to verify a fitting method in one embodiment of the present invention, wherein: The left side shows the form of the simulated brain network being tested, and the right side shows the randomly selected bifurcation parameters as the benchmark bifurcation parameters. After substituting them into the model, simulated data are generated, and then the parameters of the simulated data are estimated and compared with the benchmark bifurcation parameters. Figure 2 It is a schematic diagram of the correlation between different spectral characteristics and parameters, where: The horizontal axis represents 8 spectrum characteristic indicators, which are: low-frequency power proportion, low-frequency amplitude proportion, half-band amplitude proportion, low-frequency amplitude sum, low-frequency frequency center proportion, low-frequency bandwidth proportion, low-frequency spectrum flatness proportion, low-frequency spectrum roll-off ratio and low-frequency spectrum entropy proportion; The ordinate represents the correlation between the bifurcation parameter value and the indicator value; Figure 3 FIG. 1 is a schematic diagram of selecting different frequency bands to calculate the gradient index fALFF during parameter update iteration in another embodiment of the present invention, wherein: Figure 3 (a) shows the goodness of fit of the data characteristics of the simulated data (wherein, in the left part of (a), the ordinate represents the KS distance), the difference of parameter fitting (wherein, in the middle part of (a), the ordinate represents the difference between the fitted parameters and the benchmark parameters) and the correlation of parameter fitting (wherein, in the right part of (a), the ordinate represents the correlation between the fitted parameters and the benchmark parameters) under the indicators of five frequency bands (0.01-0.07Hz, 0.01-0.10Hz, 0.01-0.15Hz, 0.01-0.20Hz, 0.01-0.25Hz) when TR of the simulated data is 1; Figure 3 (b) shows the goodness of fit of the data characteristics of the simulated data (wherein, in the left part of (b), the ordinate represents the KS distance), the difference of parameter fitting (wherein, in the middle part of (b), the ordinate represents the difference between the fitted parameters and the benchmark parameters) and the correlation of parameter fitting (wherein, in the right part of (b), the ordinate represents the correlation between the fitted parameters and the benchmark parameters) under the indicators of five frequency bands (0.01-0.05Hz, 0.01-0.15Hz, 0.01-0.25Hz, 0.01-0.35Hz, 0.01-0.45Hz) when TR of the simulated data is 2; Figure 4 FIG. 1 is a schematic diagram comparing the fitting performance of different initial value selection methods in another embodiment of the present invention, wherein: The left figure shows the error comparison of parameter fitting at each node under three initial value methods; The left ordinate of the upper right figure represents the correlation comparison of parameter fitting under the three initial value methods, and the right ordinate represents the error comparison of parameter fitting under the three initial value methods; The lower right figure shows the comparison of the goodness of fit of the data characteristics of the simulated data under the three initial value methods; Figure 5 is a schematic diagram of a process for calculating an approximate loss function in another embodiment of the present invention, wherein: Approximate loss function (see Figure 5 The lower part) combines the three elements of KS distance, indicator convergence and similarity with the actual loss function (see Figure 5 When compared with the upper part, a high degree of consistency was shown; Figure 6 It is a schematic diagram of the overall process, improvement effect and application prospect of the Landau-Stewart oscillator network in a simulation data embodiment of the present invention, wherein: The upper left part shows the structural connectivity, data morphology, data features (FCD), and gradient descent indicators of the extracted simulation data; The middle left part shows the initial value selection process for fitting simulation data; The lower left part shows the process of fitting the gradient descent of the simulated data and calculating the approximate loss function; The lower right side shows the performance of fitting the simulated data; The upper right part shows the application prospect of the present invention; Figure 7 is a schematic diagram of the improvement effect obtained by verifying the present invention in the simulation data verification, wherein: The left figure shows the difference in parameter fitting (the ordinate is the difference between the fitted parameters and the reference parameters (after standardization)); The middle figure shows the correlation of parameter fitting (the ordinate is the correlation between the fitted parameters and the benchmark parameters); The right figure shows the goodness of fit of the data features of the simulated data (the ordinate represents the KS distance); Figure 8 is a schematic diagram of fitting performance when the network scale is N=264 in another embodiment of the present invention; Fig. 9 It is a schematic diagram of parameter distribution and fitting of patients with depression (MDD) and a healthy control group (HC) in a real data set application embodiment of the present invention; Fig.10 is a schematic diagram showing the results of MDD parameter analysis in a real data set application embodiment of the present invention; Fig.11 It is a schematic diagram of the correlation between the parameters of MDD and the Hamilton Depression Rating Scale (HAMD) score in an embodiment of the present invention applied to a real data set. DETAILED DESCRIPTION

[0046] The following will refer to the attached Figures 1 to 11 Specific embodiments of the present invention are described in detail. Although specific embodiments of the present invention are shown in the accompanying drawings, it should be understood that the present invention can be implemented in various forms and should not be limited by the embodiments described herein. On the contrary, these embodiments are provided to enable a more thorough understanding of the present invention and to fully convey the scope of the present invention to those skilled in the art.

[0047] It should be noted that certain words are used in the specification and claims to refer to specific components. Those skilled in the art should understand that technicians may use different nouns to refer to the same component. This specification and claims do not use the difference in nouns as a way to distinguish components, but use the difference in the functions of the components as the criterion for distinction. As mentioned throughout the specification and claims, "including" or "comprising" is an open term, so it should be interpreted as "including but not limited to". The subsequent description of the specification is a preferred embodiment of the present invention, but the description is based on the general principles of the specification and is not intended to limit the scope of the present invention. The scope of protection of the present invention shall be determined by the attached claims.

[0048] To facilitate understanding of the embodiments of the present invention, further explanation will be given below by taking specific embodiments as examples in conjunction with the accompanying drawings, and each of the accompanying drawings does not constitute a limitation on the embodiments of the present invention.

[0049] In one embodiment, the present invention discloses a whole-brain dynamic parameter fitting method based on a Landau-Stewart oscillator network, comprising the following steps: S100, generating time series data of blood oxygen level-dependent BOLD signals in simulated real brain imaging data as simulation data, wherein the simulation data is used for global coupling parameters of the Landau-Stewart oscillator network and bifurcation parameters The fitting of S200, extracting data features of simulation data; S300, select global coupling parameters based on the data characteristics of the simulation data The initial value of S400, select bifurcation parameters based on the data characteristics of the simulation data The initial value of S500, calculating a gradient descent index of the simulation data according to the data features of the simulation data; S600, calculating an adaptive learning rate according to data features of the simulation data; S700, for the data features of the simulation data, introducing a data feature correction value for gradient correction of the data features; wherein, Introduced during gradient descent The data feature correction value is used as a correction measure for data feature fitting to accurately match the data features: ; in, The subscript t represents the bifurcation parameter Iterate the tth step. The superscript sim indicates the simulation value, and the superscript emp indicates the empirical value. Represents the bifurcation parameter The data feature correction value at the t-th iteration step, Represents the bifurcation parameter At the t-th iteration, the average value of the empirical FCD is Represents the bifurcation parameter At the t-th iteration, the average value of the simulated FCD is calculated; S800, construct an approximate loss function; S900, implement the gradient descent process, and use the data feature correction value And approximate loss function, perform gradient correction of data features during gradient descent, and find parameters and to complete the parameters of the Landau-Stewart oscillator network and The fitting of

[0050] Preferably, in step S100, The simulation data is generated using a Landau-Stewart oscillator network and randomly generated connection structures and parameter values ​​of various simulated brain networks.

[0051] Preferably, the method further comprises the following steps: S1000: Further use various simulated brain networks and random parameters to conduct a series of tests to evaluate the universality and accuracy of the fitting; when the universality and accuracy need to be improved, repeat S100 to S900 until the universality and accuracy meet the requirements.

[0052] Preferably, step S100 includes: (1) Setting up the simulated brain network: The connection density of the simulated brain network is set according to the connection density of the real brain network; the connection weight is uniformly distributed The simulated brain network is randomly selected from the network and is an undirected network.

[0053] Preferably, When the simulated brain network is a symmetric undirected network with M nodes, the SC matrix is ​​an M×M symmetric matrix, in which each row and column represents a node in the network; each element SC in the SC matrix ij represents the connection weight from node i to node j.

[0054] Preferably, In order to ensure the diversity and unpredictability of the simulated brain network, two network models were selected during the testing process: a connection structure based on the Watts-Strogatz (WS) network model, and a connection structure based on the Barabási-Albert (BA) network model; The fitting parameters were confirmed by comparing the connection structures of the two network models with the connection structures based on the real brain network. and effectiveness; among them, For the Watts-Strogatz (WS) network model, to change the reconnection probability That is, the probability of cutting off and randomly reconnecting the edges of the network is introduced, and additional random elements are introduced to experiment with multiple probability values; For the Barabási-Albert (BA) network model, scale-free networks with different connection densities are constructed, in which the number of connections between nodes is unevenly distributed: some nodes have a large number of connections, while other nodes have fewer.

[0055] Preferably, step S100 further includes: (2) Generate simulated data: A bifurcation parameter is randomly generated for each node in the simulated brain network As the benchmark bifurcation parameter value for evaluating the goodness of fit, the parameter set As a verification benchmark, it is input into the Landau-Stewart oscillator network, and the number of generated nodes is , the time length is The simulated data size is ; The generated simulation data is subjected to feature fitting and parameter estimation, wherein: the accuracy and robustness of the adopted parameter estimation method are evaluated by comparing the fitting results with the original parameters used to generate the simulation data, and the method is continuously improved; if the fitted parameters are highly consistent with the parameters used to generate the simulation data, the simulation data is trusted, wherein the trusted simulation data is not only used for the parameters of the Landau-Stewart oscillator network and It is suitable for clinical and other scenarios.

[0056] Preferably, step S200 specifically includes: The following key data features are determined based on the simulation data, and the parameters of the Landau-Stewart oscillator network are adjusted based on these key data features: and The fitting is performed, and these multiple key data features are also used as quantitative measures of the goodness of data fitting; The multiple key data features are determined in the following manner: (1) Determine functional connectivity (FC): The Pearson correlation coefficient between the simulated data of each pair of nodes in the simulated brain network is calculated to form a symmetric matrix, in which, Each element in the matrix is ​​the correlation coefficient r between the node pair xy , calculated as follows:

[0057] in, and Is a node and nodes in The value of the simulation data at time, and The nodes are and nodes The mean of the simulated data; is the time length of the simulated data; Using the correlation coefficient r xy represents the size of functional connectivity (FC); (2) Determine the functional connectivity matrix (dFC) at each moment: Using the Hilbert Transform Compute Node , p respectively at time The phase of the simulated data and , and calculate the phase difference between the node pair of node n and p at each time point, and then use the cosine function to calculate the phase coherence matrix Elements : ; in, The matrix size is , is the number of nodes, is the time length of the simulated data; (3) Determine dynamic functional connectivity (FCD): To compare the empirical and simulated FCDs, the Kolmogorov-Smirnov (KS) distance, which is used to quantify the maximum difference between the cumulative distribution functions of two samples, is used to measure the difference between the elements of the upper triangular part of the empirical FCD matrix calculated from the data and the simulated FCD matrix, where The range of Kolmogorov-Smirnov or KS distance is , the smaller the KS distance is, the better the fitting degree of the data features is. The empirical FCD matrix is ​​a matrix based on phase consistency; wherein, The empirical FCD matrix is ​​a A symmetric matrix, where is the length of the sliding window time window; The elements of the empirical FCD matrix is the bifurcation parameter Iteration Step and The functional connection matrix corresponding to the step and The upper triangle part and The cosine similarity between is expressed as follows: ; In this way, the empirical FCD matrix is ​​calculated from the empirical data; wherein the empirical data is obtained by generating simulation data in S100; (4) Calculate the probabilistic metastable substate (PMS): Step a. Feature vector extraction: From the bifurcation parameters The first feature vector ( ); The steps to extract the first eigenvector of the matrix are as follows: 1. Definition of matrix: Suppose there is a matrix , this is a The matrix represents samples and Features 2. Calculate the covariance matrix: First, calculate the matrix The covariance matrix of ;if Each column represents a feature and the row represents a sample, then the covariance matrix Calculated by: ,in yes The transposed matrix of 3. Calculate eigenvalues ​​and eigenvectors: Solve the covariance matrix The eigenvalues ​​and eigenvectors of satisfy the following equations: ,in is the eigenvalue, is the corresponding eigenvector; 4. Extract the first eigenvector: Sort the eigenvalues ​​from large to small. The eigenvector corresponding to the largest eigenvalue is the first eigenvector. ,in, Represents the direction with the largest directional variance in the data set; 5. Standardization: The feature vector Normalize it so that its modulus is 1, that is, its length is 1; Step b. Cluster analysis: Performing cluster analysis using the feature vectors extracted in step a to discover common and recurring connection patterns between different time points and different samples; wherein the cluster analysis uses a k-means clustering algorithm to cluster data points, i.e., dFC matrices at different time points, into multiple clusters based on the similarity of the feature vectors; Step c. Define the metastable substate: Considering that each cluster in the clustering result represents a unique metastable substate, the centroid of each cluster, i.e., the cluster center, is defined as the metastable substate; Step d. Probability and transition probability calculation: For each metastable substate, calculate its probability of appearing at all time points, and analyze the transition probability between different substates, that is, the probability of converting from one substate to another; Step e. Measure goodness of fit: The KL distance is used to compare the empirical data and the fitted data to obtain the probability of the metastable state, where Metastable substate The empirical probability of is the simulation probability on the same metastable substate; The KL distance between the empirical data and the fitted data is calculated to evaluate the fit of the model to the actually observed brain state probability distribution: .

[0058] Preferably, step S500 specifically includes: The frequency spectrum characteristics of each node extracted from the simulation data are selected as key indicators to characterize the numerical value of the bifurcation parameters of each node; wherein, This indicator is used to formulate the basis and evaluation criteria for parameter updating during the gradient descent process, so as to effectively guide the adjustment of bifurcation parameters during the gradient descent optimization process; This indicator is the node heterogeneity and is expressed as , where each element f of the indicator vector i is the spectrum feature of the i-th brain region; The process of calculating this indicator using simulated data includes: (1) Demean: The simulated data of each node is de-meaned so that the mean of the processed simulated data at all time points on each node is zero, which helps to remove the fixed deviation in the data and better analyze its fluctuations and changes; (2) Filtering processing: Filter the simulated data after removing the mean value of each node; (3) Detrend: Perform de-linear trend processing on the simulated data after de-averaging and filtering to eliminate possible linear changes in the data and ensure the accuracy of the analysis; (4) Fast Fourier Transform (FFT): Perform fast Fourier transform on the simulated data after de-linearization to convert it from the time domain to the frequency domain, thereby obtaining the spectrum characteristics of each frequency component; (5) Calculate the proportion of low-frequency spectrum features of each node in the simulated data: A low-frequency frequency band range where the blood oxygen level dependent (BOLD) signal time series is located is selected, and the value of the spectral characteristics that can characterize the simulated data within this range is calculated. The proportion of low-frequency spectral characteristics of each node of the simulated data is calculated by dividing the sum of the spectral characteristics of this frequency band by the sum of the spectral characteristics of the full frequency band range; Among them, the low-frequency frequency band range is generally concentrated between 0.01Hz and 0.1Hz. This frequency range is considered to be most relevant to neural activity in a resting state, mainly because the changes in blood oxygen caused by neural activity are relatively slow in time, and there is a large amount of physiological noise in the high-frequency band. Therefore, the BOLD signal in fMRI analysis mainly focuses on low-frequency fluctuations to effectively reflect the functional state of the brain and its network activity; Among them, the proportion f of any spectrum feature is calculated by the following formula: , Among them, in digital signal processing, the full frequency range is half of the sampling frequency; according to the Nyquist theorem, in order to unambiguously capture and reconstruct all frequency components in the signal, the sampling must be performed at least at twice the rate of the highest frequency component in the signal; when TR=2s, the full frequency range is defined as 0.01-0.25Hz; Among them, TR stands for repetition time, which refers to the interval between two adjacent radio frequency pulses during magnetic resonance imaging; specifically, TR is the time interval from the start of one imaging signal acquisition to the start of the next imaging signal acquisition, usually expressed in seconds or milliseconds; The spectrum feature index f is derived from the spectrum features of each node extracted from the simulation data, which is used to characterize the numerical value of the bifurcation parameters of each node; the spectrum feature index f is used as a gradient descent index to formulate the basis and evaluation criteria for parameter updating during the gradient descent process, so as to effectively guide the adjustment of the bifurcation parameters during the gradient descent optimization process; Exemplarily, the characteristic index or feature is the amplitude, and accordingly, a low frequency range (e.g. ), calculate the amplitude and within this range, and use this as the characteristic index and of the sub-band. Further, in the full frequency band, it is defined as In the example, the amplitude and within this range are calculated and used as the characteristic index and of the full frequency band.

[0059] (6) Calculate the proportion of spectrum features in different frequency bands for each node in the simulated data: In addition to the low-frequency band where the BOLD signal time series is located, an arbitrary frequency band is selected, and the ratio of its spectral characteristics to the spectral characteristics of the full frequency band is calculated as a way to characterize the node dynamics; wherein the ratio is also calculated according to the proportion f of any spectral feature described above; Furthermore, this ratio is processed by standardization, exponential or multiple scaling to obtain an index more suitable for specific analysis or application, so as to better characterize the size of each node bifurcation parameter; (7) Calculate the spectrum characteristics of the full frequency band of each node of the simulated data: The frequency center, bandwidth, spectrum flatness, spectrum roll-off, and spectrum entropy of the entire frequency band are calculated to characterize the size of the bifurcation parameters of each node.

[0060] Preferably, The low-frequency amplitude fraction (fALFF) is preferentially selected as the index for calculating the gradient in the gradient descent process to achieve parameter update and improve the goodness of fitting of data spectrum characteristics, so as to effectively guide the adjustment of bifurcation parameters in the gradient descent optimization process. fALFF focuses on the amplitude of the low-frequency signal, that is, the ratio of the amplitude of the signal in the low-frequency signal band where the BOLD signal time series is located to the amplitude of the full-band signal, which is the node heterogeneity and is expressed as .

[0061] Preferably, in step S300, Fix all bifurcation parameters to , a parameter sweep method was used to determine the optimal global coupling constant that minimizes the Kolmogorov-Smirnov (KS) distance of the FCD matrix fit of the data .

[0062] Preferably, in step S400, Fixed optimal Value, yes exist The parameter was scanned with a step size of 0.01 in the range to find the best fit to the data FCD. The initial value of .

[0063] Preferably, In order to overcome the possible loss of individual information caused by group averaging, a set of individualized parameters is estimated for each subject.

[0064] Preferably, Select the initial value suitable for each subject , ensuring the personalization and accuracy of data processing.

[0065] Preferably, A convergence evaluation criterion was developed based on the low-frequency amplitude fraction (fALFF) metric extracted from the data to guide the gradient descent optimization of the bifurcation parameters.

[0066] Preferably, in step S800, During the gradient descent process, an approximate loss function is constructed based on the true loss function of the benchmark bifurcation parameters to evaluate the degree of fit between the parameters and the data features, thereby obtaining the parameter combination with the minimum loss function value and ensuring that the model exhibits good fitting ability and wide applicability in different application scenarios.

[0067] Preferably, in step S900, Select the parameter combination that minimizes the approximate loss function, and this set of parameters is the optimal fitting parameters.

[0068] Preferably, step S600 specifically includes: According to the calculated standard deviation, mean, and amplitude characteristics of the indicators, an individualized adaptive learning rate strategy is selected for gradient descent, and the learning rate is adjusted according to the characteristics of the specific empirical data of the sample to improve the efficiency and universality of model learning; Since the parameter range is proportional to the standard deviation of the spectral feature index f, the standard deviation is selected to adjust the iterative learning rate of gradient descent:

[0069] in, is the coefficient of the learning rate, which is selected according to the characteristics of the simulated data; for example, , The spectral feature index f is derived from the spectral features of each node extracted from the simulation data, which is used to characterize the numerical value of the bifurcation parameters of each node; the spectral feature index f is used as a gradient descent index to formulate the basis and evaluation criteria for parameter updating during the gradient descent process, so as to effectively guide the adjustment of the bifurcation parameters during the gradient descent optimization process.

[0070] Preferably, step S800 specifically includes: (1) Introducing the true loss function: In order to more accurately determine whether the iteration of gradient descent has reached the target value, a loss function is introduced as an evaluation criterion; Bifurcation parameter value is known and serves as the benchmark bifurcation parameter value for evaluating the goodness of fit. Each node randomly generates a bifurcation parameter , get the parameter set As a validation benchmark, it is used to construct a true loss function to verify the performance of parameter fitting in terms of quantitative differences and qualitative correlations; Based on this, the true loss function consists of two parts: a. Differences between fitting parameters and benchmark parameters: ; b. Similarity between the fitted parameters and the benchmark parameters: ; in, It is step The fitting parameter value when ; The actual loss function calculation formula is defined as follows: ; in, refers to the true loss function, is the bifurcation parameter The difference between the fitted parameters and the true parameters at the t-th iteration step, is the bifurcation parameter The similarity between the fitting parameters and the benchmark parameters at the t-th iteration step, std means the standard deviation, and mean means the mean; (2) Construct an approximate loss function: An approximate loss function consisting of the goodness of fit of data features, the convergence degree of spectral features and the goodness of fit of parameters is constructed to evaluate the degree of fit between parameters and data features, thereby obtaining the parameter combination with the minimum loss function value: ; in, refers to the approximate loss function, is the bifurcation parameter The difference in data characteristics between the fitted data and the empirical data at the t-th iteration step, is the bifurcation parameter The difference between the spectral characteristic index of the fitting data at the t-th iteration step and the spectral characteristic index of the empirical data, yes The correlation between the moment fitting parameters and the spectral characteristic indicators of the empirical data; a. Goodness of fit of data characteristics: According to the fitting requirements of data characteristics, the corresponding indicators that can quantify the goodness of fit are selected, including: the correlation between the simulated and benchmark FC matrices, the KS distance between FCD matrices, and the KL distance between PMS vectors; Preferably, for the dynamic functional connectivity characteristics of the data, the KS distance used to quantify the maximum difference between the cumulative distributions of the empirical and simulated FCD matrices is used to measure the difference between the elements of the upper triangular part of the empirical matrix and the simulated FCD matrix: , represents the data feature distance between the empirical data and the fitted data at step t, Represents the dynamic functional connectivity features of empirical data, represents the dynamic functional connectivity features of the fitted data at step t; b. Convergence degree of spectrum characteristics: According to the fitting parameter values ​​G and Calculate the corresponding fitted data value; then, calculate At each node, the empirical data The maximum absolute value of the difference is used as a quantitative measure of the goodness of fit of the spectral feature. : , in, For parameters At the t-th iteration, the spectral feature index extracted from the fitting data is: is the spectral characteristic index of the empirical data; max means taking the maximum value; c. Parameter fit goodness: Simulations show that the bifurcation parameters Spectral characteristic indicators with empirical data There is a strong correlation between the two, and excessive or insufficient iterations in gradient descent may lead to a significant decrease in the spectral characteristic index of the fitting parameters and the empirical data. The correlation between the parameters and The correlation between and Maintaining high consistency during the iterations of gradient descent means that: In the absence of known actual parameter values, qualitative fitting performance is evaluated using the correlation between the fitted parameters and an empirical metric as part of the loss function: ; Based on the true loss function of the benchmark bifurcation parameter and the exploration of ways to quantify the goodness of fit from various angles, an approximate loss function of the gradient descent process is constructed to evaluate the degree of fit between the parameters and the data features, thereby obtaining the parameter combination with the minimum loss function value to ensure that the model exhibits good fitting ability and wide applicability in different application scenarios:

[0071] in:

[0072]

[0073]

[0074] d. Other standardization methods: When calculating the approximate loss function, different standardization methods can be used to optimize the performance of the model by adjusting the mean and standard deviation of the data features to have a uniform scale and distribution, including: i. Max-min normalization: scaling all features to a specific range (usually 0 to 1):

[0075] in, is the maximum value operation, is the minimum operation; ii. Mean standardization: Divide by the mean value of the feature to change the scale of the data, making the data scale uniform and easier for algorithm processing:

[0076] in, is the mean operation; iii. Robust scaling: Use the median and interquartile range (IQR) of the data to scale the data, which is more robust to outliers:

[0077] in, is the median operation; iv. Mean normalization: .

[0078] Preferably, step S900 specifically includes: During the gradient descent process, the parameters of each generation are used Generate new simulated data and update the tuning parameters by calculating the spectral characteristics , to better match the benchmark bifurcation parameter values, including: (1) Initial parameter setting: The optimal parameter combination obtained from the initial value selection in the gradient descent process and Start to proceed; (2) Substitute the dynamic equations of the brain model (such as the Watts-Strogatz (WS) network model or the Barabási-Albert (BA) network model) to generate simulation data: The parameters for the iteration Substitute the dynamic equations and solve them through numerical integration (such as Euler method, Runge-Kutta method) to obtain the time evolution of the system under given parameters and generate simulation data corresponding to each node. ;in , ; (3) Calculate the spectral characteristics of the simulated data: The generated simulated data Perform spectrum analysis and calculate ; The process includes performing a fast Fourier transform (FFT) on the analog data and then extracting the relevant spectral features fALFF (it should be noted that frequency center, bandwidth, spectral peak, spectral flatness, spectral roll-off and spectral entropy can also be used); (4) Calculate the spectral feature difference: calculate and The difference between them is used as the updated gradient of the parameters: ; (5) Calculate the data feature difference: ; (6) Parameter update: Using the difference in spectral features To update the parameters , so that it gradually approaches the value that can better fit the real data; Among them, the parameter update rule is based on a simple gradient descent method, and the updated parameters are expressed as:

[0079]

[0080] in, is the learning rate, which controls the step size of parameter updates and is calculated by the standard deviation of the indicator; (7) Iteration process: Repeat the above steps and calculate the value of the approximate loss function at each iteration; (8) Select the parameter value at the minimum loss function value: The bifurcation parameter value at the position where the approximate loss function value is the smallest is the value closest to the benchmark bifurcation parameter, that is, the better fitting parameter obtained.

[0081] In another embodiment, the present invention discloses a whole-brain dynamic parameter fitting method based on a Landau-Stewart oscillator network, comprising the following steps: S100, using the Landau-Stewart oscillator network to generate simulated data for subsequent model improvement and verification; wherein the simulated data is a blood oxygen level-dependent BOLD signal time series simulating real brain imaging data; Construct simulated brain networks, such as the connection structure based on the Watts-Strogatz (WS) network model, the connection structure based on the Barabási-Albert (BA) network model, and the connection structure based on the real brain network; It can be understood that the simulated brain network and simulated data are used for scientific research purposes such as verification and improvement of the model; S200, extracting data features of simulation data: extracting statistical features such as functional connectivity (FC), functional connectivity matrix at each moment (tFC), dynamic functional connectivity (FCD), Kolmogorov-Smirnov (KS) distance and probabilistic metastable substate (PMS) from the simulation data; S300, select global coupling parameters for fitting simulation data The initial value of Fix all bifurcation parameters to , a parameter sweep method was used to determine the optimal global coupling constant that minimizes the Kolmogorov-Smirnov (KS) distance of the FCD matrix fit of the data ; S400, select bifurcation parameters for fitting simulation data The initial value of Fix the optimal value determined in the previous step Value, yes exist The parameter was scanned with a step size of 0.01 in the range to find the best fit to the data FCD. In order to overcome the loss of individual information that may be caused by group averaging, this method estimates a set of individualized parameters for each subject. Therefore, a unique initial value is selected for each subject. , ensuring the personalization and accuracy of data processing.

[0082] S500, calculating a gradient descent index for fitting of simulation data; After detailed simulation analysis and evaluation, a convergence evaluation criterion was developed based on the fALFF index extracted from the data to guide the gradient descent optimization of the bifurcation parameters.

[0083] S600, calculating an adaptive learning rate for fitting of simulation data; According to the standard deviation of the fALFF indicator, an individualized adaptive learning rate strategy is selected for gradient descent, and the learning rate is adjusted according to the specific data characteristics of the subjects, thereby improving the efficiency and universality of model learning.

[0084] S700, implementing a gradient descent process and a gradient correction based on data features for fitting the simulation data; After thousands of simulations, we found the average value and parameters of the FCD characteristics of the data There is a positive correlation between the average value of and the correction measures for data feature fitting are introduced during the gradient descent process to accurately match the data features.

[0085] S800, constructing an approximate loss function for fitting the simulated data; During the gradient descent process, an approximate loss function is constructed based on the true loss function of the fitting benchmark bifurcation parameters to evaluate the degree of fit between the parameters and the data features, thereby obtaining the parameter combination with the minimum loss function value, ensuring that the model exhibits good fitting ability and wide applicability in different application scenarios.

[0086] S900, find the optimal parameters for fitting the simulation data: select a parameter combination that minimizes the approximate loss function, and this group of parameters is the optimal fitting parameters.

[0087] This means that the model parameters that perform best in the entire training process are selected. At this time, the model has the best fitting effect in parameter values ​​and correlations, data spectrum characteristics and data dynamic characteristics. Ultimately, the present invention can obtain the optimal parameters as the result of parameter fitting, which can be further used for prediction or analysis of the model.

[0088] For the above embodiment, specifically, in order to more accurately characterize and preserve the differences between data samples, we introduced several innovative individualization strategies, which not only enhanced the model's ability to capture the unique brain dynamics of each sample, but also improved the overall explanatory power of the model and the accuracy of the reconstruction of the BOLD signal time series through more refined parameter adjustments, among which, In the process of gradient descent, the present invention constructs an approximate loss function based on the characteristics of the real loss function of the similarity of parameter fitting and its accuracy that can be calculated when fitting the benchmark bifurcation parameter; the approximate loss function can effectively evaluate the degree of fit between the parameter and the data feature, thereby obtaining the parameter combination that minimizes the loss function value, so that the fitting method can improve the fitting accuracy and robustness of the parameters under different samples, more complex parameter value ranges and distributions. At the same time, the present invention introduces a gradient correction mechanism in the specific method, enhances the goodness of fit to the data features, thereby significantly improving the overall performance of the model, and greatly making up for the shortcomings of the existing fitting method.

[0089] It has been verified that the method of the present invention has shown higher accuracy, stability and generalization ability through comparison before and after the method is improved and comprehensive evaluation of more than 1,000 simulations, and has improved the individual adaptability of the model. The present invention not only demonstrates the comprehensiveness of the fitting method, but also highlights its innovation and practicality in precision modeling and parameter estimation, enabling it to exert greater potential in research on mental illness and the like.

[0090] The present invention has brought a major breakthrough in the parameter fitting of the Landau-Stewart oscillator network. Through simulation verification and data experiments, the method of the present invention can not only significantly improve the fitting performance of the model, but also has a wider applicability, and can open up new avenues for clinical and other applications.

[0091] It has been verified that the method disclosed by the present invention is applicable to brain network models of different scales (such as Figure 8 As shown, the fitting effect of the brain network model with 264 nodes is still better). The Landau-Stewart oscillator network is shown below:

[0092] The BOLD signal time series of each brain region was simulated, and the simulated data were generated in the model validation; In the real dataset, the brain region and brain regions The structural connection between is the structural connection matrix data of the human brain; in model verification, it is the node and nodes The structural connection between To generate simulated brain network data; represents the global coupling strength; is the bifurcation parameter, which is the node heterogeneity ; is a standard deviation Additive Gaussian noise; is the natural frequency of each brain region, located at Within the bandwidth, it is calculated from empirical data; In another embodiment, the parameter is and The fitting method includes the following steps: (1) Setting up the simulated brain network: The connection density of the simulated brain networks was set according to the connection density of the real brain networks; in these networks, the connection weights were uniformly distributed from Randomly selected from When the network is a 5-node network,

[0093] In order to ensure the diversity and unpredictability of the network, two network construction methods were selected during the test: a connection structure based on the Watts-Strogatz (WS) network model, a connection structure based on the Barabási-Albert (BA) network model, and a connection structure based on the real brain network; for the WS network model, the reconnection probability was changed by The authors introduced an additional random element by means of the probability of cutting and randomly reconnecting network edges, and experimented with multiple probability values. For the BA network model, they constructed scale-free networks generated at different connection densities, which are characterized by an uneven distribution of the number of connections between nodes, with some nodes having many connections and others having few. These simulated brain network architectures were used for method validation to confirm the effectiveness of the fitting strategy in dealing with complex networks. (2) Generate simulated data: A bifurcation parameter is randomly generated for each node as the benchmark bifurcation parameter value for evaluating the goodness of fit, and the corresponding parameter set is obtained as the verification benchmark, which is input into the Landau-Stewart oscillator network to generate the number of nodes. , the time length is The simulated data size is , simulate the scenario of inferring unknown parameters from experimental data in practical applications; Continue to connect the simulated brain network structure of the 5 nodes shown above and randomly generate bifurcation parameters: , and then generate simulated data: The generated simulated data are subjected to feature fitting and parameter estimation. By comparing the fitting results with the original parameters used to generate the simulated data, the accuracy and robustness of the parameter estimation method adopted are evaluated and continuously improved. If the fitted parameters are highly consistent with the parameters used to generate the data, it can be confirmed that the adopted method performs well in the simulation environment and is equally reliable when applied to actual data. This process not only provides a solid foundation for the development and improvement of parameter estimation technology, but also ensures its reliability and universal applicability in different research and clinical environments.

[0094] (3) Extracting data features The following key data features are extracted from the simulated data, and the model is used to fit the data features and used as a quantitative measure of the goodness of data fit: a. Functional connectivity (FC): The Pearson correlation coefficient between the simulated data of each pair of nodes is calculated to form a symmetric matrix. Each element in the matrix is ​​the correlation coefficient between the node pair, that is, the size of the functional connection:

[0095] in, and Is a node and nodes in The simulation data value at time, and The nodes are and nodes The mean of the simulated data; is the length of the simulated data time points; b. Functional connectivity matrix at each moment (dFC): Use the Hilbert transform to calculate each node ( ) Every moment ( ) of the simulated data , calculate the phase difference between each pair of nodes at each time point, and then use the cosine function to calculate the phase coherence matrix to obtain a size of The dFC matrix, where is the number of nodes, is the length of the simulated data time point: ; c. Dynamic functional connectivity (FCD) and KS distance: The FCD matrix based on phase congruency is a A symmetric matrix, where is the length of the sliding window time window; the elements in this matrix It is a time window and Functional connectivity matrix and The upper triangle part and The cosine similarity between: ; In order to compare the empirical and simulated FCDs, the Kolmogorov-Smirnov (KS) distance, which is used to quantify the maximum difference between the cumulative distribution functions of two samples, is used to measure the difference between the upper triangular elements of the empirical matrix and the simulated FCD matrix, where the range of the KS distance is , the smaller the KS distance is, the better the fitting degree of the data features is; d. Calculate the probabilistic metastable substate (PMS), including the following steps: i. Feature vector extraction: Extract the first feature vector from the functional connectivity (dFC) matrix at each moment ( ); The steps to extract the first eigenvector of the matrix are as follows: 1. Definition of matrix: Suppose there is a matrix , this is a The matrix represents samples and Features 2. Calculate the covariance matrix: First, calculate the matrix The covariance matrix of ;if Each column represents a feature and the row represents a sample, then the covariance matrix Calculated by: ,in yes The transposed matrix of 3. Calculate eigenvalues ​​and eigenvectors: Solve the covariance matrix The eigenvalues ​​and eigenvectors of satisfy the following equations: ,in is the eigenvalue, is the corresponding eigenvector; 4. Extract the first eigenvector: Sort the eigenvalues ​​from large to small. The eigenvector corresponding to the largest eigenvalue is the first eigenvector. ,in, Represents the direction with the largest directional variance in the data set; 5. Standardization: The feature vector Normalize it so that its modulus is 1, that is, its length is 1.

[0096] For example: Suppose there is a matrix:

[0097] 1. Calculate the covariance matrix: ; 2. Calculate eigenvalues ​​and eigenvectors: For , the eigenvalues ​​are 2 and 3. The corresponding eigenvectors are and ;

[0098] 3. Extract the first eigenvector: The maximum eigenvalue is 3, and the corresponding eigenvector is .

[0099] Therefore, the matrix The first eigenvector of yes .

[0100] ii. Cluster analysis: Cluster analysis is performed using the extracted feature vectors to discover common and recurring connectivity patterns across different time points and samples. The k-means clustering algorithm is used to cluster data points (i.e., dFC matrices at different time points) into several clusters based on the similarity of feature vectors. iii. Define metastable substates: Each cluster in the clustering result represents a unique metastable substate, and the centroid of the cluster (i.e., the cluster center) is defined as the representative metastable substate; iv. Probability and transition probability calculation: For each metastable substate, calculate its probability of appearing at all time points, and analyze the transition probability between different substates, that is, the probability of converting from one substate to another; v. Measure goodness of fit: Use KL distance to compare the probability of metastable states obtained from simulation and empirical (experimental) data. Metastable substate The empirical probability of is the simulation probability on the same metastable substate; the KL distance between the empirical data and the simulated data is calculated to evaluate the fit of the simulation model to the actually observed brain state probability distribution: ; (4) Select global coupling parameters The initial value of: Fix all bifurcation parameters to , a parameter sweep method was used to determine the optimal global coupling constant that minimizes the Kolmogorov-Smirnov (KS) distance of the FCD matrix fit of the data .

[0101] (5) Select bifurcation parameters The initial value of: The optimal Based on the value, exist The parameter was scanned with a step size of 0.01 in the range to find the best fit to the data FCD. The initial value of is the parameter value with the smallest KS distance. In order to overcome the loss of individual information that may be caused by group averaging, this method estimates a set of individualized parameters for each subject. Therefore, a unique initial value is selected for each subject. , ensuring the personalization and accuracy of data processing.

[0102] (6) Calculate the gradient descent index The spectral characteristics of each node extracted from the data are selected as the key indicator to characterize the numerical value of the bifurcation parameter of each node; this indicator is used to formulate the basis and evaluation criteria for parameter updating during the gradient descent process to effectively guide the adjustment of the bifurcation parameter during the gradient descent optimization process; this indicator is node heterogeneous, that is, each element of the indicator vector is the spectral characteristic of each brain region, expressed as The process of calculating this indicator using simulated data includes: a. Demeaning: Demean the simulated data of each node so that the mean of the processed simulated data at all time points on each node is zero; this process helps to remove fixed deviations in the data, thereby better analyzing its fluctuations and changes; b. Filtering: Filter the simulated data after removing the mean value of each node to eliminate possible linear changes in the data and ensure the accuracy of the analysis.

[0103] c. Detrend: Detrend the simulated data after removing the mean value to eliminate possible linear changes in the data and ensure the accuracy of the analysis; d. Fast Fourier Transform (FFT): Perform fast Fourier transform on the simulated data after de-linear trending to convert it from the time domain to the frequency domain, thereby obtaining the spectrum characteristics of each frequency component; e. Calculate the proportion of low-frequency spectrum features of each node of the simulated data: Select the low-frequency band range where the BOLD signal is located, calculate the value that can characterize the spectrum features of the simulated data within this range, and calculate it by dividing the sum of the spectrum features of this band by the sum of the spectrum features of the entire band:

[0104] In digital signal processing, the full frequency range is usually defined by the Nyquist frequency, which is half the sampling frequency. According to the Nyquist theorem, in order to unambiguously capture and reconstruct all frequency components in a signal, the signal must be sampled at least twice the rate of the highest frequency component. When TR = 2s, the full frequency range is defined as ; f. Calculate the proportion of spectrum features in different frequency bands for each node of the simulation data: In addition to the low-frequency band where the BOLD signal is located, any frequency band can be selected and the ratio of its spectrum features to the spectrum features of the full frequency band can be calculated as a way to characterize the node dynamics. In addition, this ratio can be processed in various ways, including standardization, exponential or multiple scaling, to obtain an indicator that is more suitable for specific analysis or application, so as to better characterize the size of the bifurcation parameters of each node. g. Calculate the full-band spectrum characteristics of each node of the simulated data: calculate the frequency center, bandwidth, spectrum flatness, spectrum roll-off, and spectrum entropy of the full-band to characterize the size of the bifurcation parameters of each node. Among them, the low-frequency amplitude ratio (fALFF) extracted from the data is preferentially selected as the indicator for calculating the gradient during the gradient descent process to achieve parameter update and improve the goodness of fitting of the data spectrum characteristics, so as to effectively guide the adjustment of the bifurcation parameters during the gradient descent optimization process. Among them, fALFF focuses on the amplitude of the low-frequency signal, that is, the ratio of the amplitude of the signal in the low-frequency signal band where the BOLD signal is located (generally 0.01-0.1Hz, this frequency range is considered to be most relevant to neural activity in the resting state) to the amplitude of the full-band signal (0-0.25Hz when TR=2s), which is the node heterogeneity and is expressed as ; For example, taking a network with 5 nodes as an example, the generated structural connection matrix is:

[0105] In the parameters The simulated data generated for: Follow the steps below to calculate fALFF for this simulation data: a. Demeaning: First, the original simulated data of each node is demeaned so that the mean of the processed simulated data at all time points on each node is zero, that is, the original simulated data after demeaning is:

[0106] b. Filtering: Filter the simulated data after removing the mean value of each node to eliminate possible linear changes in the data and ensure the accuracy of the analysis; c. Detrend: Detrend the simulated data after removing the mean of each node to eliminate possible linear changes in the data and ensure the accuracy of the analysis:

[0107] d. Fast Fourier Transform (FFT): Perform fast Fourier transform on the simulated data after de-linearization to convert the simulated data from the time domain to the frequency domain, and take the modulus of the FFT result to obtain the amplitude of each frequency component, that is, convert it into:

[0108] e. Calculate the amplitude characteristics of the low-frequency range: Select a low-frequency range (e.g. ), calculate the amplitude and value within this range: ; f. Calculate the full-band amplitude characteristics: The full-band is usually defined by the Nyquist frequency, which is half the sampling frequency. According to the Nyquist theorem, in order to unambiguously capture and reconstruct all frequency components in the data, the data must be sampled at least twice the rate of the highest frequency component in the data. That is, when TR = 2s, the full-band is defined as : ; g. Calculate fALFF: It is calculated by dividing the sum of the low-frequency amplitudes by the sum of the full-band amplitudes. The formula is:

[0109] Finally, the simulated data is calculated as follows: [0.8254, 0.8714, 0.4526, 0.4930, 0.5927], at this time the index and parameter The correlation between ; In addition to the 0.01-0.1Hz frequency band in step (5) of S500, "Calculate the proportion of low-frequency spectrum features of each node in the simulated data", any low-frequency band can also be selected, and the amplitude and ratio with the full frequency band can be calculated as the index. When other frequency bands are selected, the remaining steps are the same, only the frequency band selection in step (5) is changed. After simulation verification on the generated model, it is found that when the low-frequency band is selected as half of the full frequency band, the correlation between the index and the parameter may be higher, and the fitting performance is better. When TR=2s, the full frequency band is , the low frequency band can be selected as , the indicator at this time Calculated as: For this frequency band, the simulated data in the example calculates: , at this time, the indicator and parameter The correlation between ; In addition to the commonly used fALFF indicator in imaging research, the following indicators can also reflect the dynamic characteristics of the node from different angles, and can characterize the value of the bifurcation parameter to a certain extent. Therefore, it can also be used as an indicator to guide the direction of gradient descent and the learning rate. When using other indicators that characterize spectral features, steps (1) to (4) in the S500 step remain unchanged, and the sum of the amplitude features in steps (5) to (7) is replaced by the sum of other features, and the method remains unchanged. Other indicators that can be used are: a. Low-frequency amplitude fluctuation (ALFF): It is a functional magnetic resonance imaging (fMRI) data analysis method used to measure the amplitude of the low-frequency band of the brain at rest. It is an important indicator for studying spontaneous neural activity in the brain and reflects the intensity of low-frequency fluctuations in the brain at rest. The calculation process is to first select a specific low-frequency band, such as , calculate the amplitude and value in this frequency range:

[0110] in, It is The amplitude of the frequency component, is the total number of frequency components in the selected low-frequency band; since the overall signal strength of different nodes may vary greatly, directly comparing the ALFF values ​​of different nodes may introduce bias and mask the actual functional differences. After replacing the amplitude in steps (5) to (7) in step S500 with this variable, the simulated data in the example calculates the low-frequency amplitude fluctuation as the gradient descent indicator, that is: , at this time, the indicator and parameter The correlation between them is poor and unstable: ; b. Low-frequency power ratio: The power of each frequency directly reflects the contribution of the frequency component in the signal; this analysis provides the distribution of signal energy in the frequency domain, which helps to understand the characteristics of the signal from different angles; calculate a certain low-frequency band, such as , the ratio of the power sum to the power sum of the entire frequency band is used as an indicator: ,

[0111] in, is the total number of frequency components in the selected low-frequency band, It is frequency components, is the power on this frequency component; after replacing the amplitude in steps (5) to (7) in step S500 with this variable, the simulated data in the example calculates the low-frequency power ratio as the gradient descent indicator, that is: , with the parameter Relevance ; c. SpectralCentroid: It is the "center of gravity" of the signal spectrum, used to indicate the central position of the signal spectrum weight distribution; it indicates the average frequency, or "center of gravity" of the oscillator motion spectrum; a high frequency center may indicate that the oscillator is more active in a higher frequency range, which should be consistent with the bifurcation parameter Inversely proportional; therefore, when used to calculate the gradient in gradient descent, the indicator should be taken as the negative value of the frequency center, and the calculation method is as follows: ,

[0112] in, It is frequency components, is the power at that frequency component, It is After replacing the amplitude in steps (5) to (7) in step S500 with this variable, the frequency center is calculated as the gradient descent indicator for the simulated data in the example, that is: , with the parameter Relevance ; d. Spectral Bandwidth: It is used to describe the spectral distribution width of the signal. It is the distribution width of the energy around the center of the spectrum and represents the width of the oscillator motion spectrum. The larger the spectral bandwidth of the oscillator motion, the wider the distribution of its frequency components, the more complex the motion, and the vibrations of multiple frequencies. The smaller the bandwidth value of the oscillator, the more it vibrates within a narrow frequency range (such as simple harmonic vibration). The larger the bandwidth value, the wider the spectrum (such as forced vibration or nonlinear vibration). It should be related to the bifurcation parameter. is inversely proportional, so The indicator should be taken as the negative value of the bandwidth; the bandwidth calculation method is based on the frequency center: ,

[0113] in, It is frequency components, is the power at that frequency component, It is The frequency center of the node, It is The bandwidth of the nodes; after replacing the amplitude in steps (5) to (7) in step S500 with this variable, the bandwidth is calculated as the gradient descent indicator for the simulated data in the example, that is: , with the parameter Relevance ; e. Spectral Flatness: reflects the shape of the spectrum and determines whether the spectrum of the oscillator is close to noise. Low spectral flatness indicates that the energy of the oscillator is concentrated on certain specific frequencies, which is usually related to resonance or forced response. High spectral flatness (close to 1) indicates that the energy is more evenly distributed at each frequency and the vibration is closer to random noise. Low spectral flatness (close to 0) indicates that the signal contains significant tones or harmonic components. Spectral flatness is defined as the ratio of the geometric mean to the arithmetic mean of the amplitude distribution of the signal, which is related to the bifurcation parameter. is inversely proportional, so The indicator should take its negative value: ,

[0114] in, It is The power of the frequency components, It is The spectral flatness of the node; after replacing the amplitude in steps (5) to (7) in step S500 with this variable, the spectral flatness of the simulated data in the example is calculated as the gradient descent indicator, that is: [-0.0812, -0.0342, -0.8881, -0.8228, -0.6097], with parameters Relevance ; f. Spectral Roll-off: Spectral roll-off refers to the frequency point below a certain percentage (usually 85% or 90%, 85% is selected here), below which most of the signal energy is contained, indicating the upper limit of the frequency where most of the energy is concentrated; this helps to determine the energy concentration area of ​​the oscillator response and understand in which frequency range the oscillator is mainly active. The calculation method is that the cumulative energy reaches the roll-off frequency of 85% or 90% of the total energy, and the bifurcation parameter is inversely proportional, so The indicator should take its negative value: ,

[0115] in, It is The power of a frequency component, is the roll-off percentage, It is After replacing the amplitude in steps (5) to (7) in step S500 with this variable, the spectrum roll-off is calculated as the gradient descent indicator for the simulated data in the example, that is: [-0.0625, -0.0575, -0.1850, -0.1600, -0.1250], with parameters Relevance ; g. Spectral Entropy: Spectral Entropy is based on the concept of information theory and is used to measure the uniformity of the signal's spectral distribution. It is the entropy of the normalized distribution of the signal's spectral amplitude. A high spectral entropy indicates that the signal is more evenly distributed on the spectrum, which is usually related to noise. A low spectral entropy indicates that the signal's energy is concentrated on a specific frequency, which is also related to the bifurcation parameter. is inversely proportional, so The indicator should take its negative value and bifurcation parameter is inversely proportional, so The indicator should take its negative value: ,

[0116] in, It is The power of a frequency component, It is The spectral entropy of the node; after replacing the amplitude in steps (5) to (7) in step S500 with this variable, the spectral entropy of the simulated data in the example is calculated as the gradient descent indicator, that is: [-3.1257, -2.9633, -4.5079, -4.4402, -4.1635], with the parameters Relevance .

[0117] (7) Calculate adaptive learning rate According to the standard deviation, mean, amplitude and other characteristics of the indicators calculated in step (6), an individualized adaptive learning rate strategy is selected for gradient descent, and the learning rate is adjusted according to the specific data characteristics of the subjects, thereby improving the efficiency and universality of model learning; The inventors found in simulation that the parameter range is proportional to the standard deviation of the spectral characteristic index. Therefore, in the method of the present invention, the standard deviation is selected to adjust the iterative learning rate of the gradient descent:

[0118] in, is the coefficient of the learning rate, which can be selected according to the data characteristics. For example, .

[0119] (8) Gradient correction of data features: After thousands of simulations, we found the average value and parameters of the FCD characteristics of the data The average value of is positively correlated, and the correction measures of data feature fitting are introduced in the gradient descent process to accurately match the data features:

[0120] (9) Based on the parameters of the generated simulated data and the parameters obtained by fitting the simulated data, construct the true loss function: In order to more accurately determine whether the iteration of gradient descent has reached the target value, the present invention introduces a loss function as an evaluation criterion; in the experiment, the bifurcation parameter value of the generated data It is known that as the benchmark bifurcation parameter value for evaluating the goodness of fit, it can be used to construct the true loss function to verify the performance of parameter fitting in quantitative (difference) and qualitative (correlation) aspects; based on this, the true loss function consists of two parts: a. Difference between fitted bifurcation parameters and benchmark bifurcation parameters: ; b. Similarity between the fitted bifurcation parameters and the benchmark bifurcation parameters: ; in, represents the base bifurcation parameter value for generating data, It is step The fitting parameter value at the time of ; In order to ensure that the two components contribute equally to the overall loss, so as to make a more comprehensive and accurate evaluation of the performance of the model, the loss value of each part is converted into the deviation ratio relative to its best performance in the iteration process; The true loss function calculation formula is defined as follows: ; (10) Based on the characteristics of the simulated data, construct an approximate loss function: In the process of gradient descent, the approximate loss function can quantify the goodness of fit between parameters and data features, and parameter features, so that the parameter combination with the minimum loss function value is closest to the benchmark bifurcation parameter value, ensuring that the model exhibits good fitting ability and wide applicability in different application scenarios; Based on the real loss function when fitting the benchmark bifurcation parameters and the exploration of ways to quantify the goodness of fit from various angles, an approximate loss function consisting of three parts: the goodness of fit of data features, the degree of convergence of spectral features, and the goodness of fit of parameters is constructed to evaluate the degree of fit between parameters and data features, thereby obtaining the parameter combination with the minimum loss function value: ; a. Goodness of fit of data characteristics: According to the fitting requirements of the data characteristics, the corresponding quantifiable goodness of fit indicators are selected, such as the correlation between the simulated and real FC matrices, the KS distance between the FCD matrices, and the KL distance between the PMS vectors. t ; b. Convergence degree of spectrum characteristics: Substitute the fitting parameter values ​​obtained in each iteration into the dynamic equation of the brain model and calculate the corresponding fitting Then, calculate its difference with the actual The difference between the values, and the maximum absolute value of the difference is taken as a quantitative measure of the goodness of fit of the spectral feature: ; c. Parameter fit goodness: Simulation results show that: simulation parameters and spectrum characteristic indicators There is a strong correlation between the two, and excessive or insufficient iterations in gradient descent may lead to a significant decrease in the parameters and the empirical index Further simulations show that the simulation parameters are consistent with the benchmark parameters. The correlation between and The high consistency maintained during the iterations of gradient descent means that in the absence of known actual parameter values, the correlation between simulated parameters and empirical metrics can be used as part of the loss function to evaluate qualitative fitting performance: ; Based on the true loss function when fitting the benchmark bifurcation parameters and the exploration of ways to quantify the goodness of fit from various angles, a gradient descent process approximate loss function is constructed to evaluate the degree of fit between the parameters and the data features, thereby obtaining the parameter combination with the minimum loss function value, ensuring that the model exhibits good fitting ability and wide applicability in different application scenarios:

[0121] in:

[0122]

[0123]

[0124] The approximate loss function can better approximate the changing trend of the real loss function; d. Other standardization methods: When calculating the approximate loss function, different standardization methods can be used to adjust the mean and standard deviation of data features to make them have a uniform scale and distribution, thereby optimizing the performance of the model, including: i. Max-min normalization: scaling all features to a specific range (usually 0 to 1):

[0125] ii. Mean standardization: Divide by the mean value of the feature to change the scale of the data, making the data scale uniform and easier for algorithm processing:

[0126] iii. Robust scaling: Use the median and interquartile range (IQR) of the data to scale the data, which is more robust to outliers:

[0127] iv. Mean normalization:

[0128] (11) Gradient descent process: In this gradient descent process, the present invention uses the parameters of each generation The dynamic model generates time series and updates the adjustment parameters by calculating the spectral characteristics. , to better match the benchmark bifurcation parameter value. The flow chart is as follows Figure 6 The detailed description process is as follows: a. Initial parameter setting: The optimal parameter combination obtained from the initial value selection during the gradient descent process and Get started.

[0129] b. Substitute the dynamic equations of the brain model to simulate the time series: The parameters for the iteration Substitute the dynamic equations of the brain model into the differential equations, solve the differential equations through numerical integration (such as Euler method, Runge-Kutta method) and other techniques, obtain the time evolution of the system under given parameters, and generate the time series corresponding to each node .in ,

[0130] c. Calculate the spectral characteristics of the simulated data: generate the time series Perform spectrum analysis and calculate its spectrum characteristic indicators The process involves performing a Fast Fourier Transform (FFT) on the time series and then extracting the relevant spectral features fALFF (frequency center, bandwidth, spectral peak, spectral flatness, spectral roll-off, spectral entropy, etc. can also be used).

[0131] d. Calculate the difference in spectral characteristics: Calculate the spectral characteristics of the simulated data Spectral characteristics of real data The difference between them is used as the updated gradient of the parameters: .

[0132] e. Calculate the data feature difference:

[0133] f. Parameter update: using the difference of spectral features To update the parameters , so that it gradually approaches the value that can better fit the real data.

[0134] The update rule is based on a simple gradient descent method, and the updated parameters are expressed as:

[0135]

[0136] in, is the learning rate, which controls the step size of parameter update and is calculated by the standard deviation of the indicator.

[0137] g. Iteration process: Repeat the above steps and calculate the value of the approximate loss function at each iteration.

[0138] h. Select the parameter value at the minimum loss function value: The bifurcation parameter value at the position where the approximate loss function value is the minimum is the closest to the benchmark bifurcation parameter value, that is, the better fitting parameter we get. The fitting effect is as follows Figure 7 shown.

[0139] It can be understood that the model of the present invention models the BOLD signal time series in brain imaging data. Therefore, when facing real brain imaging data, before performing parameter fitting on the subjects to establish an individualized model, the functional magnetic resonance imaging (fMRI) image data must be processed first to extract the blood oxygen level-dependent BOLD signal time series data required for model fitting.

[0140] Therefore, in another embodiment, the present invention also discloses a method for extracting blood oxygen level dependent BOLD signal time series data, comprising the following steps: S1, data preprocessing; Wherein, the data preprocessing is used to correct various potential influencing factors in fMRI data to make the data suitable for further analysis; S2, spatial standardization; Wherein, the spatial normalization is used to convert fMRI data into standard spatial coordinates to ensure that signals of different samples can be compared; S3, define template; S4, extract BOLD time series; S5. Remove confounding factors; It is used to remove possible confounding factors from the extracted time series.

[0141] Preferably, in step S1, the data preprocessing includes: S11. Timing correction: Correct the timing of different slices collected at different times during the fMRI scan; S12, motion correction: correcting or compensating for the signal deviation caused by the head movement; S13, spatial registration: aligning functional images with anatomical images and mapping different samples to the same standard template; S14, spatial smoothing: Improve the signal-to-noise ratio by smoothing the data.

[0142] Preferably, in step S2, Standard templates used include the Montreal Neurological Institute (MNI) template or the Talairach template to map sample images to the standard template.

[0143] Preferably, step S3 specifically includes: According to the research objectives of interest, a brain region template is defined or selected, and each brain region in the template corresponds to each node in the subsequent model; among them, the brain region segmentation templates used include AAL template, Harvard-Oxford template, and Dosenbach template.

[0144] Preferably, step S4 specifically includes: Use the brain region mask in the template to perform segmentation in the preprocessed and normalized fMRI images; The average BOLD signal time series of all voxels within each brain region defined in the template was calculated to obtain a time series representing the region, corresponding to the time series of each node in the subsequent model.

[0145] Preferably, Confounding factors include head movement, whole-brain average signal, CSF signal, etc. Regression analysis was used to remove confounding factors to ensure that the extracted time series more accurately reflected neural activity.

[0146] Preferably, Fix all bifurcation parameters to , a parameter scanning method is used to determine the optimal global coupling constant that minimizes the Kolmogorov-Smirnov (KS) distance of the FCD matrix fitting of the time series .

[0147] Preferably, in step S400, Fixed optimal Value, yes exist The parameter scan with a step size of 0.01 was performed in the range to find the best fit for the BOLD signal time series FCD. The initial value of .

[0148] Preferably, In order to overcome the possible loss of individual information caused by group averaging, a set of individualized parameters is estimated for each subject.

[0149] Preferably, Select the initial value suitable for each subject , ensuring the personalization and accuracy of BOLD signal time series processing.

[0150] Preferably, in step S500, A convergence evaluation criterion was developed based on the low-frequency amplitude fraction (fALFF) indicator extracted from the BOLD signal time series to guide the gradient descent optimization of the bifurcation parameters.

[0151] Preferably, in step S600, According to the standard deviation of the fALFF index, an individualized adaptive learning rate strategy is selected for gradient descent, and the learning rate is adjusted according to the specific BOLD signal time series characteristics of the subjects to improve the efficiency and universality of model learning.

[0152] Preferably, in step S800, During the gradient descent process, an approximate loss function is constructed based on the true loss function of the benchmark bifurcation parameters to evaluate the degree of fit between the parameters and the time series characteristics of the BOLD signal, thereby obtaining the parameter combination with the minimum loss function value and ensuring that the model exhibits good fitting ability and wide applicability in different application scenarios.

[0153] Preferably, in step S900, Select the parameter combination that minimizes the approximate loss function, and this set of parameters is the optimal fitting parameters.

[0154] In another embodiment, the present invention also shows how to apply it in a real data set, that is, using a method for extracting blood oxygen level dependent BOLD signal time series data from fMRI imaging, comprising the following steps: S11. Data preprocessing: The data preprocessing is used to correct various potential influencing factors in fMRI data to make the data suitable for further analysis; Exemplarily, the data preprocessing mainly includes: S111. Timing correction: During fMRI scanning, different slices are not collected at the same time and need to be corrected.

[0155] S112. Motion correction: Since head movement can cause signal deviation, this movement needs to be corrected or compensated.

[0156] S113. Spatial registration: Align functional images with anatomical images and map different samples to the same standard template.

[0157] S114, spatial smoothing: Improve the signal-to-noise ratio by smoothing the data.

[0158] S12. Spatial standardization: The spatial normalization is used to convert fMRI data to standard spatial coordinates to ensure that signals from different samples can be compared. Commonly used templates include MNI (Montreal Neurological Institute) or Talairach templates, and sample images are mapped to these standard templates.

[0159] S13. Define template: Define or select a brain region template according to the research objectives of interest. Each brain region in the template corresponds to each node in the subsequent model. You can use ready-made brain region segmentation templates, such as AAL (90 brain regions in the cerebral cortex, and subsequent explanations are based on this template), Harvard-Oxford (48 cortical regions and 21 subcortical regions), Dosenbach (160 brain regions), and other atlases.

[0160] S14. Extract BOLD time series: Using the brain region mask in the template, segmentation is performed in the preprocessed and normalized fMRI images. The average BOLD signal of all voxels within each brain region defined in the template is calculated to obtain a time series representing the region, which corresponds to the time series of each node in the subsequent model.

[0161] S15. Remove confounding factors: Remove possible confounding factors (such as head motion, whole-brain average signal, CSF signal, etc.) from the extracted time series. Regression analysis is used to remove these confounding factors to ensure that the extracted time series more accurately reflects neural activity.

[0162] The fitting process of the real data set follows the above process, wherein the structural connection between the nodes used in the equation is replaced by the real human brain structural connection data to simulate the structure of the brain network, and the simulated data is replaced by the BOLD signal time series extracted after steps S11-S15, specifically the following steps: S16. Extract data features for fitting of BOLD signal time series: extract features from the time series data of the sample: functional connectivity (FC), functional connectivity matrix at each moment (tFC), dynamic functional connectivity (FCD) and Kolmogorov-Smirnov (KS) distance and probabilistic metastable substate (PMS); S17. For the fitting of BOLD signal time series, select the global coupling parameters The initial value of S18. Select bifurcation parameters for fitting the BOLD signal time series The initial value of S19, calculating the gradient descent index for the fitting of the BOLD signal time series; S20, calculating the adaptive learning rate for the fitting of the BOLD signal time series; S21, for the fitting of BOLD signal time series, implement the gradient descent process and gradient correction based on data characteristics; S22, construct an approximate loss function for fitting the BOLD signal time series; S23. Find the optimal parameters for fitting the BOLD signal time series; S24. Analyze and discuss the optimal parameters.

[0163] Here, the complete application and analysis process of the present invention in a real data set is demonstrated: It should be noted that the example shows the public dataset "Depression (MDD) Dataset" from the REST-meta-MDD project of the DIRECT Alliance; This dataset is composed of resting-state functional magnetic resonance imaging (R-fMRI) data of patients with major depressive disorder (MDD) and matched normal controls (HC) collected by 25 research groups from 17 hospitals in China. There are 2428 samples in total, including 1300 MDD samples and 1128 HC samples. After excluding samples that did not meet the age standard (18-65 years old), incomplete information, low image quality (by visual inspection), poor spatial normalization, insufficient coverage (<0.9), short sampling time (time step <190), low depression score (HAMD <7), inappropriate TR parameters, excessive head movement (mean FD>0.2mm), and abnormal fALFF index, 666 MDD and 648 HC were selected for in-depth parameter fitting and statistical analysis; Among these MDD patients, females accounted for the majority, with a total of 427 cases, compared with 239 male patients. Among them, there were 138 recurrent MDD and 317 first-episode MDD (159 first-episode unmedicated (FEDN) and 100 first-episode medicated (FEMed) MDD). Regarding medication status, 207 MDD patients were not taking medication, and 179 MDD patients were taking medication. The specific groups and numbers are shown in Table 1, where Table 1 illustrates the specific groups and numbers of MDD in a real data set application embodiment of the present invention:

[0164] Table 1 Different groups and their number of people In addition, the shared information included: subject ID, gender, age, education level, seizure status, medication status, disease course, Hamilton Depression Rating Scale (HAMD) and Hamilton Anxiety Rating Scale (HAMA); After performing steps S11-S15 on this MDD data set, the BOLD signal time series is obtained; The improved and verified model parameter estimation method of the present invention is applied to fitting the BOLD signal time series, first completing step S16 to extract the distribution of the data feature FCD of the BOLD signal time series; Then proceed to step S17: fix the bifurcation parameters to be Through parameter scanning, we determined the optimal global coupling parameters for each sample as ,like Fig. 9 As shown in the left figure of (a); Step S18: Fixing Then, according to the initial value selection method, the bifurcation parameters of each sample are scanned The initial value of Fig. 9 As shown in the right figure of (a); Using these initial values, the adaptive method was applied to iterate for 666 MDD and 648 HC samples through steps S19-S23, and finally an individualized model was established, and the bifurcation parameter values ​​of the brain region heterogeneity of each sample were fitted; The parameters of the individualized model finally obtained according to the present invention approximate the target value under the constraint of the loss function, effectively fit the BOLD data characteristics, and can reveal the dynamic characteristics of the brain regions of each patient under different pathological conditions. For example, when fitting the functional connectivity density (FCD) matrix, the average KS distance of the MDD group was 0.0820 (standard deviation = 0.0319), and the average KS distance of the HC group was 0.0812 (standard deviation = 0.0307) (e.g. Fig. 9 (b) shows a better fitting performance; To analyze the differences between the two groups in bifurcation parameters, a t-test was performed on the KS fit index of the HC group and the MDD group to ensure that the differences were not caused by differences in fitting accuracy. The results showed that there was no significant difference in goodness of fit between the two groups (p = 0.7525). Fig. 9 Based on this bifurcation parameter, we conduct relevant statistical analysis, classification analysis and discussion as described in S24, such as: a. Classification analysis: According to the specific condition information of each MDD patient, they were divided into different groups (see Table 1 for the groups and their numbers), and 10 classification experiments were carried out. In each experiment, the bifurcation parameter vectors of the participants were used as input data for linear and radial basis function (RBF) support vector machines (SVMs), and appropriate hyperparameters were selected for each model (Note: Appropriate hyperparameters can be selected based on existing technologies). By using the GridSearchCV tool in the Scikit-learn software package (Note: This tool is a hyperparameter adjustment tool), all possible parameter combinations were evaluated, with the aim of finding the best parameter configuration during the cross-validation process to ensure that the model achieves the highest balanced accuracy. After determining the best hyperparameters, the SVM model was retrained on the entire training dataset, and its performance on the test dataset was evaluated. Key performance indicators such as accuracy, precision, recall, and F1 score of the model were calculated, and a five-fold cross-validation strategy was used to calculate the mean and standard deviation of these performance indicators; The final classification results show that the bifurcation parameter generally provides higher classification performance, as shown in Table 2 below, wherein Table 2 illustrates the classification accuracy of medicated and unmedicated depression (MDD) patients using functional connectivity (FC) or bifurcation parameters as features in the same classifier in an embodiment of the present invention applied to a real data set:

[0165] Table 2 Comparison of classification results b. Test and analysis of significant differences between patients and healthy controls: A t-test was performed for each brain region in HC and MDD, and FDR multiple comparison correction was performed to identify brain regions with significant differences. In addition to p values ​​less than 0.05, the effect size index (Cohen's d) was calculated, and finally 11 brain regions with effect sizes d greater than 0.25 were identified, such as Fig.10 As shown in the asterisk part in (a). Among them, three brain regions in the top five feature weights in the SVM classifier are very significant in terms of statistics and effect size (such as Fig.10 (a) shows the overlap of the asterisk and shaded area), specifically the right hippocampus ( , ), left supplementary motor area ( , ) and the right supplementary motor area ( , ). This indicates that these brain regions have great potential as biomarkers for distinguishing HC from MDD; In general, the study found that in patients with depression, the bifurcation parameters of the hippocampus were more negative, indicating that this area has lower activity levels and is insensitive to changes in response to external or internal stimuli. The hippocampus is a key area in the brain responsible for memory, learning, and emotional regulation. Its reduced activity level may manifest as functional inhibition, which may impair its normal regulation of emotional and cognitive processes, thereby potentially triggering depressive symptoms.

[0166] Therefore, it can be speculated that this phenomenon may be due to atrophy or volume reduction of brain regions. In fact, many MRI and fMRI studies have also confirmed that patients with depression have significantly reduced gray matter volume in the hippocampus and altered functional brain responses to emotional stimuli, which is associated with impaired processing of negative emotions and cognitive functions.

[0167] c. Significant correlation analysis with severity scores: The widely accepted scoring criteria of the Hamilton Depression Rating Scale (HAMD) are: a total score below 7 indicates normal; 7 to 17 points indicate possible depression; 17 to 24 points confirm depression; and more than 24 points indicate severe depression. For this, linear mixed models (LMMs) were used to deeply explore the correlation between the bifurcation parameters of each brain region and the HAMD score. Among them, the score of each sample was set as the dependent variable, the bifurcation parameter of the node was used as a fixed effect, the intercept affected by site differences was used as a random effect, and age, gender, education, and head motion were used as covariates to control potential confounding factors in the form of: , For this, the fixed effect produced a value that was less than the threshold after multiple comparison correction (FDR) and other covariates, and it was considered that there was a significant correlation between the bifurcation parameter value of the node and the HAMD score; In the severely depressed group (HAMD>24), there was a significant correlation between the bifurcation parameters of the left thalamus and the HAMD score, and this correlation remained after correction for multiple comparisons (FDR) of 90 brain regions (LMM: , , (for the entire dataset), (For the average data of each scoring point, FDR was taken into account). In addition, these parameter values ​​showed significant differences compared with the HC group (t-test: , Cohen's 0.2577, after FDR), such as Fig.11 Middle (a), Fig.11 This finding highlights the key role of the thalamus in MDDs. Detailed trends in thalamic bifurcation parameter values ​​as scores change can be seen in Fig.11 middle (d); In addition, in order to identify differences in the dynamic characteristics of brain regions between MDD patients with different severities and HCs, a comprehensive comparative analysis was performed for each brain region using t-tests and FDR correction. The study found that in more severe MDD patients (HAMD>24), there were more significant differences and larger effect sizes in the nodes of the limbic network and frontal regions of the brain, such as the rectal gyrus, hippocampus, and thalamus. Among them, the limbic system plays a vital role in emotion regulation, memory formation, and coping with stress. These findings reveal a potential link between severe depression and dysfunction of the brain's limbic system, such as Fig.11 Middle (a), Fig.11 (c) shows the nodes with significant differences between the two groups after t-test and FDR correction; Fig.11 In (e), the effect sizes of these differences are visualized, providing a clearer understanding of changes in different brain regions.

[0168] About the application of the present invention: It can be understood that the scheme disclosed by the present invention can be used to set personalized parameters for each subject in a real data set, so as to better analyze the differences in brain regions between patients and healthy controls, capture the dynamic changes in brain function, and control or guide clinical practice.

[0169] Although the embodiments of the present invention are described above in conjunction with the accompanying drawings, the present invention is not limited to the above specific embodiments and application fields, and the above specific embodiments are only illustrative and instructive, rather than restrictive. A person of ordinary skill in the art can also make many forms under the guidance of this specification and without departing from the scope of protection of the claims of the present invention, all of which belong to the protection of the present invention.

Claims

1. A whole-brain dynamic parameter fitting method based on a Landau-Stewart oscillator network, comprising the following steps: S100, generating time series data of blood oxygen level dependence on BOLD signal in simulated real brain imaging data as simulation data, wherein: The simulation data are used for the global coupling parameters of the Landau-Stewart oscillator network and bifurcation parameters The fitting of S200, extracting data features of simulation data; S300, select global coupling parameters based on the data characteristics of the simulation data The initial value of S400, select bifurcation parameters based on the data characteristics of the simulation data The initial value of S500, calculating a gradient descent index of the simulation data according to the data features of the simulation data; S600, calculating an adaptive learning rate according to data features of the simulation data; S700, for the data features of the simulation data, introducing a data feature correction value for gradient correction of the data features; wherein, Introduced during gradient descent The data feature correction value is used as a correction measure for data feature fitting to accurately match the data features: ； in, The subscript t represents the bifurcation parameter Iterate the tth step. The superscript sim indicates the simulation value, and the superscript emp indicates the empirical value. Represents the bifurcation parameter The data feature correction value at the t-th iteration step, Represents the bifurcation parameter At the t-th iteration, the average value of the empirical FCD is Represents the bifurcation parameter At the t-th iteration, the average value of the simulated FCD is calculated; S800, constructing an approximate loss function; S900, implement the gradient descent process, and use the data feature correction value And approximate loss function, perform gradient correction of data features during gradient descent, and find parameters and to complete the parameters of the Landau-Stewart oscillator network and The fitting of 2. The method according to claim 1, wherein: In step S100, The simulation data is generated using a Landau-Stewart oscillator network and randomly generated connection structures and parameter values ​​of various simulated brain networks.

3. The method according to claim 1, wherein: The method further comprises the steps of: S1000: Further use various simulated brain networks and random parameters to conduct a series of tests to evaluate the universality and accuracy of the fit; when the universality and accuracy need to be improved, repeat S100 to S900 until the universality and accuracy meet the requirements.

4. The method according to claim 1, wherein: Step S100 includes: (1) Setting up the simulated brain network: The connection density of the simulated brain network is set according to the connection density of the real brain network; the connection weight is uniformly distributed The simulated brain network is randomly selected from the network and is an undirected network.

5. The method according to claim 4, wherein: When the simulated brain network is a symmetric undirected network with M nodes, the SC matrix is ​​an M×M symmetric matrix, in which each row and column represents a node in the network; each element SC in the SC matrix ij represents the connection weight from node i to node j.

6. The method according to claim 3, wherein: In order to ensure the diversity and unpredictability of the simulated brain network, two network models were selected during the testing process: a connection structure based on the Watts-Strogatz (WS) network model, and a connection structure based on the Barabási-Albert (BA) network model; The fitting parameters were confirmed by comparing the connection structures of the two network models with the connection structures based on the real brain network. and effectiveness; among them, For the Watts-Strogatz (WS) network model, to change the reconnection probability That is, the probability of cutting off and randomly reconnecting the edges of the network is introduced, and additional random elements are introduced to experiment with multiple probability values; For the Barabási-Albert (BA) network model, scale-free networks with different connection densities are constructed, in which the number of connections between nodes is unevenly distributed: some nodes have a large number of connections, while other nodes have fewer.

7. The method according to claim 5, wherein: Step S100 also includes: (2) Generate simulated data: A bifurcation parameter is randomly generated for each node in the simulated brain network As the benchmark bifurcation parameter value for evaluating the goodness of fit, the parameter set As a verification benchmark, it is input into the Landau-Stewart oscillator network, and the number of generated nodes is , the time length is The simulated data size is ; The generated simulation data is subjected to feature fitting and parameter estimation, wherein: the accuracy and robustness of the adopted parameter estimation method are evaluated by comparing the fitting results with the original parameters used to generate the simulation data, and the method is continuously improved; if the fitted parameters are highly consistent with the parameters used to generate the simulation data, the simulation data is trusted, wherein the trusted simulation data is not only used for the parameters of the Landau-Stewart oscillator network and It is suitable for clinical and other scenarios.

8. The method according to claim 1, wherein: Step S200 specifically includes: The following key data features are determined based on the simulation data, and the parameters of the Landau-Stewart oscillator network are adjusted based on these key data features: and The fitting is performed, and these multiple key data features are also used as quantitative measures of the goodness of data fitting; The multiple key data features are determined in the following manner: (1) Determine functional connectivity (FC): The Pearson correlation coefficient between the simulated data of each pair of nodes in the simulated brain network is calculated to form a symmetric matrix, in which, Each element in the matrix is ​​the correlation coefficient r between the node pair xy , calculated as follows: in, and Is a node and nodes in The value of the simulation data at time, and The nodes are and nodes The mean of the simulated data; is the time length of the simulated data; Using the correlation coefficient r xy represents the size of functional connectivity (FC); (2) Determine the functional connectivity matrix (dFC) at each moment: Using the Hilbert Transform Compute Node , p respectively at time The phase of the simulated data and , and calculate the phase difference between the node pair of node n and p at each time point, and then use the cosine function to calculate the phase coherence matrix Elements : ； in, The matrix size is , is the number of nodes, is the time length of the simulated data; (3) Determine dynamic functional connectivity (FCD): To compare the empirical and simulated FCDs, the Kolmogorov-Smirnov (KS) distance, which is used to quantify the maximum difference between the cumulative distribution functions of two samples, is used to measure the difference between the elements of the upper triangular part of the empirical FCD matrix calculated from the data and the simulated FCD matrix, where The range of Kolmogorov-Smirnov or KS distance is , the smaller the KS distance is, the better the fitting degree of the data features is. The empirical FCD matrix is ​​a matrix based on phase consistency; wherein, The empirical FCD matrix is ​​a A symmetric matrix, where is the length of the sliding window time window; The elements of the empirical FCD matrix is the bifurcation parameter Iteration Step and The functional connection matrix corresponding to the step and The upper triangle part and The cosine similarity between is expressed as follows: ； In this way, the empirical FCD matrix is ​​calculated from the empirical data; wherein the empirical data is obtained by generating simulation data in S100; (4) Calculate the probabilistic metastable substate (PMS): Step a. Feature vector extraction: From the bifurcation parameters The first feature vector ( ); The steps to extract the first eigenvector of the matrix are as follows: 1). Definition of matrix: Suppose there is a matrix , this is a The matrix represents samples and Features 2). Calculate the covariance matrix: First, calculate the matrix The covariance matrix of ;if Each column represents a feature and the row represents a sample, then the covariance matrix Calculated by: ,in yes The transposed matrix of 3). Calculate eigenvalues ​​and eigenvectors: Solve the covariance matrix The eigenvalues ​​and eigenvectors of satisfy the following equations: ,in is the eigenvalue, is the corresponding eigenvector; 4). Extract the first eigenvector: Sort the eigenvalues ​​from large to small, and the eigenvector corresponding to the largest eigenvalue is the first eigenvector ,in, Represents the direction with the largest directional variance in the data set; 5). Standardization: The feature vector Normalize it so that its modulus is 1, that is, its length is 1; Step b. Cluster analysis: Performing cluster analysis using the feature vectors extracted in step a to discover common and recurring connection patterns between different time points and different samples; wherein the cluster analysis uses a k-means clustering algorithm to cluster data points, i.e., dFC matrices at different time points, into multiple clusters based on the similarity of the feature vectors; Step c. Define the metastable substate: Considering that each cluster in the clustering result represents a unique metastable substate, the centroid of each cluster, i.e., the cluster center, is defined as the metastable substate; Step d. Probability and transition probability calculation: For each metastable substate, calculate its probability of appearing at all time points, and analyze the transition probability between different substates, that is, the probability of converting from one substate to another; Step e. Measure goodness of fit: The KL distance is used to compare the empirical data and the fitted data to obtain the probability of the metastable state, where Metastable substate The empirical probability of is the simulation probability on the same metastable substate; The KL distance between the empirical data and the fitted data is calculated to evaluate the fit of the model to the actually observed brain state probability distribution: 。 9. The method according to claim 1, wherein: Step S500 specifically includes: The frequency spectrum characteristics of each node extracted from the simulation data are selected as key indicators to characterize the numerical value of the bifurcation parameters of each node; wherein, This indicator is used to formulate the basis and evaluation criteria for parameter updating during the gradient descent process, so as to effectively guide the adjustment of bifurcation parameters during the gradient descent optimization process; This indicator is the node heterogeneity and is expressed as , where each element f of the indicator vector i is the spectrum feature of the i-th brain region; The process of calculating this indicator using simulated data includes: (1) Demean: The simulated data of each node is de-meaned so that the mean of the processed simulated data at all time points on each node is zero, which helps to remove the fixed deviation in the data and better analyze its fluctuations and changes; (2) Filtering processing: Filter the simulated data after removing the mean value of each node; (3) Detrend: Perform de-linear trend processing on the simulated data after de-averaging and filtering to eliminate possible linear changes in the data and ensure the accuracy of the analysis; (4) Fast Fourier Transform (FFT): Perform fast Fourier transform on the simulated data after de-linearization to convert it from the time domain to the frequency domain, thereby obtaining the spectrum characteristics of each frequency component; (5) Calculate the proportion of low-frequency spectrum features of each node in the simulated data: A low-frequency frequency band range where the blood oxygen level dependent (BOLD) signal time series is located is selected, and the value of the spectral characteristics that can characterize the simulated data within this range is calculated. The proportion of low-frequency spectral characteristics of each node of the simulated data is calculated by dividing the sum of the spectral characteristics of this frequency band by the sum of the spectral characteristics of the full frequency band range; Among them, the low-frequency frequency band range is generally concentrated between 0.01Hz and 0.1Hz. This frequency range is considered to be most relevant to neural activity in a resting state, mainly because the changes in blood oxygen caused by neural activity are relatively slow in time, and there is a large amount of physiological noise in the high-frequency band. Therefore, the BOLD signal in fMRI analysis mainly focuses on low-frequency fluctuations to effectively reflect the functional state of the brain and its network activity; Among them, the proportion f of any spectrum feature is calculated by the following formula: ， Among them, in digital signal processing, the full frequency range is half of the sampling frequency; according to the Nyquist theorem, in order to unambiguously capture and reconstruct all frequency components in the signal, the sampling must be performed at least at twice the rate of the highest frequency component in the signal; when TR=2s, the full frequency range is defined as 0.01-0.25Hz; Among them, TR stands for repetition time, which refers to the interval between two adjacent radio frequency pulses during magnetic resonance imaging; specifically, TR is the time interval from the start of one imaging signal acquisition to the start of the next imaging signal acquisition, usually expressed in seconds or milliseconds; (6) Calculate the proportion of spectrum features in different frequency bands for each node in the simulated data: In addition to the low-frequency band where the BOLD signal time series is located, an arbitrary frequency band is selected, and the ratio of its spectral characteristics to the spectral characteristics of the full frequency band is calculated as a way to characterize the node dynamics; wherein the ratio is also calculated according to the proportion f of any spectral feature described above; Furthermore, this ratio is processed by standardization, exponential or multiple scaling to obtain an index more suitable for specific analysis or application, so as to better characterize the size of each node bifurcation parameter; (7) Calculate the spectrum characteristics of the full frequency band of each node of the simulated data: The frequency center, bandwidth, spectrum flatness, spectrum roll-off, and spectrum entropy of the entire frequency band are calculated to characterize the size of the bifurcation parameters of each node.

10. The method according to claim 9, wherein: The low-frequency amplitude fraction (fALFF) is preferentially selected as the index for calculating the gradient in the gradient descent process to achieve parameter update and improve the goodness of fitting of data spectrum characteristics, so as to effectively guide the adjustment of bifurcation parameters in the gradient descent optimization process. fALFF focuses on the amplitude of the low-frequency signal, that is, the ratio of the amplitude of the signal in the low-frequency signal band where the BOLD signal time series is located to the amplitude of the full-band signal, which is the node heterogeneity and is expressed as .

Images