CN115736950B - 一种基于多脑区协同相幅传递的睡眠动力学分析方法 - Google Patents
一种基于多脑区协同相幅传递的睡眠动力学分析方法 Download PDFInfo
- Publication number
- CN115736950B CN115736950B CN202211386244.9A CN202211386244A CN115736950B CN 115736950 B CN115736950 B CN 115736950B CN 202211386244 A CN202211386244 A CN 202211386244A CN 115736950 B CN115736950 B CN 115736950B
- Authority
- CN
- China
- Prior art keywords
- amplitude
- phase
- frequency
- signal
- cross
- 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.)
- Active
Links
- 238000012546 transfer Methods 0.000 title claims abstract description 54
- 238000004458 analytical method Methods 0.000 title claims abstract description 19
- 210000004556 brain Anatomy 0.000 claims abstract description 45
- 238000010168 coupling process Methods 0.000 claims abstract description 34
- 238000005859 coupling reaction Methods 0.000 claims abstract description 34
- 230000008878 coupling Effects 0.000 claims abstract description 32
- 230000006870 function Effects 0.000 claims abstract description 28
- 239000011159 matrix material Substances 0.000 claims abstract description 24
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 22
- 238000001228 spectrum Methods 0.000 claims abstract description 18
- 238000004364 calculation method Methods 0.000 claims abstract description 12
- 238000000528 statistical test Methods 0.000 claims abstract description 7
- 238000012360 testing method Methods 0.000 claims abstract description 6
- 238000000034 method Methods 0.000 claims description 45
- 239000013598 vector Substances 0.000 claims description 17
- 238000005070 sampling Methods 0.000 claims description 12
- 238000013507 mapping Methods 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 9
- 238000007781 pre-processing Methods 0.000 claims description 5
- 230000004424 eye movement Effects 0.000 claims description 3
- 238000012880 independent component analysis Methods 0.000 claims description 3
- 238000012847 principal component analysis method Methods 0.000 claims description 3
- 230000033764 rhythmic process Effects 0.000 claims description 3
- 238000012937 correction Methods 0.000 claims description 2
- 230000009466 transformation Effects 0.000 claims description 2
- 238000012800 visualization Methods 0.000 claims description 2
- 230000003993 interaction Effects 0.000 abstract description 9
- 238000005516 engineering process Methods 0.000 abstract description 3
- 230000001149 cognitive effect Effects 0.000 abstract description 2
- 230000001364 causal effect Effects 0.000 description 5
- 230000007246 mechanism Effects 0.000 description 5
- 230000008667 sleep stage Effects 0.000 description 5
- 230000003247 decreasing effect Effects 0.000 description 3
- 238000011002 quantification Methods 0.000 description 3
- 230000004044 response Effects 0.000 description 3
- 230000007704 transition Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 201000010099 disease Diseases 0.000 description 2
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 210000001652 frontal lobe Anatomy 0.000 description 2
- 210000002569 neuron Anatomy 0.000 description 2
- 210000000869 occipital lobe Anatomy 0.000 description 2
- 230000010355 oscillation Effects 0.000 description 2
- 230000037007 arousal Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000007177 brain activity Effects 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000005056 memory consolidation Effects 0.000 description 1
- 230000004630 mental health Effects 0.000 description 1
- 230000001537 neural effect Effects 0.000 description 1
- 230000007472 neurodevelopment Effects 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 230000035790 physiological processes and functions Effects 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000003956 synaptic plasticity Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Abstract
本发明涉及一种基于多脑区协同相幅传递的睡眠动力学分析方法，属于认知神经科学和信息技术交叉领域。本发明提出一种基于交叉相幅传递熵的多维脑电相位‑幅值耦合分析方法。首先采集多通道睡眠脑电信号并对其进行预处理，去除杂讯和干扰。对预处理后的信号进行噪声辅助多元经验模式分解，得到一系列尺度相互对准的多元本征模函数。计算每两组脑电通道间相对低频本征模函数分量与相对高频本征模函数分量间的交叉相幅传递熵，进而构造相幅耦合特征矩阵。对计算得到的特征矩阵内数值进行显著性检验，提高结果准确性。进一步将统计检验显著的相幅耦合数值平均分配到双频率耦合图谱中，直观反映大脑不同区域间的信息交互作用。
Description
技术领域
本发明涉及一种基于多脑区协同相幅传递的睡眠动力学分析方法，特别涉及一种基于交叉相幅传递熵的多维脑电信号跨频耦合分析方法，属于认知神经科学和信息技术交叉领域。
背景技术
睡眠作为一种不可或缺的生理状态，被认在对神经元发育、突触可塑性、记忆巩固、心理健康等功能中发挥作用，且在特定的睡眠阶段具有典型的大脑活动，该活动不仅在细胞层面，也存在于网络层面，研究大脑在睡眠期间不同区域神经元间的交互作用有助于深入理解大脑睡眠运行机制。目前，相幅耦合(PAC)是睡眠活动分析最关键的技术之一，该技术关注低频信号的相位与高频信号的幅值之间的交互调制，在睡眠动力学领域，常被应用于生理学或病理学研究。
睡眠脑电信号中通常包含多个频带的信号，为了进行基于神经振荡的多频率分量之间的耦合分析，需要将原始信号的不同频率成分进行分离，并选取合适的指标量化耦合。然而，现有方法侧重于单通道脑电信号的相幅调制，选用的分解方法无法同时对多通道信号频段进行有效可靠的划分，且现有的量化指标常具有统计依赖性，不能兼顾交互强度与信息流向的探讨，有关多通道间跨频交互与因果推理的研究尚未充分展开。因此，如何提出一种量化大脑睡眠状态下的多区域睡眠动力学特征方案，可以在完成多通道非线性非平稳信号的不同频率成分提取之后，实现因果特征的量化是本领域技术人员亟待解决的技术问题。
发明内容
本发明的技术解决问题是：克服现有技术的不足，提出一种基于多脑区协同相幅传递的睡眠动力学分析方法。
本发明的技术解决方案是：
一种基于多脑区协同相幅传递的睡眠动力学分析方法，该方法的步骤包括：
步骤一，多导睡眠监测仪采集受测者整晚睡眠数据，得到多通道脑电信号；
步骤二，对步骤一采集到的多通道脑电信号进行预处理；
步骤三，对步骤二预处理后的多通道脑电信号进行噪声辅助多元经验模态分解，得到本征模函数；
步骤四，对步骤三中分解得到的本征模函数进行希尔伯特变换，提取本征模函数的瞬时相位信号和瞬时幅值信号；
步骤五，根据步骤四得到的本征模函数的瞬时相位信号和瞬时幅值信号构造通道间跨层相位-幅值数据组，并计算数据组内数据间的交叉相幅传递熵；
步骤六，对步骤五计算得到的数据组内的交叉相幅传递熵进行显著性检验，保留具有统计学意义的结果；
步骤七，将步骤六得到的计算得到的统计检验显著的跨脑区交叉相幅传递熵值整理成矩阵并对矩阵进行可视化，完成基于多脑区协同相幅传递的睡眠动力学分析。
所述的步骤一中，多通道脑电信号的采样频率为fs；
所述的步骤二中，进行预处理的具体方法为：应用线性插值法去除原始脑电信号中的尖峰干扰，使用50Hz或60Hz陷波滤波器去除工频干扰，使用截止频率为0.3Hz高通滤波器消除伪影对低频的污染，采用主成分分析方法和快速独立成分分析方法估计脑电信号的独立成分，排除眼动，眨眼和心律以及电源线噪声引起的脑电图污染；
所述的步骤三中，进行噪声辅助多元经验模态分解的具体步骤如下：
(31)向N通道脑电信号x(t)＝{x1(t),x2(t),…,xN(t)}中加入M通道高斯白噪声g(t)＝{g1(t),g2(t),…,gM(t)}，构成(N+M)通道重建信号v(t)＝{x1(t),x2(t),…,xN(t),g1(t),g2(t),…,gM(t)}，t＝1,2,…,T，信号长度为T；
(32)采用Hammersley序列采样法，在(N+M-1)维球面上获得均匀采样点集，采样点集的方向角集为{θk＝[θk(1),θk(2),…,θk(N+M-1)]}，k＝1,2,…,K，K为方向角的总数量；并构建对应于方向角集{θk}的(N+M)维空间的方向向量集{zθk＝[zk(1),zk(1),…,zk(N+M)]}，k＝1,2,…,K；
(33)计算步骤(31)构成的(N+M)通道重建信号v(t)在步骤(32)构建的方向向量集{zθk}中每个方向向量zθk上的映射，得到映射集合{pθk(t)}，k＝1,2,…,K，并确定映射集合{pθk(t)}中每个映射信号的极值所对应的瞬时时刻，得到瞬时时刻集合{tΔ θk}，k＝1,2,…,K，tΔ为极值点的位置，取值范围为[1,T]；
(34)对极值点[tΔ θk,pθk(tΔ θk)]用多元样条插值函数进行插值，得到K个多元包络{eθk(t)}，k＝1,2,…,K，并计算得到的K个多元包络{eθk(t)}的均值m(t)＝[eθ1(t)+eθ2(t)+…+eθK(t)]/K；
(35)通过h(t)＝v(t)-m(t)提取固有模式函数h(t)；
(36)判断h(t)是否满足多元本征模函数(简称IMF)的判断标准，若不满足则将h(t)的内容赋值给v(t)重复步骤(32)-(35)，直至h(t)满足多元IMF判断标准，进入步骤(37)；
(37)将v(t)-h(t)的内容赋值给v(t)重复步骤(32)-(36)，直至v(t)被减至为单调或常值序列停止；
经过多次分解过程，原始(N+M)通道重建信号v(t)被分解为多个频率由高到低排列的IMFs分量{hi(t)}和余量r(t)的加和形式，i＝1,2,…,q，即
v(t)＝h1(t)+ h2(t)+…+ hq(t)+r(t) (5)
式(1)中，q表示分解出来的多IMFs层数，hi(t)＝{h1(i)(t),h2(i)(t),…,hN+M(i)(t)}，i＝1,2,…,q，和r(t)＝{r1(t),r2(t),…,rN+M(t)}，t∈[1,T]，分别对应于(N+M)通道信号的(N+M)通道IMFs分量和(N+M)通道余量，最后从(N+M)元IMFs中删除M通道噪声对应的IMFs，保留有用信号的N通道IMFs；
所述的步骤四中，对步骤三中分解得到的N通道IMFs进行希尔伯特变换，提取瞬时相位信号和瞬时幅值信号An(i)(t)＝sqrt(HT[hn(i)(t)]2+hn(i)(t)2)，其中HT(·)表示希尔伯特变换，i＝1,2,…,q，q表示分解出来的多元IMFs的层数，n＝1,2,…,N，N为脑电信号的通道总数，t∈[1,T]；
所述的步骤五中，由于经验模式分解具有倍频特性，所以分解所得本征模函数IMFs的均值频率逐层递减，构造通道间跨层相位-幅值数据组，使得数据组内相位信号所对应的IMF频率永远小于幅值信号所对应的IMF频率，即其中n＝1,2,…,N，m＝1,2,…,N，i＝1,2,…,q，j＝1,2,…,q。计算每对相位-幅值数据组中相位与幅值信号的交叉相幅传递熵(简称相cross-PATE)为：
其中，H(·)表示香农熵，t为离散时间变量，u为预测时刻，u>0，At+u为(t+u)时刻的幅值，和At d(A)分别为相位和幅值信号的相空间重构向量，At d(A)＝(At,At-τ,At-2τ,…,At-(d(A)-1)τ)，/>τ表示时间延迟；这样定义的cross-PATE以转移概率的方式，通过量化相位信号的过去信息对幅值信号未来值的贡献作为从相位到幅值方向的信息流动；
由于传递熵具有非对称性，即TEX→Y≠TEY→X，进而定义交叉幅相传递熵(简称cross-APTE)：
比较cross-PATE与cross-APTE的大小，如果cross-PATE>cross-APTE，则说明低频相位向高频幅值传递了更多的信息，可认为两者之间存在因果关系；
交叉相幅传递熵与交叉幅相传递熵的具体计算过程如下：
(51)根据Takens延时嵌入理论重建信号的状态空间，得到At d(A)＝(At,At-τ,At-2τ,…,At-(d(A)-1)τ)，其中d(A)，/>为嵌入维数，τ为嵌入延迟；
(52)依据Ragwitz准则确定这三个参数d(A)，τ的取值，具体为：Ragwitz准则通过寻找对未来状态的最佳预测值确定嵌入维数和嵌入延迟，预先给定嵌入维数d和嵌入延迟τ的取值范围，使用局部常量预测方法，最小化均方预测误差，寻找时间序列未来样本的最佳预测值，以时间序列y(t)为例，重构相空间后得到的矢量yt d(y)未来值yt+u的最大似然估计通过取其所有邻居的均值而得到，即/>U(t’)是yt d(y)邻域范围内的所有向量yt’ d(y)的集合，U(t’)＝{yt’ d(y):||yt’ d(y)-yt d(y)||≤ε}，邻域通过确定邻域半径ε或邻域内的邻居数来决定，当估计值/>与真实值yt+u的均方误差/>最小时，对应的嵌入参数(嵌入维数d和嵌入延迟τ)为最优值，其中T’为重构相空间后的信号长度；
(53)熵项的计算通过最近邻和Kraskov--Grassberger估计器，其中，最近邻技术通过统计嵌入给定空间的相邻数据点之间的距离来估计香农熵。然而，各熵项的空间维数相差较大，直接应用最近邻技术会造成错误。若固定邻居数量，式中各项的空间尺度相差较大，而每项的误差偏差取决于这些尺度，因此误差不会相互抵消，而是累积。为了解决这个问题，Kraskov-/>-Grassberger估计器在最近邻基础上，固定最高维空间中的邻居数量λ(λ取值可由用户自行设置，默认λ＝4)，可获得第λ近邻的距离L，将该距离投影到低维空间得L’，作为低维空间项的邻居范围，进而寻找各低维项的邻居数量η。因此交叉相幅传递熵通过下式计算：
其中ψ为digamma函数，最高维空间由At+u，At d(A)和所张成，η(At d(A))、η(At+u,At d(A))、/>分别代表各低维空间中在距离L’内的邻居数，<·>t表示在不同时间点上取均值；熵值越大，说明该方向上的耦合程度越强，传递的信息量越大；
所述的步骤六中，进行显著性检验的方法为：使用混洗方法生成100个独立的替代数据，比较实际脑电信号的传递熵TE值与生成替代数据的传递熵TEshuff值，以mean(TEshuff)和σshuff分别表示替代数据传递熵值集合的均值和标准差，计算z分数：
z-score＝[TE-mean(TEshuff)]/σshuff (9)
保留具有统计学意义(Bonferroni校正后p<0.05)的传递熵TE数值以及所对应的相位-幅值数据组；
所述的步骤七中，进行可视化的方法为：将步骤六得到的统计检验显著的跨脑区cross-PATE整理成矩阵，矩阵中行与相位-幅值数据组中的幅值相对应，列与相位相对应，考虑相位与幅值之间具有频率差，因此TE矩阵是一个上三角矩阵，将矩阵中的TE值与瞬时频率做匹配，均匀分散到双频耦合图谱上，双频耦合图谱的横轴对应于相位频率，即低频，频率范围设为[f1,f2]，频率分辨率设为α；纵轴对应于幅值频率，即高频，频率范围设为[F1,F2]，频率分辨率设为β，将带有瞬时频率标签的TE矩阵中数值均匀分散到频率分辨率为(f2-f1)/α×(F2-F1)/β的双频耦合图谱中，即得到从A脑区到B脑区的cross-PATE图谱，同样的，计算得到从A脑区到B脑区的cross-APTE图谱，比较同一相位-幅值数据组中cross-PATE和cross-APTE大小，以区分从A脑区方向到B脑区方向的驱动和响应成分，另外，还可以计算从B脑区到A脑区方向的cross-PATE和cross-APTE，进而判断A，B两脑区之间的信息流动；
一般情况下，同一相位-幅值数据组中cross-PATE>cross-APTE，这与所讨论的实际生理意义相一致，所以通常我们只需要比较两脑区间两个方向上的cross-PATE大小，以判断跨脑区的信息流动方向和强度。
有益效果
(1)本发明的步骤三中使用噪声辅助经验模式分解(NA-MEMD)分解多通道脑电信号，该方法是一种面向多通道脑电信号多尺度分解的方法，能够自适应将脑电信号分解为频率相互对准的多维本征模态分量，使得多变量脑电信号经由分解之后，所得到的各分量在相同信号源下分解为多个固有振荡，在空间上，各本征模态分量的频率分布相互对准，同时保留了脑电信号本身的非线性、非平稳以及动态变化属性。
(2)本发明的步骤五中使用传递熵(TE)作为相幅耦合的量化方式，基于信息论的传递熵具有模型无关性和非参数特性，其通过量化相关信源对于被观测信源不确定度减少的贡献，作为该信源与被观测信源间的信息交互，在量化交互作用强度的同时，能够提供信息交互作用的方向，可用作因果机制量化指标。
(3)本发明的步骤六中对计算得到的特征矩阵内数值进行显著性检验，通过将原数据生成的TE与替代数据生成的TEshuff之间的差异转化为z分数的方式进行显著性分析能够有效去除虚假的传递熵TE，提高结果准确性。且步骤七中进一步将统计检验显著的相幅耦合数值平均分配到双频率耦合图谱中，通过可视化的形式展现脑电数据的相幅耦合分析结果，直观反映大脑不同区域间的信息交互作用。
(4)本发明将NA-MEMD与传递熵方法的优势相结合，提出交叉相幅传递熵(cross-PATE)方法，用于量化大脑睡眠状态下不同区域神经元间的交互作用，探讨多脑区协同睡眠动力学机制，适用于提取脑电信号耦合特征，高效可靠，易于软件化。
附图说明
图1为本发明的方法流程示意图；
图2为睡眠脑电走势图及NA-MEMD分解后IMFs走势图，(a)6通道睡眠脑电的原始走势图及其幅频谱；(b)NA-MEMD分解后得到5组本征模函数(IMFs)及其幅频谱；
图3为10名健康被试在五个睡眠阶段(Awake、N1、N2、N3和REM)中额叶(F)和枕叶(O)通道间的平均定向相幅双频耦合图。
具体实施方式
为了量化大脑不同区域间的动态耦合以及信息流动方向，本发明提出基于噪声辅助多元经验模式分解和传递熵理论的交叉相幅传递熵脑电耦合分析方法，图一为本发明的技术流程图，详细过程如下：
步骤一，受试者按照实验要求完成实验操作，应用到多导睡眠监测仪采集受测者整晚睡眠数据，脑电(EEG)信号采样频率为fs，存储在计算机中。
步骤二，对采集到的多通道脑电信号进行如下预处理：应用线性插值法去除原始信号中的尖峰干扰，使用50Hz陷波滤波器去除工频干扰，使用截止频率为0.3Hz高通滤波器消除伪影对低频的污染，采用主成分分析方法和快速独立成分分析方法估计脑电信号的独立成分，排除眼动，眨眼和心律以及电源线噪声引起的脑电图污染。
步骤三，对预处理后的多维脑电信号进行噪声辅助多元经验模态分解(NA-MEMD)，将信号分解为多个多元本征模函数和余量分量，具体步骤如下：
1)向N通道脑电信号x(t)＝{x1(t),x2(t),…,xN(t)}中加入M通道高斯白噪声g(t)＝{g1(t),g2(t),…,gM(t)}，构成(N+M)通道重建信号v(t)＝{x1(t),x2(t),…,xN(t),g1(t),g2(t),…,gM(t)}，t＝1,2,…,T，信号长度为T。
2)采用Hammersley序列采样法，在(N+M-1)维球面上获得均匀采样点集，方向角为{θk＝[θk(1),θk(2),…,θk(N+M-1)]}，k＝1,2,…,K，K为方向角的总数量；构建对应于角集{θk}的(N+M)维空间的方向向量集{zθk＝[zk(1),zk(1),…,zk(N+M)]}，k＝1,2,…,K。
3)计算输入信号v(t)在每个方向向量zθk上的映射{pθk(t)}，k＝1,2,…,K。确定各方向的投影信号{pθk(t)}的极值所对应的瞬时时刻{tΔ θk}，k＝1,2,…,K，tΔ为极值点的位置，取值范围为[1,T]，T为信号长度。
4)用多元样条插值函数插值极值点[tΔ θk,pθk(tΔ θk)]，得到K个多元包络{eθk(t)}，k＝1,2,…,K。对球空间K个方向，计算多元包络线的均值m(t)＝[eθ1(t)+eθ2(t)+…+eθK(t)]/K。
5)通过h(t)＝v(t)-m(t)提取固有模式函数h(t)，判断h(t)是否满足多元本征模函数(IMF)的判断标准，若不满足则将h(t)当做第(2)步的输入信号，继续第(2)～(5)步迭代；若h(t)满足多元IMF判断标准，则将v(t)-h(t)的结果当做第(2)步的输入信号，继续第(2)～(5)步迭代，提取新的多元IMF分量。上述过程直至原信号被减至为单调或常值序列停止。
经过多次分解过程，原始(N+M)通道重建信号v(t)被分解为多个频率由高到低排列的IMFs分量{hi(t)}和余量r(t)的加和形式，i＝1,2,…,q，即
v(t)＝h1(t)+ h2(t)+…+ hq(t)+r(t) (10)
式(1)中，q表示分解出来的多元IMFs层数，hi(t)＝{h1(i)(t),h2(i)(t),…,hN+M(i)(t)},i＝1,2,…,q，和r(t)＝{r1(t),r2(t),…,rN+M(t)}，t∈[1,T]，分别对应于(N+M)通道信号的(N+M)通道IMFs分量和(N+M)通道余量，最后从(N+M)元IMFs中删除M通道噪声对应的IMFs，保留有用信号的N通道IMFs；
步骤四，对步骤三中分解得到的N通道IMFs进行希尔伯特变换，提取瞬时相位信号和瞬时幅值信号An(i)(t)＝sqrt(HT[hn(i)(t)]2+hn(i)(t)2)，其中HT(·)表示希尔伯特变换，i＝1,2,…,q，q表示分解出来的多元IMFs层数，n＝1,2,…,N，N为脑电信号的通道总数，t∈[1,T]。
步骤五，由于经验模式分解具有倍频特性，所以分解所得本征模函数IMFs的均值频率逐层递减。构造通道间跨层相位-幅值数据组，使得组内相位信号所对应的IMF频率永远小于幅值信号所对应的IMF频率，即 其中n＝1,2,…,N，m＝1,2,…,N，i＝1,2,…,q，j＝1,2,…,q。计算每对相位-幅值数据组中相位与幅值信号的交叉相幅传递熵：
其中，H(·)表示香农熵，t为离散时间变量，u为预测时刻，u>0，因此At+u为t+u时刻的幅值预测值。和At d(A)分别为相位和幅值信号的相空间重构向量At d(A)＝(At,At-τ,At-2τ,…,At-(d(A)-1)τ)，/>这样定义的cross-PATE以转移概率的方式，通过量化相位信号的过去信息对幅值信号未来值预测的贡献作为从相位到幅值方向的信息流动。
由于传递熵具有非对称性，即TEX→Y≠TEY→X，进而定义幅值相位传递熵：
比较cross-PATE与cross-APTE的大小，如果cross-PATE>cross-APTE，则说明低频相位向高频幅值传递了更多的信息，可认为两者之间存在因果关系。
(2)式和(3)式的具体计算过程如下：
1)根据Takens延时嵌入理论重建信号的状态空间，得到At d(A)＝(At,At-τ,At-2τ,…,At-(d(A)-1)τ)，其中d(A)，/>为嵌入维数，τ为嵌入延迟，依据如下Ragwitz准则确定这三个参数的取值。
2)Ragwitz准则通过寻找对未来状态的最佳预测值而确定嵌入维数和嵌入延迟。预先给定嵌入维数d和嵌入延迟τ的取值范围，使用局部常量预测方法，最小化均方预测误差，寻找时间序列未来样本的最佳预测值。以时间序列y(t)为例，重构相空间后得到的矢量yt d(y)未来值yt+u的最大似然估计通过取其所有邻居的均值而得到，即U(t’)是yt d(y)邻域范围内的所有向量yt’ d(y)的集合，U(t’)＝{yt’ d(y):||yt’ d(y)-yt d(y)||≤ε}，邻域通过确定邻域半径ε或邻域内的邻居数来决定，当估计值/>与真实值yt+u的均方误差/> 最小时，对应的嵌入参数(嵌入维数d和嵌入延迟τ)为最优值，其中T’为重构相空间后的信号长度；
3)式(2)和式(3)中联合熵和边缘熵的计算可以通过最近邻和Kraskov--Grassberger估计器。其中，最近邻技术通过统计嵌入给定空间的相邻数据点之间的距离来估计香农熵。然而，交叉相幅传递熵和交叉幅相传递熵计算公式中各项的空间维数相差较大，直接应用最近邻技术会造成错误。若固定邻居数量，式中各项的空间尺度相差较大，而每项的误差偏差取决于这些尺度，因此误差不会相互抵消，而是累积。为了解决这个问题，Kraskov-/>-Grassberger估计器在最近邻基础上，只固定最高维空间中的邻居数量λ(λ取值可由用户自行设置，默认λ＝4)，可获得第λ近邻的距离L，将该距离投影到低维空间，作为低维空间项的邻居范围L’，进而寻找各低维项的邻居数量η。因此交叉相幅传递熵计算公式的传递熵可通过下式计算：
其中ψ为digamma函数，最高维空间由At+u，At d(A)和所张成，η(At d(A))、η(At+u,At d(A))、/>分别代表各低维空间中在距离L’内的邻居数，<·>t表示在不同时间点上取均值。熵值越大，说明该方向上的耦合程度越强，传递的信息量越大。
步骤六，考虑到实际计算时数据量是有限的，因此进行熵值估计时会出现偏差，为了保证上一步计算得到的TE值的真实性和有效性，同时剔除无效值，使用混洗方法生成100个独立的替代数据：在相位变化为2π的一个周期内随机打乱时间序列。比较实际脑电信号的TE值与生成替代数据的TEshuff值，以mean(TEshuff)和σshuff分别表示替代数据TE值集合的均值和标准差，计算z分数：
z-score＝[TE-mean(TEshuff)]/σshuff (14)
保留具有统计学意义(Bonferroni校正后p<0.05)的TE数值以及所对应的相位-幅值数据组。
步骤七，将上一步计算得到的统计检验显著的跨脑区cross-PATE整理成矩阵，矩阵中行与相位-幅值数据组中的幅值相对应，列与相位相对应，考虑相位与幅值之间具有频率差，因此TE矩阵是一个上三角矩阵。将矩阵中的TE值与瞬时频率做匹配，均匀分散到双频耦合图谱上(如图1最下面图形所示)。双频耦合图谱的横轴对应于相位频率，即低频，频率范围设为[f1,f2]，频率分辨率设为α；纵轴对应于幅值频率，即高频，频率范围设为[F1,F2]，频率分辨率设为β，将带有瞬时频率标签的TE矩阵中数值均匀分散到频率分辨率为(f2-f1)/α×(F2-F1)/β的双频耦合图谱中，即得到从A脑区到B脑区的cross-PATE图谱。同样的也可以计算得到从A脑区到B脑区的cross-APTE图谱。比较同一相位-幅值数据组中cross-PATE和cross-APTE大小，以区分从A脑区方向到B脑区方向的驱动和响应成分。另外，还可以计算从B脑区到A脑区方向的cross-PATE和cross-APTE，进而判断A，B两脑区之间的信息流动。
一般情况下，同一相位-幅值数据组中cross-PATE>cross-APTE，这与所讨论的实际生理意义相一致，所以通常我们只需要比较两脑区间两个方向上的cross-PATE大小，以判断跨脑区的信息流动方向和强度。当然，同一相位-幅值数据组中cross-PATE>cross-APTE的关系需要被验证，只有在该不等式真实成立的条件下，才能舍弃该部分的讨论，否则相位与幅值的驱动与响应关系也需要被讨论。
对所有多通道脑电信号分解得到的IMFs重复以上步骤，得到多脑区协同的交叉频率相幅耦合结果。
以睡眠数据库ISRUC-Sleep的研究结果为例，该数据库中脑电信号采样率为200Hz，包含六个通道(F3-A2、C3-A2、O1-A2、F4-A1、C4-A1、O2-A1)，脑电图走势及频谱如图例2(a)所示。经由NA-MEMD分解后，得到5组IMFs(图2(b))，每行的IMFs频率相互对准，不同通道IMFs之间的幅频谱高度重合，从上到下IMFs的频率逐渐降低。图3为ISRUC-Sleep数据库中10名健康被试在五个睡眠阶段(Awake、N1、N2、N3和REM)中额叶(F)和枕叶(O)通道间的平均定向相幅双频耦合图。可以发现，健康成年人在睡眠觉醒期和睡眠N1期存在从脑枕叶到脑额叶方向的相幅耦合传递，而在睡眠的其他阶段，相幅耦合传递方向转变为从脑额叶到脑枕叶方向。同时，随着睡眠进入深睡期，δ-θ/α相幅耦合逐渐增强。通过分析不同睡眠阶段的跨脑区信息流动，得出的显著特征可用于睡眠分期的模型构造，且根据以上分析结果，我们可以对健康人睡眠期间的大脑动力学情况进行归纳，进而反演睡眠期间大脑运行机制，作为证据支持睡眠的潜在功能。同时可将本发明的方法用于患病人群，以健康人数据作为对照，对疾病的发病原因、治疗方法进行探究。基于交叉相幅传递熵的耦合分析方法可以作为提取大脑耦合特征进而反演大脑运行机制的可靠分析工具。
Claims (1)
1.一种基于多脑区协同相幅传递的睡眠动力学分析方法，其特征在于该方法的步骤包括：
步骤一，多导睡眠监测仪采集受测者整晚睡眠数据，得到多通道脑电信号；
步骤二，对步骤一采集到的多通道脑电信号进行预处理；
步骤三，对步骤二预处理后的多通道脑电信号进行噪声辅助多元经验模态分解，得到本征模函数；
步骤四，对步骤三中分解得到的本征模函数进行希尔伯特变换，提取本征模函数的瞬时相位信号和瞬时幅值信号；
步骤五，根据步骤四得到的本征模函数的瞬时相位信号和瞬时幅值信号构造通道间跨层相位-幅值数据组，并计算数据组内数据间的交叉相幅传递熵；
步骤六，对步骤五计算得到的数据组内的交叉相幅传递熵进行显著性检验，保留具有统计学意义的结果；
步骤七，将步骤六得到的计算得到的统计检验显著的跨脑区交叉相幅传递熵值整理成矩阵并对矩阵进行可视化，完成基于多脑区协同相幅传递的睡眠动力学分析；
所述的步骤一中，多通道脑电信号的采样频率为fs；
所述的步骤二中，进行预处理的具体方法为：应用线性插值法去除原始脑电信号中的尖峰干扰，使用50Hz或60Hz陷波滤波器去除工频干扰，使用截止频率为0.3Hz高通滤波器消除伪影对低频的污染，采用主成分分析方法和快速独立成分分析方法估计脑电信号的独立成分，排除眼动，眨眼和心律以及电源线噪声引起的脑电图污染；
所述的步骤三中，进行噪声辅助多元经验模态分解的具体步骤如下：
(31)向N通道脑电信号x(t)＝{x1(t),x2(t),…,xN(t)}中加入M通道高斯白噪声g(t)＝{g1(t),g2(t),…,gM(t)}，构成(N+M)通道重建信号v(t)＝{x1(t),x2(t),…,xN(t),g1(t),g2(t),…,gM(t)}，t＝1,2,…,T，信号长度为T；
(32)采用Hammersley序列采样法，在(N+M-1)维球面上获得均匀采样点集，采样点集的方向角集为{θk＝[θk(1),θk(2),…,θk(N+M-1)]}，k＝1,2,…,K，K为方向角的总数量；并构建对应于方向角集{θk}的(N+M)维空间的方向向量集{zθk＝[zk(1),zk(1),…,zk(N+M)]}，k＝1,2,…,K；
(33)计算步骤(31)构成的(N+M)通道重建信号v(t)在步骤(32)构建的方向向量集{zθk}中每个方向向量zθk上的映射，得到映射集合{pθk(t)}，k＝1,2,…,K，并确定映射集合{pθk(t)}中每个映射信号的极值所对应的瞬时时刻，得到瞬时时刻集合{tΔ θk}，k＝1,2,…,K，tΔ为极值点的位置，取值范围为[1,T]；
(34)对极值点[tΔ θk,pθk(tΔ θk)]用多元样条插值函数进行插值，得到K个多元包络{eθk(t)}，k＝1,2,…,K，并计算得到的K个多元包络{eθk(t)}的均值m(t)＝[eθ1(t)+eθ2(t)+…+eθK(t)]/K；
(35)通过h(t)＝v(t)-m(t)提取固有模式函数h(t)；
(36)判断h(t)是否满足多元本征模函数(简称IMF)的判断标准，若不满足则将h(t)的内容赋值给v(t)重复步骤(32)-(35)，直至h(t)满足多元IMF判断标准，进入步骤(37)；
(37)将v(t)-h(t)的内容赋值给v(t)重复步骤(32)-(36)，直至重建信号v(t)被减至为单调或常值序列停止；
经过多次分解过程，原始(N+M)通道重建信号v(t)被分解为多个频率由高到低排列的IMFs{hi(t)}和余量r(t)的加和形式，i＝1,2,…,q，即
v(t)＝h1(t)+h2(t)+…+hq(t)+r(t) (1)
式(1)中，q表示分解出来的IMF层数，hi(t)＝{h1(i)(t),h2(i)(t),…,hN+M(i)(t)}，i＝1,2,…,q，和r(t)＝{r1(t),r2(t),…,rN+M(t)}，t∈[1,T]，分别对应于(N+M)通道信号的(N+M)元本IMFs分量和(N+M)通道余量，最后从(N+M)元IMFs中删除M通道噪声对应的IMFs，保留有用信号的N通道IMFs；
所述的步骤四中，对步骤三中分解得到的N通道IMFs进行希尔伯特变换，提取瞬时相位信号和瞬时幅值信号An(i)(t)＝sqrt(HT[hn(i)(t)]2+hn(i)(t)2)，其中HT(·)表示希尔伯特变换，i＝1,2,…,q，q表示分解出来的多元IMFs的层数，n＝1,2,…,N，N为脑电信号的通道总数，t∈[1,T]；
所述的步骤五中，数据组内相位信号所对应的IMF频率小于幅值信号所对应的IMF频率，即其中n＝1,2,…,N，m＝1,2,…,N，i＝1,2,…,q，j＝1,2,…,q；
所述的步骤五中，每对相位-幅值数据组中相位与幅值信号的交叉相幅传递熵，为：
其中，H(·)表示香农熵，t为离散时间变量，u为预测时刻，u>0，At+u为(t+u)时刻的幅值，和At d(A)分别为相位和幅值信号的相空间重构向量，At d(A)＝(At,At-τ,At-2τ,…,At-(d(A)-1)τ)，/>τ表示时间延迟；
交叉幅相传递熵为：
交叉相幅传递熵与交叉幅相传递熵的具体计算过程如下：
(51)根据Takens延时嵌入理论重建信号的状态空间，得到At d(A)＝(At,At-τ,At-2τ,…,At-(d(A)-1)τ)，其中d(A)，/>为嵌入维数，τ为嵌入延迟；
(52)依据Ragwitz准则确定这三个参数d(A)，τ的取值，具体为：Ragwitz准则通过寻找对未来状态的最佳预测值确定嵌入维数和嵌入延迟，预先给定嵌入维数d和嵌入延迟τ的取值范围，使用局部常量预测方法，最小化均方预测误差，寻找时间序列未来样本的最佳预测值，以时间序列y(t)为例，重构相空间后得到的矢量yt d(y)未来值yt+u的最大似然估计通过取其所有邻居的均值而得到，即/>U(t’)是yt d(y)邻域范围内的所有向量yt’ d(y)的集合，U(t’)＝{yt’ d(y):||yt’ d(y)-yt d(y)||≤ε}，邻域通过确定邻域半径ε或邻域内的邻居数来决定，当估计值/>与真实值yt+u的均方误差/>最小时，对应的嵌入参数(嵌入维数d和嵌入延迟τ)为最优值，其中T’为重构相空间后的信号长度；
(53)熵项的计算通过最近邻和Kraskov--Grassberger估计器，Kraskov--Grassberger估计器在最近邻基础上，固定最高维空间中的邻居数量λ，获得第λ近邻的距离L，将该距离投影到低维空间得L’，作为低维空间项的邻居范围，进而寻找各低维项的邻居数量η；
所述的步骤六中，进行显著性检验的方法为：使用混洗方法生成100个独立的替代数据，比较实际脑电信号的传递熵TE值与生成替代数据的传递熵TEshuff值，以mean(TEshuff)和σshuff分别表示替代数据传递熵值集合的均值和标准差，计算z分数：
z-score＝[TE-mean(TEshuff)]/σshuff (4)
保留具有统计学意义，即Bonferroni校正后p<0.05的传递熵值以及所对应的相位-幅值数据组；
所述的步骤七中，进行可视化的方法为：将步骤六得到的统计检验显著的跨脑区cross-PATE整理成矩阵，矩阵中行与相位-幅值数据组中的幅值相对应，列与相位相对应，将矩阵中的传递熵值与瞬时频率做匹配，均匀分散到双频耦合图谱上，双频耦合图谱的横轴对应于相位频率，即低频，频率范围设为[f1,f2]，频率分辨率设为α；纵轴对应于幅值频率，即高频，频率范围设为[F1,F2]，频率分辨率设为β，将带有瞬时频率标签的传递熵矩阵中数值均匀分散到频率分辨率为(f2-f1)/α×(F2-F1)/β的双频耦合图谱中，即得到从A脑区到B脑区的cross-PATE图谱。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211386244.9A CN115736950B (zh) | 2022-11-07 | 2022-11-07 | 一种基于多脑区协同相幅传递的睡眠动力学分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211386244.9A CN115736950B (zh) | 2022-11-07 | 2022-11-07 | 一种基于多脑区协同相幅传递的睡眠动力学分析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115736950A CN115736950A (zh) | 2023-03-07 |
| CN115736950B true CN115736950B (zh) | 2024-02-09 |
Family
ID=85357937
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211386244.9A Active CN115736950B (zh) | 2022-11-07 | 2022-11-07 | 一种基于多脑区协同相幅传递的睡眠动力学分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115736950B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115982574B (zh) * | 2023-03-20 | 2023-06-09 | 北京理工大学 | 基于频率限定掩模因果分解的脑电信息流向特征提取方法 |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104571533A (zh) * | 2015-02-10 | 2015-04-29 | 北京理工大学 | 一种基于脑机接口技术的装置和方法 |
| CN106491083A (zh) * | 2016-10-11 | 2017-03-15 | 天津大学 | 用于脑状态监测的头戴式智能穿戴电极数量优化法及应用 |
| WO2017205734A1 (en) * | 2016-05-26 | 2017-11-30 | University Of Washington | Reducing sensor noise in multichannel arrays using oversampled temporal projection and associated systems and methods |
| CN110367980A (zh) * | 2019-07-10 | 2019-10-25 | 南京邮电大学 | 基于多元经验模态分解的脑电信号情绪识别方法 |
| WO2021034351A1 (en) * | 2019-08-22 | 2021-02-25 | Hecox Kurt E | Systems and methods for seizure detection based on changes in electroencephalogram (eeg) non-linearities |
| CN112617862A (zh) * | 2020-01-15 | 2021-04-09 | 燕山大学 | 一种基于多通道神经信号分析的信号间耦合强度判断方法 |
| CN112827039A (zh) * | 2021-03-01 | 2021-05-25 | 中山大学 | 麻醉呼吸机的控制方法、装置及可调控麻醉呼吸机系统 |
| CN113158793A (zh) * | 2021-03-15 | 2021-07-23 | 东北电力大学 | 一种基于多特征融合的多类运动想象脑电信号识别方法 |
| CN113558637A (zh) * | 2021-07-05 | 2021-10-29 | 杭州电子科技大学 | 一种基于相位传递熵的音乐感知下脑网络构建方法 |
| CN113712571A (zh) * | 2021-06-18 | 2021-11-30 | 陕西师范大学 | 一种基于Rényi相位传递熵和轻量级卷积神经网络的异常脑电信号检测方法 |
| CN114343635A (zh) * | 2021-12-06 | 2022-04-15 | 北京理工大学 | 一种基于变分相幅耦合的情绪识别方法及装置 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9545219B2 (en) * | 2012-11-16 | 2017-01-17 | University Of Manitoba | Acoustic system and methodology for identifying the risk of obstructive sleep apnea during wakefulness |
| JP2016520375A (ja) * | 2013-04-23 | 2016-07-14 | ザ ジェネラル ホスピタル コーポレイション | 脳のコヒーレンスおよびシンクロニーの指標を使用して麻酔および鎮静を監視するためのシステムおよび方法 |
| US11406316B2 (en) * | 2018-02-14 | 2022-08-09 | Cerenion Oy | Apparatus and method for electroencephalographic measurement |
| CN110101384B (zh) * | 2019-04-22 | 2022-01-28 | 自然资源部第一海洋研究所 | 用于复杂网络的功能性网络分析系统及分析方法 |
-
2022
- 2022-11-07 CN CN202211386244.9A patent/CN115736950B/zh active Active
Patent Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104571533A (zh) * | 2015-02-10 | 2015-04-29 | 北京理工大学 | 一种基于脑机接口技术的装置和方法 |
| WO2017205734A1 (en) * | 2016-05-26 | 2017-11-30 | University Of Washington | Reducing sensor noise in multichannel arrays using oversampled temporal projection and associated systems and methods |
| CN106491083A (zh) * | 2016-10-11 | 2017-03-15 | 天津大学 | 用于脑状态监测的头戴式智能穿戴电极数量优化法及应用 |
| CN110367980A (zh) * | 2019-07-10 | 2019-10-25 | 南京邮电大学 | 基于多元经验模态分解的脑电信号情绪识别方法 |
| WO2021034351A1 (en) * | 2019-08-22 | 2021-02-25 | Hecox Kurt E | Systems and methods for seizure detection based on changes in electroencephalogram (eeg) non-linearities |
| CN112617862A (zh) * | 2020-01-15 | 2021-04-09 | 燕山大学 | 一种基于多通道神经信号分析的信号间耦合强度判断方法 |
| CN112827039A (zh) * | 2021-03-01 | 2021-05-25 | 中山大学 | 麻醉呼吸机的控制方法、装置及可调控麻醉呼吸机系统 |
| CN113158793A (zh) * | 2021-03-15 | 2021-07-23 | 东北电力大学 | 一种基于多特征融合的多类运动想象脑电信号识别方法 |
| CN113712571A (zh) * | 2021-06-18 | 2021-11-30 | 陕西师范大学 | 一种基于Rényi相位传递熵和轻量级卷积神经网络的异常脑电信号检测方法 |
| CN113558637A (zh) * | 2021-07-05 | 2021-10-29 | 杭州电子科技大学 | 一种基于相位传递熵的音乐感知下脑网络构建方法 |
| CN114343635A (zh) * | 2021-12-06 | 2022-04-15 | 北京理工大学 | 一种基于变分相幅耦合的情绪识别方法及装置 |
Non-Patent Citations (3)
| Title |
|---|
| Wang, Yufei ; Shi, Wenbin and Yeh, Chien-Hung.Sleep Dynamic Analysis Technology Based on Cross-Phase-Amplitude Transfer Entropy in Multiple Brain Regions.《Annu Int Conf IEEE Eng Med Biol Soc.》.2022,正文第I-IV部分. * |
| 基于排序模式的脑电信号非线性特征分析方法及应用;李帅飞;《中国优秀硕士学位论文全文数据库基础科学辑》;全文 * |
| 自主动作和电刺激下的脑电/肌电相干性对比分析研究;傅炜东;《中国优秀硕士学位论文全文数据库医药卫生科技辑》;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115736950A (zh) | 2023-03-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wang et al. | Automatic epileptic seizure detection in EEG signals using multi-domain feature extraction and nonlinear analysis | |
| Wang et al. | Single slice based detection for Alzheimer’s disease via wavelet entropy and multilayer perceptron trained by biogeography-based optimization | |
| Henriques et al. | Nonlinear methods most applied to heart-rate time series: a review | |
| Betzel et al. | Structural, geometric and genetic factors predict interregional brain connectivity patterns probed by electrocorticography | |
| Mursalin et al. | Automated epileptic seizure detection using improved correlation-based feature selection with random forest classifier | |
| Subasi et al. | Classification of EEG signals using neural network and logistic regression | |
| Sunil Kumar et al. | Bio-signals Compression Using Auto Encoder | |
| Ruiz-Gómez et al. | Computational modeling of the effects of EEG volume conduction on functional connectivity metrics. Application to Alzheimer’s disease continuum | |
| Huggins et al. | Deep learning of resting-state electroencephalogram signals for three-class classification of Alzheimer’s disease, mild cognitive impairment and healthy ageing | |
| Mahjoub et al. | Epileptic seizure detection on EEG signals using machine learning techniques and advanced preprocessing methods | |
| You et al. | Automatic focal and non-focal EEG detection using entropy-based features from flexible analytic wavelet transform | |
| Ferrari et al. | Dealing with confounders and outliers in classification medical studies: The Autism Spectrum Disorders case study | |
| CN115736950B (zh) | 一种基于多脑区协同相幅传递的睡眠动力学分析方法 | |
| Khare et al. | SchizoNET: a robust and accurate Margenau–Hill time-frequency distribution based deep neural network model for schizophrenia detection using EEG signals | |
| Arizmendi et al. | Automated classification of brain tumours from short echo time in vivo MRS data using Gaussian Decomposition and Bayesian Neural Networks | |
| Yang et al. | Link prediction in brain networks based on a hierarchical random graph model | |
| CN113017627A (zh) | 一种基于双通道相位同步特征融合的抑郁症和双相障碍脑网络分析方法 | |
| Lin et al. | Utilizing transfer learning of pre-trained AlexNet and relevance vector machine for regression for predicting healthy older adult’s brain age from structural MRI | |
| CN114343635A (zh) | 一种基于变分相幅耦合的情绪识别方法及装置 | |
| Cordes et al. | Energy-period profiles of brain networks in group fMRI resting-state data: a comparison of empirical mode decomposition with the short-time fourier transform and the discrete wavelet transform | |
| Cleatus et al. | Epileptic seizure detection using spectral transformation and convolutional neural networks | |
| Liu et al. | Automatic sleep staging algorithm based on random forest and hidden markov model | |
| Santos Febles et al. | Machine learning techniques for the diagnosis of schizophrenia based on event-related potentials | |
| Gu et al. | Detection of Attention Deficit Hyperactivity Disorder in children using CEEMDAN-based cross frequency symbolic convergent cross mapping | |
| Zhang et al. | Four-classes human emotion recognition via entropy characteristic and random Forest |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |