微状态介绍及计算方法

       微状态代表了大脑在静息状态下的短暂、稳定的电压拓扑图谱。这些拓扑图谱通常持续约80-120毫秒，然后迅速转换到新的拓扑图谱。对脑电微状态的分析通常分为三步，流程如图1所示。
       一、微状态模式提取
       计算预处理后脑电信号每个时间点的全局场电位（Global Field Potential, GFP），GFP是所求时间点所有电极的标准差，抽取局部最大GFP点的脑电拓扑图作为特征进行K-Means聚类,聚类后得到K个微状态拓扑图。
       二、拟合微状态序列
       将各个时间点的脑电拓扑图与聚类后得到微状态拓扑图根据相似性进行拟合，为每个时间点分配最匹配的微状态标签，得到微状态序列。
       三、提取微状态特征
       从微状态序列中提取微状态的常见特征：
       1、微状态持续时间（Microstate Duration）：每个微状态持续的平均时间长度。
       2、微状态出现频率（Microstate Occurrence）：每个微状态在数据中出现的次数。
       3、微状态时间覆盖率（Microstate Time Coverage）：每个微状态在总分析时间中所占的比例。
       4、微状态转换概率（Microstate Transition Probability）：从一个微状态转换到另一个微状态的概率。

  下面以2秒中微状态A的计算为例进行详细说明：

原始微状态序列	AA...A(连续100个)	BB...B(连续200个)	AA...A(连续200个)	CC...C(连续1500个)
去重后微状态序列	A	B	A	C
简写	A100	B200	A200	C1500

  
  计算过程如下：

指标	计算公式	结果
(1)采样率(HZ)	1000	已知
(2)计算连续A出现的段数(次)	2	2
(3)计算连续A出现的时间点总数	A100+A200	300
(4)计算所有微状态出现时间点数	A100+A200+B200+C1500	2000
(5)计算微状态A平均持续时间(S)：	(3)/(2)/(1)	300/2/1000=0.15
(6)计算微状态A片段发生频次(次/S)：	(2)/((4)/(1))	2/(2000/1000)=1
(7)计算微状态A时间覆盖率（比率）：	(3)/(4)	300/2000=0.15
(8)计算微状态A转换到B次数	1	1
(9)计算微状态A转换到C次数	1	1
(10)计算微状态A转换到B概率:	（8）/((8)+(9))	1/2=0.5

图1：微状态分析流程（摘自[1]）

微状态分析计算举例

1、加载数据
raw = mne.io.RawArray(eeg_data, info=info)
2、抽取每个时间点GFP峰值
gfp_peaks = extract_gfp_peaks(raw)
3、聚类拟合，并可视化微状态拓扑图
ModK=ModKMeans(n_clusters=5,random_state=42)
ModK.fit(gfp_peaks, n_jobs=2)
ModK.plot()

4、提取微状态参数并可视化
segmentation = ModK.predict(raw, factor=10, half_window_size=8)
parameters = segmentation.compute_parameters()
fig, ax = plt.subplots()
x = ModK.cluster_names
y = [parameters[elt + "_timecov"] for elt in x]
ax = sns.barplot(x=x, y=y)
ax.set_xlabel("Microstates")
ax.set_ylabel("Time Coverage (ratio)")
plt.show()


注：本例使用数据为30名重度抑郁患者，SDS范围66-77分；分析方法采用python中用于计算脑电微状态的包pycrostates，详细教程见：https://pycrostates.readthedocs.io/en/latest/index.html
[1] Khanna, A., Pascual-Leone, A., Michel, C. M., & Farzan, F. (2015). Microstates in resting-state EEG: current status and future directions. Neuroscience & Biobehavioral Reviews, 49, 105-113.