CN116439726A - 一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法 - Google Patents
一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法 Download PDFInfo
- Publication number
- CN116439726A CN116439726A CN202310447881.0A CN202310447881A CN116439726A CN 116439726 A CN116439726 A CN 116439726A CN 202310447881 A CN202310447881 A CN 202310447881A CN 116439726 A CN116439726 A CN 116439726A
- Authority
- CN
- China
- Prior art keywords
- model
- power spectral
- density function
- function
- cpbm
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 230000001037 epileptic effect Effects 0.000 title claims abstract description 28
- 230000002401 inhibitory effect Effects 0.000 title claims abstract description 26
- 206010001497 Agitation Diseases 0.000 title claims abstract description 21
- 238000004364 calculation method Methods 0.000 title claims abstract description 18
- 206010010904 Convulsion Diseases 0.000 claims abstract description 56
- 230000003595 spectral effect Effects 0.000 claims abstract description 47
- 238000000034 method Methods 0.000 claims abstract description 31
- 210000002569 neuron Anatomy 0.000 claims abstract description 23
- 230000001364 causal effect Effects 0.000 claims abstract description 15
- 206010015037 epilepsy Diseases 0.000 claims abstract description 12
- 238000001228 spectrum Methods 0.000 claims abstract description 10
- 208000028329 epileptic seizure Diseases 0.000 claims abstract description 9
- 238000005070 sampling Methods 0.000 claims abstract description 9
- 238000002922 simulated annealing Methods 0.000 claims abstract description 9
- 238000000137 annealing Methods 0.000 claims abstract description 7
- 238000004138 cluster model Methods 0.000 claims abstract 3
- 230000006870 function Effects 0.000 claims description 64
- 238000004422 calculation algorithm Methods 0.000 claims description 32
- 239000013598 vector Substances 0.000 claims description 15
- 230000008859 change Effects 0.000 claims description 10
- 230000008569 process Effects 0.000 claims description 10
- 230000002964 excitative effect Effects 0.000 claims description 9
- 210000001926 inhibitory interneuron Anatomy 0.000 claims description 7
- 230000005764 inhibitory process Effects 0.000 claims description 5
- 210000002763 pyramidal cell Anatomy 0.000 claims description 5
- 230000007704 transition Effects 0.000 claims description 5
- 230000003993 interaction Effects 0.000 claims description 4
- 238000005315 distribution function Methods 0.000 claims description 3
- 230000005284 excitation Effects 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 238000012546 transfer Methods 0.000 claims description 3
- 238000001816 cooling Methods 0.000 abstract description 3
- 238000010438 heat treatment Methods 0.000 abstract description 3
- 206010012812 Diffuse cutaneous mastocytosis Diseases 0.000 description 16
- 230000000694 effects Effects 0.000 description 10
- 238000006243 chemical reaction Methods 0.000 description 6
- 210000004556 brain Anatomy 0.000 description 4
- 210000005036 nerve Anatomy 0.000 description 4
- 230000001537 neural effect Effects 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 238000010845 search algorithm Methods 0.000 description 3
- 208000012902 Nervous system disease Diseases 0.000 description 2
- 238000001541 differential confocal microscopy Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 230000000306 recurrent effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 239000000758 substrate Substances 0.000 description 2
- 206010000117 Abnormal behaviour Diseases 0.000 description 1
- 208000025966 Neurological disease Diseases 0.000 description 1
- 208000003443 Unconsciousness Diseases 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000001154 acute effect Effects 0.000 description 1
- 230000002902 bimodal effect Effects 0.000 description 1
- 230000007321 biological mechanism Effects 0.000 description 1
- 210000003710 cerebral cortex Anatomy 0.000 description 1
- 230000001054 cortical effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000010304 firing Methods 0.000 description 1
- 208000013403 hyperactivity Diseases 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 208000035824 paresthesia Diseases 0.000 description 1
- 230000001314 paroxysmal effect Effects 0.000 description 1
- 230000008447 perception Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 210000004761 scalp Anatomy 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
- A61B5/369—Electroencephalography [EEG]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
- A61B5/369—Electroencephalography [EEG]
- A61B5/372—Analysis of electroencephalograms
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
- A61B5/369—Electroencephalography [EEG]
- A61B5/372—Analysis of electroencephalograms
- A61B5/374—Detecting the frequency distribution of signals, e.g. detecting delta, theta, alpha, beta or gamma waves
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/40—Detecting, measuring or recording for evaluating the nervous system
- A61B5/4076—Diagnosing or monitoring particular conditions of the nervous system
- A61B5/4094—Diagnosing or monitoring seizure diseases, e.g. epilepsy
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Veterinary Medicine (AREA)
- Biophysics (AREA)
- Pathology (AREA)
- Engineering & Computer Science (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Public Health (AREA)
- Physics & Mathematics (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Psychiatry (AREA)
- Psychology (AREA)
- Neurology (AREA)
- Neurosurgery (AREA)
- Physiology (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明公开了一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法，主要包括建立神经元集群模型模块、建立功率谱密度函数计算模块、建立混合模拟退火原理模块以及兴奋性和抑制性平衡计算四个部分，神经元集群模型模块使用cPBM模型仿真癫痫发作各阶段EEG信号的功率谱密度函数；功率谱密度函数计算模块依据状态空间方程生成预测功率谱密度函数以及计算真实EEG信号的采样功率谱密度函数；混合模拟退火原理模块在动态因果模型中引入了模拟退火算法，提出了一种包括升温和降温的混合退火方案，用于提高模型参数估计的准确性；兴奋性和抑制性平衡计算使用模型参数估计结果中的C5和C8计算得到Epf，Epf从癫痫发作间期到癫痫发作期的增幅约达190％。
Description
技术领域
本发明涉及一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法，属于脑电信号处理技术领域。
背景技术
癫痫是最常见的神经疾病之一，影响着全世界约6500万人。这种疾病的特点是反复发作，患者通常表现为感觉、行动异常甚至失去意识。现有研究表明，癫痫是一种由大脑神经元过度或同步放电引起的急性、反复发作或阵发性的神经系统疾病，按照癫痫发作的不同阶段，癫痫发作通常分为发作间期、发作前和发作期。然而癫痫发作的生物学机制复杂、尚不清楚；但有研究指出，癫痫发作的过度或同步快速放电活动与神经元之间的兴奋性和抑制性平衡密切相关。另一方面，与癫痫发作相关的快速放电活动可以直观地反映在EEG信号种。且EEG信号具有高分辨率、采集成本低、安全便捷等优点，因此EEG信号并广泛用于癫痫发作的相关研究中。
在计算神经科学中，大脑皮层区域的电活动可以通过一组非线性微分方程生成，即神经元集群模型。有研究发现，与癫痫发作活动相关的真实EEG信号总是显示一些光谱特征，其对应的PSD在α和β频带(8～30Hz)之间通常具有一个或两个峰值。但现有的神经元集群模型很难拟合具有多个峰值的频谱。为了解决这些问题，Ursino教授提出了完全基于生理的模型(cPBM)，其通过向快速抑制中间神经元添加新的抑制性自环路和新的输入来重建具有与运动感知相关的皮层活动中具有一个或两个峰值的PSD的EEG信号。本发明基于cPBM模型对癫痫EEG信号地PSD进行仿真，以此来研究模型参数与癫痫发作之间的潜在关系。
近年来，关于神经元集群模型中参数估计的研究越来越多。这些方法可分为两类：启发式搜索算法和贝叶斯估计方法。启发式搜索算法通常通过搜索策略最小化描述真实信号与神经元集群模型输出之间偏差的目标函数。然而，启发式搜索算法适用于神经元集群模型中待估计参数较少的情况，当待估计参数较多时，算法的时间复杂度往往较大。贝叶斯估计算法将神经元集群模型中的每个估计参数视为具有一定概率分布的随机变量，常见的有包括卡尔曼滤波器和动态因果建模。其中，卡尔曼滤波技术更适合于在时域上进行参数估计，并且DCM(时域中的DCM)可以被视为卡尔曼滤波的一种变体。因此，为了在频域上对cPBM模型进行参数估计，本发明采用了频谱DCM算法。但频谱DCM在估计过程中容易受到局部最优问题。因此本发明在频谱DCM算法中引入了模拟退火算法，并采用一种两个方向(升温/降温)上的混合退火方案，即H-DCM算法。
综上，基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法主要采用cPBM模型对癫痫各个状态进行仿真，并使用H-DCM算法对cPBM模型进行参数估计，以此来研究模型参数与癫痫发作之间的潜在关系。
发明内容
本发明为解决上述问题，提出了一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法，用于研究E-I平衡与癫痫发作的潜在关系。
本发明的目的是通过以下技术方案来实现的：一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法，该方法包括：
S1：建立神经元集群模型模块，
cPBM模型由锥体细胞(Pp)、兴奋性神经元(Pe)、缓慢抑制性中间神经元(Ps)和快速抑制性中间神经元(Pf)构成。这四个神经元集群分别通过兴奋性连接{C1,C2,C3,C5}和抑制性连接{C4,C6,C7,C8}相互作用(如图2所示)。和/>是两个遵循高斯概率分布的独立白噪声的外部输入，分别作用于Pp和Pf，其中 y(t)表示模型的输出，cPBM模型采用一组连续的时间微分方程组来表示，如下公式(1)所示：
其中xi表示系统的状态变量，S(x(t))为Sigmoid函数，定义为：
其中，r＝0.56mV-1和v0＝6mV分别决定Sigmoid函数形状的陡峭程度和平移位置，e0＝6mV决定Sigmoid函数的最大值。模型的输出y(t)为：
y(t)＝x3(t) (3)
cPBM中其他参数的先验值如表1所示。研究表明，在从发作间期到发作期的过渡过程中，神经元集群中兴奋和抑制之间的平衡可能会发生变化。因此， (其中上标T表示转置运算符)为本发明中重点研究的参数信息。
表1cPBM模型中的参数：
S2：建立功率谱密度函数计算模块
由于癫痫脑电信号是一种频域特征明显的信号，因此本文从频域的角度来模拟癫痫的不同发作状态。为了实现频域DCM算法，本文将时域信号转换为频域信号，即功率谱密度函数。考虑到不同来源的EEG信号，转换过程主要包括对模型输出信号的频谱转换以及对确定EEG信号的频域转换。
S2.1预测功率谱密度函数
步骤一：建立状态空间方程，
公式(1)和公式(3)表示的cPBM模型可以用一组14个一阶微分方程和一个分量输出y(t)的状态空间方程(参见方程(4))来描述，其中f(x(t),θ)是与状态向量 和θ相关联的非线性函数向量，代表两个输入白高斯噪声的向量。除D5,1＝Aa/C2,D13,2＝Aa和Q1,3＝1外，矩阵/>和行向量/>的值皆为零。
步骤二：线性化，
公式(4)可以通过一阶泰勒公式展开(参见公式(5)),其中为雅可比矩阵。
步骤三：傅里叶变换，
公式(6)表示公式(5)的傅里叶变换，其中v表示频率变量，从而计算获得传递函数
步骤四：获得预测功率谱密度函数，
G(v,θ)＝H(v,θ)Gu(v,θ)HH(v,θ) (7)
根据公式(7)，预测功率谱密度函数由G(v,θ)定义，分别由对应的u(t)和H(v,θ)计算得出。
S2.2采样功率谱密度函数
采样功率谱密度函数基于12阶段自回归模型计算，并定义为
最后，对于每个频率v，和G(v,θ)的所有元素分别堆叠在列向量/>和g(θ)中。
S3：混合模拟退火原理模块：
为了研究癫痫发作时脑区内部在癫痫发作过程中的状态变化，本文使用频谱DCM算法来估计cPBM模型中的参数。模型输出的似然函数为p(G)，其对数表达式为：
ln p(G)＝F+KL(q(θ),p(θ|G)) (8)
其中KL为Kullback-Leible散度，F表示自由能，q(θ)为待估计参数θ的概率密度函数，p(θ|G)是参数θ的后验概率分布函数。对于任何给定的q(θ)，KL(q(θ),p(θ|g))≥0，因此lnp(G)≥F。所以lnp(G)可以通过最大化F来求解。F的数学表达式为：
最大化F通过EM算法求解。
对于每个模型参数向量θ，计算cPBM输出信号的预测功率谱函数，并将其与癫痫患者的EEG信号计算的采样功率谱函数进行比较。参数估计过程旨在通过EM算法迭代以最大化自由能。然而传统的频谱DCM可能出现的局部最优问题，本文在频谱DCM中提出了一种混合确定性退火DCM(H-DCM)定义如下:
其中，1/β是温度参数。传统的频谱DCM可以看作是一种特殊情况，即β＝1。当温度较高(相对较小)时(例如，β较小(相对较高))，目标函数F的形状是平滑的(相对陡峭)，通过改变目标函数的平滑程度，可以获得一个较好的初始值。因此，本发明本文引入了一种混合的模拟退火算法来解决局部的最优的问题，如图3所示。H-DCM算法通过在两个方向上调节温度(β从1.0变化到1.6(βnew＝β+0.2)，以获得相对良好的初始化；然后，温度逐渐升高，β从1.6变化到1.0(βnew＝β-0.2))以在每次最大化后寻找更好的估计。
S4：兴奋性和抑制性平衡计算：
通过上述S1-S3所述方法对印度新德里HauzKhas神经与睡眠中心的10名癫痫患者数据进行参数估计，并分别统计了癫痫发作间期、发作前期和发作期的参数变化情况，计算得到了能够对兴奋性和抑制性平衡进行量化的最佳指标。
相对于现有技术，本发明的优点如下：
(1)本发明采用了基于cPBM模型的癫痫EEG信号建模方法。研究表明，癫痫患者的EEG信号的功率谱密度函数总存在一个或多个峰值，且峰值位于α～β波段之间(8～30Hz)。而具有多个峰值的功率谱密度函数拟合一直是其难点。为了解决该问题，本发明引入了cPBM模型来重建癫痫发作不同阶段的EEG信号及其功率谱密度函数，并取得了较优的重建效果。
(2)本发明提出了一种基于混合动态因果模型的参数估计方法。由于cPBM模型包含了许多生理参数，其求解过程是一个多参数求解问题，而在求解模型参数过程中，频谱动态因果模型存在局部最优问题，为了提高使用频谱动态因果模型估计cPBM中模型参数的准确性，本文提出了一种混合动态因果模型算法，该算法在目标函数中引入了温度系数β，并使用一种包含升温和降温的混合退火方案，从而有效增强了准确估计模型参数的鲁棒性。
(3)本发明提出了癫痫发作过程中E-I平衡的量化方法。为了研究癫痫发作与E-I平衡之间的过渡关系，本文对cPBM模型中各个模型参数在不同癫痫状态下的变化进行了详尽分析，从神经元集群之间相互作用的角度对E-I平衡进行了量化。
附图说明
图1为本发明的算法流程图；
图2为cPBM模型的示意图；
图3为H-DCM算法的退火方案；
图4为cPBM模型的重构信号与真实信号的拟合效果；
图5为E-I balance计算结果。
具体实施方式
为了加深对本发明的认识和理解，下面结合附图和具体实施方式，进一步阐明本发明。
实施例1：
参见图1—图5，一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法，该方法包括：
S1：建立神经元集群模型模块，
cPBM模型由锥体细胞(Pp)、兴奋性神经元(Pe)、缓慢抑制性中间神经元(Ps)和快速抑制性中间神经元(Pf)构成。这四个神经元集群分别通过兴奋性连接{C1,C2,C3,C5}和抑制性连接{C4,C6,C7,C8}相互作用(如图2所示)。和/>是两个遵循高斯概率分布的独立白噪声的外部输入，分别作用于Pp和Pf，其中 y(t)表示模型的输出，cPBM模型采用一组连续的时间微分方程组来表示，如下公式(1)所示：
其中xi表示系统的状态变量，S(x(t))为Sigmoid函数，定义为：
其中，r＝0.56mV-1和v0＝6mV分别决定Sigmoid函数形状的陡峭程度和平移位置，e0＝6mV决定Sigmoid函数的最大值。模型的输出y(t)为：
cPBM中其他参数的先验值如表1所示。研究表明，在从发作间期到发作期的过渡过程中，神经元集群中兴奋和抑制之间的平衡可能会发生变化。因此， (其中上标T表示转置运算符)为本发明中重点研究的参数信息。
表1cPBM模型中的参数：
S2：建立功率谱密度函数计算模块，
由于癫痫脑电信号是一种频域特征明显的信号，因此本文从频域的角度来模拟癫痫的不同发作状态。为了实现频域DCM算法，本文将时域信号转换为频域信号，即功率谱密度函数。考虑到不同来源的EEG信号，转换过程主要包括对模型输出信号的频谱转换以及对确定EEG信号的频域转换。
S2.1预测功率谱密度函数
步骤一：建立状态空间方程，
公式(1)和公式(3)表示的cPBM模型可以用一组14个一阶微分方程和一个分量输出y(t)的状态空间方程(参见方程(4))来描述，其中f(x(t),θ)是与状态向量 和θ相关联的非线性函数向量，代表两个输入白高斯噪声的向量。除D5,1＝Aa/C2,D13,2＝Aa和Q1,3＝1外，矩阵/>和行向量/>的值皆为零。
步骤二：线性化，
公式(4)可以通过一阶泰勒公式展开(参见公式(5)),其中为雅可比矩阵。
步骤三：傅里叶变换，
公式(6)表示公式(5)的傅里叶变换，其中v表示频率变量，从而计算获得传递函数
步骤四：获得预测功率谱密度函数，
根据公式(7)，预测功率谱密度函数由G(v,θ)定义，分别由对应的u(t)和H(v,θ)计算得出。
S2.2采样功率谱密度函数
采样功率谱密度函数基于12阶段自回归模型计算，并定义为
最后，对于每个频率v，和G(v,θ)的所有元素分别堆叠在列向量/>和g(θ)中。
S3：混合模拟退火原理模块：
为了研究癫痫发作时脑区内部在癫痫发作过程中的状态变化，本文使用频谱DCM算法来估计cPBM模型中的参数。模型输出的似然函数为p(G)，其对数表达式为：
ln p(G) ＝F+KL(q(θ),p(θ|G)) (8)
其中KL为Kullback-Leible散度，F表示自由能，q(θ)为待估计参数θ的概率密度函数，p(θ|G)是参数θ的后验概率分布函数。对于任何给定的q(θ)，KL(q(θ),p(θ|g))≥0，因此lnp(G)≥F。所以lnp(G)可以通过最大化F来求解。F的数学表达式为：
最大化F通过EM算法求解。
对于每个模型参数向量θ，计算cPBM输出信号的预测功率谱函数，并将其与癫痫患者的EEG信号计算的采样功率谱函数进行比较。参数估计过程旨在通过EM算法迭代以最大化自由能。然而传统的频谱DCM可能出现的局部最优问题，本文在频谱DCM中提出了一种混合确定性退火DCM(H-DCM)定义如下:
其中，1/β是温度参数。传统的频谱DCM可以看作是一种特殊情况，即β＝1。当温度较高(相对较小)时(例如，β较小(相对较高))，目标函数F的形状是平滑的(相对陡峭)，通过改变目标函数的平滑程度，可以获得一个较好的初始值。因此，本发明本文引入了一种混合的模拟退火算法来解决局部的最优的问题，如图3所示。H-DCM算法通过在两个方向上调节温度(β从1.0变化到1.6(βnew＝β+0.2)，以获得相对良好的初始化；然后，温度逐渐升高，β从1.6变化到1.0(βnew＝β-0.2))以在每次最大化后寻找更好的估计。
S4：兴奋性和抑制性平衡计算：
通过上述S1-S3所述方法对印度新德里HauzKhas神经与睡眠中心的10名癫痫患者数据进行参数估计，并分别统计了癫痫发作间期、发作前期和发作期的参数变化情况，计算得到了能够对兴奋性和抑制性平衡进行量化的最佳指标。
实施例2：本发明所用的数据集采用印度新德里HauzKhas神经与睡眠中心的癫痫脑电数据集，该数据集从印度新德里HauzKhas神经与睡眠中心收集了10名癫痫患者的典型分段脑电图时间序列记录。采集过程中，根据10-20电极放置系统放置镀金头皮脑电图电极，以200Hz的采样频率进行信号采集。采集到的信号在0.5到70Hz之间过滤，然后分为发作间期、发作前和发作期。每个可下载文件夹包含50个脑电时间序列信号MAT文件，每个MAT文件由持续时间为5.12秒的1024个EEG时间序列数据样本组成。
图4显示了癫痫发作三个阶段((a)、(b)、(c)和(d)分别表示发作间期、发作前期(单峰)、发作前期(双峰)和发作期)的其中一条数据的仿真结果，其中图4(a-i)(b-i)(c-i)和(d-i)真实EEG信号(Sample EEG，表示为实线)和使用4阶Runge-Kutta方法以200Hz的采样频率求解带有估计值的方程(1)而重建的EEG信号(Reconstructed EEG)；图4(a-ii)(b-ii)(c-ii)和(d-ii)中显示了三个癫痫阶段的Sample PSD/>使用H-DCM获得的估计结果的计算得到的Predicted PSD/>Reconstructed PSD/>之间的拟合效果。为了衡量功率谱密度函数的拟合效果，本文RMSE来衡量PSD之前的拟合程度，其中Rp表示SamplePSD和Predicted PSD之间拟合程度，而Rr表示Sample PSD和Reconstructed PSD之间的拟合程度。由图可知，通过cPBM模型和H-DCM算法能够重建出该4条真实数据。
其次，为了进一步验证cPBM对真实数据的拟合效果，本发明计算了所有真实数据的RMSE和自由能，并计算了评价值以此来评价cPBM模型的仿真效果，统计结果如表2所示。由表2可知，无论是Predicted PSD和Sample PSD的拟合程度还是Reconstructed PSD和Sample PSD的拟合程度/>其都保持在一个较低值。因此，cPBM模型和H-DCM算法能够模拟癫痫发作的不同状态。
表2 H-DCM算法在癫痫不同发作阶段基于cPBM模型的拟合效果
为了探究癫痫发作的动态过度机制，表3中给出了通过本发明提出的H-DCM算法计算的不同癫痫发作阶段的各50条EEG信号的参数估计结果。为了衡量模型参数在癫痫发作过程中的变化幅度，本文引入了变化率I来衡量参数的变化大小，其定义为：
在从发作间期到发作期的动态变化过程中，本文观察到参数集{C1,C2,C5,C7,A,G,g}显著增加，参数集{C8,b}(|Iθ|>9％)。对于内部连接，结果表明，在cPBM中从发作间期到发作间期的转变可以通过锥体细胞和兴奋性中间神经元之间的连接强度增加锥体细胞和快速抑制中间神经元之间的连接强度增加(和/>)以及快速抑制中间神经的自环路连接强度的减少来解释。另一方面，有临床研究表明，对快速抑制性神经元的抑制的减小可以解释兴奋性和抑制性之间的失衡，从而导致了大脑的异常放电，并引起了患者从发作间期到发作期的转变。
表3通过H-DCM算法对三个癫痫阶段的50条数据平均估计参数(平均值±方差)/>
为了从内部神经元之间兴奋性和抑制性连接的连接强度比率来衡量E-I平衡，我们定义E-I balance指标：EPF＝C5/C8(表示从Pp到Pf的兴奋连接与Pf的自我抑制连接之间的比率)。图5显示了三个阶段的50条信号的三个指标的箱线图。从这张图中，我们可以看到发作前阶段E-I平衡的三个指标的变化，EPF可以很好地量化三个癫痫阶段，其在发作间期、发作前期和发作期的平均EPF分别为0.42、0.78以及1.23，在发作间期和发作期之间变化率增加了约190％。
以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的仅为发明的优选例，并不用来限制本发明，在不脱离本发明新型精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。
Claims (4)
1.一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法，其特征在于，所述方法包括以下步骤：
S1：建立神经元集群模型模块，
建立基于cPBM的神经元集群模型，通过龙格库塔法求解cPBM方程组获得时域信号，用以对癫痫发作的不同阶段进行模拟，
S2：建立功率谱密度函数计算模块，
用于将时域信号转化为频域信号，主要包括基于自回归模型计算真实EEG信号的采样功率谱密度函数；以及依据cPBM模型的微分方程组计算预测功率谱密度函数，
S3：建立混合模拟退火原理模块，
用于求解cPBM模型的模型参数，主要包括建立目标函数(自由能)，并使用EM算法估计模型参数使自由能最大，进而使采样功率谱密度函数和预测功率谱密度函数相拟合，
S4：兴奋性和抑制性平衡计算，
通过对模型参数的估计结果进行进一步分析，使用模型参数估计结果中的C5和C8计算得到兴奋性和抑制性平衡指标Epf。
2.根据权利要求1所述的基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法，其特征在于，所述步骤S1具体为：
cPBM模型由锥体细胞(Pp)、兴奋性神经元(Pe)、缓慢抑制性中间神经元(Ps)和快速抑制性中间神经元(Pf)构成，这四个神经元集群分别通过兴奋性连接{C1，C2，C3，C5}和抑制性连接{C4，C6，C7，C8}相互作用，和/>是两个遵循高斯概率分布的独立白噪声的外部输入，分别作用于Pp和Pf，其中/> y(t)表示模型的输出，cPBM模型采用一组连续的时间微分方程组来表示，如下公式(1)所示：
其中xi表示系统的状态变量，S(x(t))为Sigmoid函数，定义为：
其中，r＝0.56mV-1和v0＝6mV分别决定Sigmoid函数形状的陡峭程度和平移位置，e0＝6mV决定Sigmoid函数的最大值，模型的输出y(t)为：
y(t)＝x3(t) (3)
在从发作间期到发作期的过渡过程中，神经元集群中兴奋和抑制之间的平衡可能会发生变化，因此，其中上标T表示转置运算符。
3.根据权利要求1所述的基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法，其特征在于，所述步骤S2具体为：
S2.1预测功率谱密度函数，
步骤一：建立状态空间方程
公式(1)和公式(3)表示的cPBM模型可以用一组14个一阶微分方程和一个分量输出y(t)的状态空间方程，其中f(x(t)，θ)是与状态向量和θ相关联的非线性函数向量，/>代表两个输入白高斯噪声的向量，除D5，1＝Aa/C2，D13，2＝Aa和Q1，3＝1外，矩阵/>和行向量/>的值皆为零，
步骤二：线性化，
公式(4)可以通过一阶泰勒公式展开(参见公式(5))，其中为雅可比矩阵，
步骤三：傅里叶变换，
公式(6)表示公式(5)的傅里叶变换，其中v表示频率变量，从而计算获得传递函数
步骤四：获得预测功率谱密度函数，
G(v，θ)＝H(v，θ)Gu(v，θ)HH(v，θ) (7)
根据公式(7)，预测功率谱密度函数由G(v，θ)定义，分别由对应的u(t)和H(v，θ)计算得出，
S2.2采样功率谱密度函数，
采样功率谱密度函数基于12阶段自回归模型计算，并定义为
最后，对于每个频率v，和G(v，θ)的所有元素分别堆叠在列向量/>和g(θ)中。
4.根据权利要求1所述的基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法，所述步骤S3具体为：
使用频谱DCM算法来估计cPBM模型中的参数，模型输出的似然函数为p(G)，其对数表达式为：
ln p(G)＝F+KL(q(θ)，p(θ|G)) (8)
其中KL为Kullback-Leible散度，F表示自由能，q(θ)为待估计参数θ的概率密度函数，p(θ|G)是参数θ的后验概率分布函数，对于任何给定的q(θ)，KL(q(θ)，p(θ|g))≥0，因此lnp(G)≥F，所以ln p(G)通过最大化F来求解，F的数学表达式为：
最大化F通过EM算法求解，
对于每个模型参数向量θ，计算cPBM输出信号的预测功率谱函数，并将其与癫痫患者的EEG信号计算的采样功率谱函数进行比较，参数估计过程旨在通过EM算法迭代以最大化自由能，然而传统的频谱DCM可能出现的局部最优问题，在频谱DCM中提出了一种混合确定性退火DCM(H-DCM)定义如下：
其中，1/β是温度参数，传统的频谱DCM可以看作是一种特殊情况，即β＝1，当温度较高时，目标函数F的形状是平滑的，通过改变目标函数的平滑程度，获得一个较好的初始值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310447881.0A CN116439726A (zh) | 2023-04-24 | 2023-04-24 | 一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310447881.0A CN116439726A (zh) | 2023-04-24 | 2023-04-24 | 一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116439726A true CN116439726A (zh) | 2023-07-18 |
Family
ID=87125458
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310447881.0A Pending CN116439726A (zh) | 2023-04-24 | 2023-04-24 | 一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116439726A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117238485A (zh) * | 2023-11-14 | 2023-12-15 | 天津市环湖医院(天津市神经外科研究所、天津市脑系科中心医院) | 基于数据处理的智能管控系统 |
-
2023
- 2023-04-24 CN CN202310447881.0A patent/CN116439726A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117238485A (zh) * | 2023-11-14 | 2023-12-15 | 天津市环湖医院(天津市神经外科研究所、天津市脑系科中心医院) | 基于数据处理的智能管控系统 |
| CN117238485B (zh) * | 2023-11-14 | 2024-01-30 | 天津市环湖医院(天津市神经外科研究所、天津市脑系科中心医院) | 基于数据处理的智能管控系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Sangaiah et al. | An intelligent learning approach for improving ECG signal classification and arrhythmia analysis | |
| CN108714026B (zh) | 基于深度卷积神经网络和在线决策融合的细粒度心电信号分类方法 | |
| Übeyli | Adaptive neuro-fuzzy inference system for classification of ECG signals using Lyapunov exponents | |
| Übeyli | Recurrent neural networks employing Lyapunov exponents for analysis of ECG signals | |
| Übeyli | Combining recurrent neural networks with eigenvector methods for classification of ECG beats | |
| Übeyli | Statistics over features: EEG signals analysis | |
| Wang et al. | ECG signal denoising based on deep factor analysis | |
| Koçer et al. | Classifying epilepsy diseases using artificial neural networks and genetic algorithm | |
| CN111227827B (zh) | 一种基于社区划分算法的脑电图信号分析方法 | |
| CN110263924B (zh) | 一种神经元群模型的参数和状态估计方法 | |
| Emanet | ECG beat classification by using discrete wavelet transform and Random Forest algorithm | |
| CN112263253B (zh) | 基于深度学习与心电信号的抑郁症识别系统、介质及设备 | |
| CN112861625B (zh) | 一种堆叠去噪自编码器模型确定方法 | |
| CN116439726A (zh) | 一种基于混合动态因果模型的癫痫兴奋性和抑制性平衡计算方法 | |
| Mahmud et al. | Sleep apnea detection from variational mode decomposed EEG signal using a hybrid CNN-BiLSTM | |
| Chashmi et al. | An efficient and automatic ECG arrhythmia diagnosis system using DWT and HOS features and entropy-based feature selection procedure | |
| Slama et al. | Application of statistical features and multilayer neural network to automatic diagnosis of arrhythmia by ECG signals | |
| Übeyli | Statistics over features of ECG signals | |
| CN116211322A (zh) | 一种基于机器学习脑电信号的抑郁症识别方法和系统 | |
| CN111227826B (zh) | 一种基于网络模体的脑电图信号分析方法 | |
| CN115736950B (zh) | 一种基于多脑区协同相幅传递的睡眠动力学分析方法 | |
| CN116172576A (zh) | 一种基于多模块神经网络的脑电信号伪迹去除方法 | |
| Signorini | Nonlinear analysis of heart rate variability signal: physiological knowledge and diagnostic indications | |
| CN116027888A (zh) | 一种基于plv动态脑功能网络的p300意图识别方法 | |
| Goshvarpour et al. | Spectral and time based assessment of meditative heart rate signals |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |