A method for calculating the balance between excitability and inhibition of epilepsy based on a hybrid dynamic causal model
The power spectral density function of epileptic EEG signals was reconstructed using the cPBM model and H-DCM algorithm, solving the problem of multi-peak fitting, improving the accuracy of parameter estimation, and quantifying the changes in the balance between excitability and inhibition during epileptic seizures.
Patent Information
- Application Number
- CN202310447881.0
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2023-04-24
- Publication Date
- 2025-11-07
- Estimated Expiration
- 2043-04-24
AI Technical Summary
Existing neuronal cluster models struggle to fit the spectrum of epileptic EEG signals with multiple peaks, and the spectral DCM algorithm is prone to local optima problems in parameter estimation.
The cPBM model was used to simulate epileptic EEG signals, and the H-DCM algorithm was used for parameter estimation. The local optimum problem was solved by a hybrid annealing scheme, and the model parameters were iteratively optimized by combining the EM algorithm.
The power spectral density function of EEG signals at different stages of epileptic seizures was effectively reconstructed, improving the accuracy of model parameter estimation and quantifying the EI balance changes during epileptic seizures.
Smart Images
Figure CN116439726B_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The application relates to a kind of epilepsy excitability and inhibitory balance calculation method based on hybrid dynamic causal model, belong to electroencephalogram signal processing technical field. BACKGROUND
[0002] Epilepsy is one of the most common neurological diseases, affecting about 650 million people worldwide. The characteristics of this disease are repeated seizures, and patients usually show abnormal feelings, actions or even loss of consciousness. Existing research shows that epilepsy is an acute, recurrent or paroxysmal neurological disease caused by excessive or synchronous discharges of brain neurons. According to different stages of seizures, seizures are usually divided into interictal, preictal and ictal. However, the biological mechanism of seizures is complex and unclear; but studies have shown that excessive or synchronous fast discharge activity of seizures is closely related to the balance between excitability and inhibition between neurons. On the other hand, the fast discharge activity related to seizures can be directly reflected in the EEG signal. And EEG signal has the advantages of high resolution, low acquisition cost, safety and convenience, so EEG signal is widely used in the related research of seizures.
[0003] In computational neuroscience, the electrical activity of the cerebral cortex region can be generated by a set of nonlinear differential equations, i.e. neuron cluster model. Studies have found that real EEG signals related to seizure activity always show some spectral characteristics, and the corresponding PSD usually has one or two peaks between alpha and beta bands (8-30 Hz). But the existing neuron cluster model is difficult to fit the spectrum with multiple peaks. In order to solve these problems, Professor Ursino proposed a completely physiological model (cPBM), which reconstructed the EEG signal with PSD having one or two peaks in the PSD of the cortical activity related to motor perception by adding new inhibitory self-loops and new inputs to fast inhibitory interneurons. The application simulates the PSD of epilepsy EEG signal based on cPBM model, in order to study the potential relationship between model parameters and seizures.
[0004] In recent years, there are more and more researches on parameter estimation in neuron cluster model. These methods can be divided into two categories: heuristic search algorithm and Bayesian estimation method. Heuristic search algorithm usually minimizes the objective function which describes the deviation between the real signal and the output of neuron cluster model through search strategy. However, heuristic search algorithm is suitable for the case that the number of parameters to be estimated in neuron cluster model is small, and when the number of parameters to be estimated is large, the time complexity of the algorithm is often large. Bayesian estimation method regards each estimated parameter in neuron cluster model as a random variable with a certain probability distribution, and the common ones include Kalman filter and dynamic causal modeling. Among them, Kalman filtering technique is more suitable for parameter estimation in time domain, and DCM (DCM in time domain) can be regarded as a variant of Kalman filtering. Therefore, in order to perform parameter estimation on cPBM model in frequency domain, the present application adopts spectral DCM algorithm. However, spectral DCM is prone to local optimum problem in the estimation process. Therefore, the present application introduces simulated annealing algorithm into spectral DCM algorithm, and adopts a hybrid annealing scheme in two directions (heating / cooling), that is, H-DCM algorithm.
[0005] In summary, the calculation method of epilepsy excitability and inhibitory balance based on hybrid dynamic causal model mainly uses cPBM model to simulate each state of epilepsy, and uses H-DCM algorithm to perform parameter estimation on cPBM model, so as to study the potential relationship between model parameters and epilepsy seizure. SUMMARY
[0006] The present application is to solve the above problems, and proposes a calculation method of epilepsy excitability and inhibitory balance based on hybrid dynamic causal model, which is used to study the potential relationship between E-I balance and epilepsy seizure.
[0007] The purpose of the present application is achieved by the following technical scheme: a calculation method of epilepsy excitability and inhibitory balance based on hybrid dynamic causal model, which comprises:
[0008] S1: establishing a neuron cluster model module,
[0009] The cPBM model is composed of pyramidal cells (P p ), excitatory neurons (P e ), slow inhibitory interneurons (P s ) and fast inhibitory interneurons (P f ). The four neuron clusters interact with each other through excitatory connections {C1, C2, C3, C5} and inhibitory connections {C4, C6, C7, C8} (as shown in Figure 2 ). and are two independent white noise external inputs following Gaussian probability distribution, which act on Pp and P f wherein y(t) represents the output of the model, the cPBM model adopts a set of continuous time differential equations to represent, as shown in the following formula (1):
[0010]
[0011] where x i represents the state variable of the system, S(x(t)) is a Sigmoid function, defined as:
[0012]
[0013] where, r=0.56mV -1 and v0=6mV respectively determine the steepness and translation position of the Sigmoid function shape, e0=6mV determines the maximum value of the Sigmoid function. The output y(t) of the model is:
[0014] y(t)=x3(t) (3) The prior values of other parameters in the cPBM model are shown in Table 1. Studies have shown that the balance between excitation and inhibition in the neuron cluster may change during the transition from the interictal period to the ictal period. Therefore, (wherein the superscript T represents the transpose operator) is the parameter information focused on in the present application.
[0015] Table 1 Parameters in the cPBM model:
[0016]
[0017] S2: Establish a power spectral density function calculation module
[0018] Since the epileptic electroencephalogram signal is a signal with obvious frequency domain characteristics, this paper simulates different seizure states of epilepsy from the perspective of frequency domain. In order to realize the frequency domain DCM algorithm, this paper converts the time domain signal into the frequency domain signal, that is, the power spectral density function. Considering different sources of EEG signals, the conversion process mainly includes the frequency spectrum conversion of the model output signal and the frequency domain conversion of the determined EEG signal.
[0019] S2.1 Predicting the power spectral density function
[0020] Step one: Establish the state space equation,
[0021]
[0022] The cPBM model represented by equations (1) and (3) can be described by a set of 14 first order differential equations and a state space equation for the output y(t) (see equation (4)) where f(x(t), θ) is a vector of nonlinear functions associated with the state vector and θ, representing the vector of two input white Gaussian noise. Except for D 5,1 = Aa / C2, D 13,2 = Aaand Q 1,3 = 1, the values of the matrices and the row vector are zero.
[0023] Step two: Linearization,
[0024]
[0025] Equation (4) can be expanded by the first order Taylor formula (see equation (5)) where is the Jacobian matrix.
[0026] Step three: Fourier transform,
[0027]
[0028] Equation (6) represents the Fourier transform of equation (5) where v represents the frequency variable, such that the transfer function
[0029] Step four: Obtain the predicted power spectral density function,
[0030] G(v, θ) = H(v, θ) G u (v, θ) H H (v, θ) (7)
[0031] According to equation (7), the predicted power spectral density function is defined by G(v, θ) and is computed from the corresponding u(t) and H(v, θ), respectively.
[0032] S2.2 Sampled power spectral density function
[0033] The sampled power spectral density function is computed based on a 12th order auto-regressive model and is defined as
[0034] Finally, for each frequency v, and all elements of G(v, θ) are stacked in the column vectors and g(θ), respectively.
[0035] S3: Hybrid simulated annealing principle module:
[0036] To investigate the state changes of the brain regions during the seizure, the spectral DCM algorithm was used to estimate the parameters in the cPBM model. The likelihood function of the model output is p(G), whose logarithmic expression is:
[0037] ln p(G) = F + KL(q(θ), p(θ|G)) (8)
[0038] where KL is the Kullback-Leible divergence, F represents the free energy, q(θ) is the probability density function of the parameter θ to be estimated, and p(θ|G) is the posterior probability distribution function of the parameter θ. For any given q(θ), KL(q(θ), p(θ|g)) ≥ 0, so ln p(G) ≥ F. Therefore, ln p(G) can be solved by maximizing F. The mathematical expression of F is:
[0039]
[0040] The maximization of F is solved by the EM algorithm.
[0041] For each model parameter vector θ, the predicted power spectral function of the cPBM output signal is calculated and compared with the sample power spectral function calculated from the EEG signal of the epilepsy patient. The parameter estimation process aims to maximize the free energy by iteration with the EM algorithm. However, the local optimal problem may occur in the traditional spectral DCM. Therefore, a hybrid deterministic annealing DCM (H-DCM) is proposed in the spectral DCM as follows:
[0042]
[0043] where 1 / β is the temperature parameter. The traditional spectral DCM can be regarded as a special case, i.e., β = 1. When the temperature is high (relatively small) (for example, β is small (relatively high)), the shape of the objective function F is smooth (relatively steep), and by changing the smoothness of the objective function, a relatively good initial value can be obtained. Therefore, a hybrid simulated annealing algorithm is introduced in this paper to solve the local optimal problem, as shown in Figure 3 The H-DCM algorithm adjusts the temperature in two directions (β changes from 1.0 to 1.6 (β new = β + 0.2)) to obtain a relatively good initialization; then, the temperature gradually increases, β changes from 1.6 to 1.0 (β new = β - 0.2)) to find a better estimate after each maximization.
[0044] S4: Excitability and inhibitory balance calculation:
[0045] The data of 10 epilepsy patients in Hauz Khas, New Delhi, India, are estimated by the method described in S1-S3, and the parameter changes of the interictal, preictal and ictal stages are counted respectively, and the best index for quantifying the balance between excitability and inhibition is calculated.
[0046] Compared with the prior art, the present application has the following advantages:
[0047] (1) The epilepsy EEG signal modeling method based on the cPBM model is adopted in the present application. Studies have shown that the power spectral density function of the EEG signal of the epilepsy patient always has one or more peaks, and the peak is located between the alpha-beta wave band (8-30 Hz). The power spectral density function fitting with multiple peaks has always been a difficulty. In order to solve this problem, the cPBM model is introduced to reconstruct the EEG signal and its power spectral density function in different stages of the epilepsy seizure, and a better reconstruction effect is obtained.
[0048] (2) The present application proposes a parameter estimation method based on a hybrid dynamic causal model. Since the cPBM model contains many physiological parameters, its solving process is a multi-parameter solving problem, and in the process of solving the model parameters, the local optimal problem exists in the spectral dynamic causal model. In order to improve the accuracy of estimating the model parameters in the cPBM using the spectral dynamic causal model, a hybrid dynamic causal model algorithm is proposed in this paper, which introduces a temperature coefficient β in the objective function, and uses a hybrid annealing scheme containing heating and cooling, thereby effectively enhancing the robustness of accurately estimating the model parameters.
[0049] (3) The present application proposes a quantification method of E-I balance in the process of epilepsy seizure. In order to study the transition relationship between epilepsy seizure and E-I balance, the changes of each model parameter in the cPBM model under different epilepsy states are analyzed in detail, and the E-I balance is quantified from the perspective of the interaction between neuron clusters. BRIEF DESCRIPTION OF DRAWINGS
[0050] Figure 1 The algorithm flowchart of the present application is shown in the figure;
[0051] Figure 2 The schematic diagram of the cPBM model is shown in the figure;
[0052] Figure 3 The annealing scheme of the H-DCM algorithm is shown in the figure;
[0053] Figure 4 The fitting effect of the reconstructed signal and the real signal of the cPBM model is shown in the figure;
[0054] Figure 5 The E-I balance calculation result is shown in the figure. DETAILED DESCRIPTION
[0055] In order to deepen the understanding and understanding of the present application, the present application is further illustrated below in conjunction with the drawings and specific embodiments.
[0056] Example 1:
[0057] Reference Figure 1 — Figure 5 A method for calculating the balance of epilepsy excitability and inhibition based on a hybrid dynamic causal model, the method comprising:
[0058] S1: establishing a neuron cluster model module,
[0059] The cPBM model is composed of pyramidal cells (P p ), excitatory neurons (P e ), slow inhibitory interneurons (P s ) and fast inhibitory interneurons (P f ). The four neuron clusters interact with each other through excitatory connections {C1, C2, C3, C5} and inhibitory connections {C4, C6, C7, C8} (as shown in Figure 2 ). and are two independent white noise external inputs following Gaussian probability distribution, acting on P p and P f respectively, where y(t) represents the output of the model, and the cPBM model adopts a set of continuous time differential equations to represent, as shown in the following formula (1):
[0060]
[0061] where x i represents the state variable of the system, and S(x(t)) is a Sigmoid function, defined as:
[0062]
[0063] where, r=0.56mV -1 and v0=6mV determine the steepness and translation position of the Sigmoid function shape, and e0=6mV determines the maximum value of the Sigmoid function. The output y(t) of the model is:
[0064]
[0065] The prior values of other parameters in cPBM are shown in Table 1. Studies have shown that the balance between excitation and inhibition in neuronal clusters can change during the transition from the interictal to the ictal phase. Therefore, where the superscript T denotes the transpose operator, is the parameter information of focus in this invention.
[0066] Table 1 Parameters in cPBM model:
[0067]
[0068] S2: Establish a power spectral density function calculation module,
[0069] Since the epileptic electroencephalogram signal is a signal with obvious frequency domain characteristics, this paper simulates different seizure states of epilepsy from the perspective of frequency domain. In order to realize the frequency domain DCM algorithm, this paper converts the time domain signal into the frequency domain signal, that is, the power spectral density function. Considering different sources of EEG signals, the conversion process mainly includes the frequency spectrum conversion of the model output signal and the frequency domain conversion of the determined EEG signal.
[0070] S2.1 Predicting the power spectral density function
[0071] Step one: Establish the state space equation,
[0072]
[0073] The cPBM model represented by formula (1) and formula (3) can be described by a set of 14 first-order differential equations and a state space equation of component output y(t) (see equation (4)), where f(x(t), θ) is a nonlinear function vector associated with the state vector and θ, representing the vector of two input white Gaussian noise. In addition to D 5,1 = Aa / C2, D 13,2 = Aaand Q 1,3 = 1, the values of matrices and row vectors are all zero.
[0074] Step two: Linearization,
[0075]
[0076] Equation (4) can be expanded by the first-order Taylor formula (see equation (5)), where is the Jacobian matrix.
[0077] Step three: Fourier transform,
[0078]
[0079] Equation (6) represents the Fourier transform of equation (5), where v represents the frequency variable, such that the transfer function is calculated
[0080] Step four: obtain the predicted power spectral density function,
[0081]
[0082] According to equation (7), the predicted power spectral density function is defined by G(v,0) and is calculated from the corresponding u(t) and H(v,0), respectively.
[0083] S2.2 Sampled power spectral density function
[0084] The sampled power spectral density function is calculated based on a 12-stage autoregressive model and is defined as
[0085] Finally, for each frequency v, and all elements of G(v,0) are stacked in column vectors and g(0), respectively.
[0086] S3: Hybrid simulated annealing principle module:
[0087] To investigate the changes in the state of brain regions during a seizure, this paper uses the spectral DCM algorithm to estimate the parameters in the cPBM model. The likelihood function of the model output is p(G), and its logarithmic expression is:
[0088] ln p(G) = F + KL(q(0), p(0|G)) (8)
[0089] where KL is the Kullback-Leible divergence, F represents the free energy, q(0) is the probability density function of the parameter 0 to be estimated, and p(0|G) is the posterior probability distribution function of the parameter 0. For any given q(0), KL(q(0), p(0|g)) ≥ 0, so ln p(G) ≥ F. Therefore, ln p(G) can be solved by maximizing F. The mathematical expression of F is:
[0090]
[0091] Maximizing F is solved by the EM algorithm.
[0092] For each model parameter vector θ, the predicted power spectral function of the cPBM output signal is computed and compared with the sampled power spectral function computed from the EEG signal of the epileptic patient. The parameter estimation procedure aims to maximize the free energy by iteration using the EM algorithm. However, the local optimum problem that can occur in the traditional spectral DCM, a hybrid deterministic annealing DCM (H-DCM) is proposed in the spectral DCM as follows:
[0093]
[0094] where 1 / β is the temperature parameter. The traditional spectral DCM can be regarded as a special case, i.e. β = 1. When the temperature is high (relatively small) (e.g. β is small (relatively high)), the shape of the objective function F is smooth (relatively steep), and a better initial value can be obtained by changing the smoothness of the objective function. Therefore, a hybrid simulated annealing algorithm is introduced in this paper to solve the local optimal problem, as shown in Figure 3 The H-DCM algorithm adjusts the temperature in two directions (β varies from 1.0 to 1.6 (β new = β + 0.2) to obtain a relatively good initialization; then, the temperature gradually increases, β varies from 1.6 to 1.0 (β new = β - 0.2) to find a better estimate after each maximization.
[0095] S4: Excitability and inhibitory balance calculation:
[0096] The parameter estimation of the data of 10 epileptic patients in Hauz Khas Neuro and Sleep Center, New Delhi, India is carried out by the method described in S1-S3 above, and the parameter changes of the interictal, preictal and ictal periods are counted respectively, and the best index capable of quantifying the excitability and inhibitory balance is calculated.
[0097] Example 2: The data set used in the present application adopts the electroencephalogram data set of epilepsy in Hauz Khas Neuro and Sleep Center, New Delhi, India. The data set collects typical segmented electroencephalogram time series records of 10 epileptic patients from Hauz Khas Neuro and Sleep Center, New Delhi, India. During the collection process, according to the 10-20 electrode placement system, gold-plated scalp electroencephalogram electrodes are placed, and the signal is collected at a sampling frequency of 200 Hz. The collected signal is filtered between 0.5 and 70 Hz, and then divided into interictal, preictal and ictal periods. Each downloadable folder contains 50 EEG time series signal MAT files, and each MAT file consists of 1024 EEG time series data samples with a duration of 5.12 seconds.
[0098] Figure 4The simulation results of one of the data of the three stages of epileptic seizure ((a), (b), (c) and (d) represent interictal, preictal (single peak), preictal (double peak) and ictal, respectively) are shown, wherein Figure 4 (a-i), (b-i), (c-i) and (d-i) are the real EEG signals (Sample EEG, represented by solid lines) and the EEG signals (Reconstructed EEG) reconstructed by solving equation (1) with the estimated values using the 4th order Runge-Kutta method at a sampling frequency of 200 Hz; Figure 4 (a-ii), (b-ii), (c-ii) and (d-ii) show the Sample PSDs of the three stages of epileptic seizure The estimated results obtained using the H-DCM The calculated Predicted PSDs The Reconstructed PSDs The fitting effects between them. In order to measure the fitting effect of the power spectral density function, the RMSE is used to measure the fitting degree before the PSD in this paper, wherein R p represents the fitting degree between the Sample PSD and the Predicted PSD, and R r represents the fitting degree between the Sample PSD and the Reconstructed PSD. As can be seen from the figure, the four real data can be reconstructed by the cPBM model and the H-DCM algorithm.
[0099] Secondly, in order to further verify the fitting effect of the cPBM on the real data, the RMSE and the free energy of all real data are calculated, and the evaluation value is calculated to evaluate the simulation effect of the cPBM model, and the statistical results are shown in Table 2. As can be seen from Table 2, whether the fitting degree between the Predicted PSD and the Sample PSD or the fitting degree between the Reconstructed PSD and the Sample PSD both remain at a low value. Therefore, the cPBM model and the H-DCM algorithm can simulate different states of epileptic seizure.
[0100] Table 2 Fitting effect of H-DCM algorithm based on cPBM model in different stages of epileptic seizure
[0101]
[0102] To investigate the dynamic transition mechanism of epileptic seizures, Table 3 presents the parameter estimation results of 50 EEG signals for each of the different seizure stages calculated using the H-DCM algorithm proposed in this invention. To measure the magnitude of change in model parameters during epileptic seizures, this paper introduces the rate of change I to measure the magnitude of parameter change, which is defined as:
[0103]
[0104] During the dynamic changes from the interictal period to the ictal period, this paper observed a significant increase in the parameter set {C1,C2,C5,C7,A,G,g}, and a significant increase in the parameter set {C8,b} (|I θ (>9%). Regarding internal connectivity, the results indicate that the interictal-to-interictal transition in the cPBM can be facilitated by increased connectivity strength between pyramidal cells and excitatory interneurons. Increased connection strength between pyramidal cells and fast inhibitory interneurons ( and And the reduction in the strength of self-circuit connections of rapidly inhibiting interneurons. This can be explained by one approach. On the other hand, clinical studies have shown that a decrease in inhibition of rapidly inhibitory neurons can explain the imbalance between excitability and inhibition, leading to abnormal brain discharges and causing the patient's transition from interictal to ictal phases.
[0105] Table 3 shows the average estimated parameters of 50 data points from three epileptic stages using the H-DCM algorithm. (mean ± variance)
[0106]
[0107] To measure EI balance by the ratio of excitatory to inhibitory connections between internal neurons, we define the EI balance index: E PF =C5 / C8 (meaning from P) p To P f Excitation connection with P f (The ratio between self-inhibitory connections). Figure 5 This chart displays box plots of three indicators for 50 signals across three phases. From this chart, we can see the changes in the three indicators of EI balance during the pre-onset phase, E... PF It can effectively quantify the three stages of epilepsy, with mean E values during the interictal, preictal, and ictal periods. PF The values were 0.42, 0.78, and 1.23, respectively, with the rate of change increasing by approximately 190% between the interictal and ictal periods.
[0108] The above shows and describes the basic principles, main features and advantages of the present application. Those skilled in the art should understand that the present application is not limited to the above-mentioned embodiments, and the above-mentioned embodiments and descriptions in the specification are only preferred examples of the application and are not intended to limit the application. Without departing from the novel spirit and scope of the present application, various changes and improvements can be made to the application, and these changes and improvements all fall within the scope of the claimed application. The scope of protection of the present application is defined by the appended claims and their equivalents.
Claims
1. A method for calculating the balance between excitability and inhibition of epilepsy based on a hybrid dynamic causal model, characterized in that, The method comprises the following steps: S1: establishing a neuron cluster model module, establishing a neuron cluster model based on cPBM, solving the cPBM equation set by the Runge-Kutta method to obtain a time domain signal, and simulating different stages of a seizure, S2: establishing a power spectral density function calculation module, for converting the time domain signal into a frequency domain signal, mainly including calculating the sampling power spectral density function of the real EEG signal based on an autoregressive model; and calculating the predicted power spectral density function according to the differential equation set of the cPBM model, S3: establishing a hybrid simulated annealing principle module, for solving the model parameters of the cPBM model, mainly including establishing an objective function, and using an EM algorithm to estimate the model parameters to maximize the free energy, and then fitting the sampling power spectral density function and the predicted power spectral density function, S4: calculation of excitatory and inhibitory balance, By further analyzing the estimation results of the model parameters, the excitatory and inhibitory balance index E is calculated using C5 and C8 in the estimation results of the model parameters pf ; The step S1 is specifically: The cPBM model consists of pyramidal cells P p Excitatory neurons P e Slow inhibitory interneurons P s and rapid inhibitory interneurons P f These four neuronal clusters interact through excitatory and inhibitory connections, respectively. and There are two independent white noise external inputs that follow a Gaussian probability distribution, acting on P respectively. p and P f ,in y(t) represents the output of the model. The cPBM model is represented by a set of continuous time differential equations, as shown in the following formula (1): where A represents the average gain of excitatory synapses, B represents the average gain of slow inhibitory synapses, G represents the average gain of fast inhibitory synapses, {C1, C2, C3, C5} are excitatory connection parameters, {C4, C6, C7, C8} are inhibitory connection parameters, x(t) = [x1(t),..., xN(t)]Trepresents the state variable of the system, S(x(t)) is a Sigmoid function defined as: 14 (t)] T , where A represents the average gain of excitatory synapses, B represents the average gain of slow inhibitory synapses, G represents the average gain of fast inhibitory synapses, {C1, C2, C3, C5} are excitatory connection parameters, {C4, C6, C7, C8} are inhibitory connection parameters, x(t) = [x1(t),..., xN(t)]Trepresents the state variable of the system, S(x(t)) is a Sigmoid function defined as: where r = 0.56 mV -1 and v0= 6 mV determine the steepness and the shift of the sigmoid function, respectively, and e0= 6 mV determines the maximum of the sigmoid function, and the output of the model y(t) is: y(t)=x3(t) (3) In the transition process from the interictal period to the ictal period, the balance between excitation and inhibition in the neuron cluster may change, where the superscript T denotes the transpose operator; The step S3 is specifically: The spectral DCM algorithm is used to estimate the parameters in the cPBM model, and the likelihood function of the model output is p(P), and the logarithmic expression is: ln p(P)=F+KL(q(θ),p(θ|P)) (8) Where KL is the Kullback-Leible divergence, F represents the free energy, q(θ) is the probability density function of the parameter θ to be estimated, and p(θ|P) is the posterior probability distribution function of the parameter θ. For any given q(θ), KL(q(θ), p(θ|P))≥0, so ln p(P)≥F, so ln p(P) is solved by maximizing F, and the mathematical expression of F is: Maximizing F is solved by the EM algorithm, For each model parameter vector θ, the predicted power spectral function of the cPBM output signal is calculated and compared with the sampling power spectral function calculated from the EEG signal of the patient with epilepsy. The parameter estimation process aims to maximize the free energy through EM algorithm iteration. However, the local optimum problem may occur in the traditional spectral DCM. In the spectral DCM, a hybrid deterministic annealing DCM, namely H-DCM, is proposed, which is defined as follows: Where 1 / β is the temperature parameter, and the traditional spectral DCM is regarded as a special case, that is, β=1. When the temperature is high, the shape of the objective function F is smooth. By changing the smoothness of the objective function, a better initial value is obtained. The H-DCM algorithm adjusts the temperature in both directions, β varying from 1.0 to 1.6, β new = β + 0.2, and then gradually increasing the temperature, β varying from 1.6 to 1.0, β new = β - 0.2, to find better estimates after each maximization.