python实现跨频耦合：PLV、MV L、MI

目的：介绍跨频耦合的背景及常用的计算方法

1、背景

1. 神经科学的核心问题之一是神经活动如何在不同的空间和时间尺度上协调。
2. 跨频率耦合(Cross-frequency coupling, CFC) 已被提出在时间和空间尺度上协调神经动力学。
   1. 时间尺度：神经振荡之间的相幅耦合。例如，低频振荡的相位可以影响高频振荡的振幅，从而调节神经活动的时序，实现信息处理和集成。
   2. 空间尺度：神经振荡之间的空间耦合。例如，不同频率范围内的神经振荡可以在不同脑区之间相互调节，形成神经网络，实现信息的传递和整合。
3. EEG信号的CFC是指不同频率范围内的神经振荡之间存在相互调节和交互作用的现象。换句话说，CFC实际上是不同频带EEG的幅度/频率/相位之间的关系。常见的耦合关系如下：【慢震荡（Slow Oscillation，SO）；快震荡（ Fast Oscillation，FO ）】
   1. 相位-幅度耦合(phase-amplitude coupling，PAC)：SO的相位与FO幅度之间的耦合，也称为“嵌套”。目前研究比较多的是θ - γ和α - γ之间的PAC。
   2. 相位-频率耦合(phase- frequency coupling，PFC)：SO的相位和FO的频率之间的耦合。
   3. 相位-相位耦合(phase-phase coupling，PPC)：SO和FO相位之间的耦合。

2、CFC在神经疾病中的应用

跨频耦合作为一种神经机制，在多种神经疾病中都得到了广泛研究：

1. 癫痫：研究表明[1]， 跨频频率耦合是定位癫痫组织的有用生物标志物。
2. 帕金森病：研究表明[2]，theta-gamma跨频耦合可以作为帕金森疾病的生物标志，并有助于对帕金森疾病的早期诊断和治疗进行监测。
3. 精神分裂症：研究表明[3]， theta-gamma跨频耦合减弱与认知障碍症状的严重程度有关。
4. 抑郁症：研究表明[4]，抑郁症患者与健康对照组相比，前额叶皮层的theta-gamma跨频耦合存在显著差异。

3、CFC常用计算方法 [5]

1. 锁相值（phase-locking value，PLV）
2. 平均向量长度（mean vector length，MVL）
3. 调制指数（modulation index，MI）

3.1 锁相值（phase-locking value，PLV）

1. 定义：PLV通过比较两个信号的相位差来评估它们之间的同步性。如果两个信号存在相位-幅度耦合，那么高频时间序列的振幅将在低频振荡。

2. 计算步骤
   

3. 计算公式
   $$PLV=|\frac{{\sum }_{t=1}^n{e}^{i({\theta }_{lt}-{\theta }_{ut})}}{n}|$$
   其中
   n：数据点总数
   t：时间数据点
   θlt：t时刻的低频相位角
   θut：希尔伯特变换后的高频振幅时间序列在t时刻的相位角

4. 优点和缺点：
   优点：受噪声影响小。
   缺点：希尔伯特转换后的振幅时间序列不一定是窄带振荡。

5. python代码实现：

# 利用PLV计算相幅耦合
def PLV_CFC_WSH(phase, amp):
   # phase：低频信号的相位
   # amp：高频信号的振幅
   amp_phi = np.angle(hilbert(amp))        # 高频信号振幅的相位
   observed_plv = np.abs(np.sum(np.exp(1j*(phase - amp_phi)))/len(phase))
   return observed_plv

3.2 平均向量长度（mean vector length，MVL）

1. 定义：MVL是一个向量集合中各个向量长度的平均值，它提供了对向量集合整体大小或幅度的一种度量。在CFC的计算中，将对应时刻的低频信号的相位角作为向量的方向，高频信号的振幅作为向量的幅度。
2. 计算步骤：
   
3. 计算公式：
   $$MVL=\left|\frac{{\sum }_{t=1}^n{a}_t{e}^{i{\theta }_t}}{n}\right|$$

其中:
n：数据点总数
t：时间数据点
at：高频振幅时间序列在t时刻的幅值
θt：低频相位时间序列在t时刻的相位角

1. 优点和缺点：
   优点：对于长数据段，高采样率，高信噪比的信号，建议使用MVL。
   缺点：
   1. MVL的计算受到频率振幅的整体幅度影响。如果某个频率的振幅较大，那么它将在计算MVL时有更大的权重。这可能会导致对整个向量集合大小或幅度的评估存在一定的偏差。
   2. 振幅异常值可能会强烈影响平均向量长度。
   3. 不同向量之间的相位角差异可能会影响平均向量长度的计算结果。
2. python代码实现：

# 利用MVL计算相幅耦合
def MVL_CFC_WSH(phase, amp):
    # phase：低频信号的相位
    # amp：高频信号的振幅
    amp_nor = (amp - np.min(amp)) / (np.max(amp) - np.min(amp))  # 振幅归一化
    observed_mvl = np.abs(np.mean(amp_nor * np.exp(1j*phase)))
    return observed_mvl

3.3 调制指数（modulation index，MI）

1. 定义：MI通过计算幅度和相位的熵来估计幅度调制指数。

2. 计算步骤：
   

3. 计算公式：

   $$\begin{array}{r}p(j)=\frac{\bar{a}}{{\sum }_{k=1}^N{\bar{a}}_k}\end{array}$$

   ā为单个bin的平均振幅，k为bins的运行指数，N为bins的总数；p是一个有N个值的向量；H§为香农熵。
   $$\begin{array}{r}H(p)=-{\sum }_{j=1}^{N}p(j)logp(j)\\ KL(U,X)=logN-H(p)\\ MI=\frac{KL(U,X)}{logN}\end{array}$$

   U为均匀分布，X为数据分布。
   KL(U, K): KL散度描述了用均匀分布U来估计数据的真实分布X的编码损失。

4. 优点和缺点：
   优点：

   1. 归一化熵是一种常用的信号复杂性度量，可以反映信号的不确定性和信息丰富度。
   2. 对于更短的数据，更低的采样率，更高噪声的数据，建议使用MI。
   3. MI是使用最多的CFC方法。

   缺点：
   4. 依赖于特定的熵测量方法；
   5. 对数据预处理的要求较高；

5. python代码实现：

# 利用MI计算相幅耦合
def MI_CFC_WSH(phase, amp, nbins):
    # phase：低频信号的相位
    # amp：高频信号的振幅
    # nbins：划分的箱子数
    bins = np.linspace(-np.pi, np.pi, nbins+1)   
    # 计算每个相位区间内的振幅均值 
    meanAmp = np.zeros(nbins)
    for i in range(nbins):
        indices = np.logical_and(phase >= bins[i], phase < bins[i+1])
        meanAmp[i] = np.mean(amp[indices])
    Pj = meanAmp/np.sum(meanAmp)
    H_Pj = -np.sum((Pj) * np.log(Pj))
    KL = np.log(nbins) - H_Pj
    observed_mi = KL / np.log(nbins)
    return [observed_mi, meanAmp]

4、置换检验

1. 定义：置换检验是一种统计方法，通过对数据进行随机重新排序或重组，产生一组随机分布，从而提供了对观察到的效应进行比较的基准。在计算CFC时，进行置换检验是为了评估观察到的耦合效应是否显著，即是否超出了随机预期。
2. 计算步骤：
   1. 随机在幅值时间序列某数据点上将数据分为两部分，并将两部分时间序列的顺序颠倒，构造出置换幅值时间序列。
   2. 通过计算原始相位时间序列和置换幅值时间序列(反之亦然)之间的耦合值来构造洗牌耦合值。
   3. 这种生成替代数据的方法是最保守的，因为除了所研究的相位角和幅值之间的时间关系，它保留了脑电图数据的所有特征，
   4. 洗牌通常重复200-1000次。
   5. 将观察到的耦合值 $CF{C}_{observed}$ 标准化为洗牌耦合值 $CF{C}_{shuffled}$ 的分布，公式如下：

$$\begin{array}{r}CF{C}_z=\frac{CF{C}_{observed}-\mu CF{C}_{shuffled}}{\sigma CF{C}_{shuffled}}\end{array}$$

其中$CFC$为耦合值，$ \mu $ 为 $CF{C}_{shuffled}$ 均值，$\sigma $ 为 $CF{C}_{shuffled}$ 标准差。只有当 $CF{C}_{observed}$ 大于95%的$CF{C}_{shuffled}$时，才被定义为显著。因为z值等于1.64对应的p值为5%，所以当$CF{C}_z$的值大于1.64 才被定义为显著。

1. python代码实现

npnts = len(eeg)
n_iter = 1000
PLV_shuffled = np.zeros(n_iter)
MVL_shuffled = np.zeros(n_iter)
MI_shuffled = np.zeros(n_iter)
for ni in range(n_iter):
    cutpoint = random.randrange(np.round(npnts/10), np.round(npnts*.9))
    new_amp = np.concatenate((amp[cutpoint:] , amp[:cutpoint]))

    PLV_shuffled[ni] = PLV_CFC_WSH(phase, new_amp)
    MVL_shuffled[ni] = MVL_CFC_WSH(phase, new_amp)
    [MI, meanAmp] = MI_CFC_WSH(phase, new_amp, n_hist_bins)
    MI_shuffled[ni] = MI
    
PLV_CFCz = (PLV_observed - np.mean(PLV_shuffled))/np.std(PLV_shuffled)
MVL_CFCz = (MVL_observed- np.mean(MVL_shuffled))/np.std(MVL_shuffled)
MI_CFCz = (MI_observed - np.mean(MI_shuffled))/np.std(MI_shuffled)

print('PLV_CFCz: ',PLV_CFCz)
print('MVL_CFCz: ',MVL_CFCz)
print('MI_CFCz: ',MI_CFCz)

5、实例分析

5.1 模拟信号的CFC

1. 构建耦合信号

fs = 1000                            # sampling frequency (samples per second/Hz) 
dt = 1/fs                           # seconds per sample 
stopTime = 5                        # length of signal in seconds 
t = np.arange(0, stopTime-dt, dt).T # time vector in seconds 
stopTime_plot = 2                   #limit time axis for improved visualization 

# 生成调制后的信号
f_P = 4       # 低频相位
f_A = 32      # 高频幅度
f_AM = 0.5    # f_AM表示振幅调制的频率
Afp = 1       # AfP是固定的标量，决定SO（慢波）的峰值振幅、
K = 2         # K是固定的标量，决定FO（快波）的峰值振幅、
beta = 0.3

Xf_P = Afp * np.sin(2 * np.pi * f_P * t)       # 慢波（低频相位）
Afa = K * beta * np.sin(2 * np.pi * f_AM * t)  
Xf_A = Afa * np.sin(2 * np.pi * f_A * t)       # 快波（高频振幅）

signal = Xf_A + Xf_P

2. 计算CFC

phasefilt = filterFGx(signal, fs, f_P, 0.5)
phase = np.angle(hilbert(phasefilt))

ampfilt = filterFGx(signal, fs, f_A, 20)
amp = np.abs(hilbert(ampfilt))**2

# 计算观测CFC
PLV_observed = PLV_CFC_WSH(phase, amp)
MVL_observed = MVL_CFC_WSH(phase, amp)
n_hist_bins = 18
[MI_observed, meanAmp] = MI_CFC_WSH(phase, amp, n_hist_bins)

# 置换检验，计算洗牌CFC
npnts = len(signal)
n_iter = 1000
PLV_shuffled = np.zeros(n_iter)
MVL_shuffled = np.zeros(n_iter)
MI_shuffled = np.zeros(n_iter)
for ni in range(n_iter):
    cutpoint = random.randrange(np.round(npnts/10), np.round(npnts*.9))
    new_amp = np.concatenate((amp[cutpoint:] , amp[:cutpoint]))

    PLV_shuffled[ni] = PLV_CFC_WSH(phase, new_amp)
    MVL_shuffled[ni] = MVL_CFC_WSH(phase, new_amp)
    [MI, meanAmp] = MI_CFC_WSH(phase, new_amp, n_hist_bins)
    MI_shuffled[ni] = MI
    
PLV_CFCz = (PLV_observed - np.mean(PLV_shuffled))/np.std(PLV_shuffled)
MVL_CFCz = (MVL_observed- np.mean(MVL_shuffled))/np.std(MVL_shuffled)
MI_CFCz = (MI_observed - np.mean(MI_shuffled))/np.std(MI_shuffled)

print('PLV_CFCz: ',PLV_CFCz)
print('MVL_CFCz: ',MVL_CFCz)
print('MI_CFCz: ',MI_CFCz)

5.2 EEG信号的CFC

1. 加载EEG信号

data = scipy.io.loadmat('accumbens_eeg.mat')
srate = 1000
npnts = len(data)
eeg = data['eeg'].reshape(-1)
fig = plt.figure(figsize=(10, 4),dpi=600)
ax1 = fig.add_subplot(1, 1, 1)  # 子图1
ax1.plot(eeg, label='eeg_data')
ax1.legend()
plt.show()

2. 计算CFC

# define frequencies for phase and for amplitude
phas_freqs =  np.arange(5, 20, 2)
ampl_freqs = np.arange(30, 150, 5)

# number of iterations used for permutation testing
n_iter = 500

# initialize output phase-amplitude matrix
phaseamp = np.zeros([len(phas_freqs),len(ampl_freqs)])

# loop over frequencies for phase
for lower_fi in range(len(phas_freqs)):
    
    # get phase values
    phasefilt = filterFGx(eeg, srate, phas_freqs[lower_fi], phas_freqs[lower_fi]*.4)
    phase = np.angle(hilbert(phasefilt))
    # phase = angle(hilbert(phasefilt))
    
    for upper_fi in range(len(ampl_freqs)):     
        ampfilt = filterFGx(eeg, srate, ampl_freqs[upper_fi], ampl_freqs[upper_fi]*.78)
        amplit = np.abs(hilbert(ampfilt))**2
        
        # calculate observed modulation index(MVL)
        [modidx, _] = MI_CFC_WSH(phase, amplit, 30)
        
        # now use permutation testing to get Z-value
        bm = np.zeros(n_iter)
        for bi in range (n_iter):
            cutpoint = random.randrange(np.round(npnts/10), np.round(npnts*.9))
            new_amp = np.concatenate((amplit[cutpoint:] , amplit[:cutpoint]))
            [mi, _] = MI_CFC_WSH(phase, new_amp, 30)
            bm[bi] = mi
        # the value we use is the normalized distance away from the mean of
        # boot-strapped values
        phaseamp[lower_fi, upper_fi] = (modidx - np.mean(bm)) / np.std(bm)

3. 可视化

fig, ax = plt.subplots()
CFC_WSH = ax.contourf(phas_freqs, ampl_freqs, phaseamp.T)
cbar = fig.colorbar(CFC_WSH)
plt.show()

完整代码：

https://github.com/weisihong9/Cross_Frequency_Coupling.git

参考文献

[1] Li, C., et al., 2021. High-frequency hubs of the ictal cross-frequency coupling network predict surgical outcome in epilepsy patients. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 29, pp.1290-1299. (2区)
[2] Salimpour, Y. and Anderson, W.S., 2019. Cross-frequency coupling based neuromodulation for treating neurological disorders. Frontiers in neuroscience, 13, p.125. (3区)
[3] Musaeus, C.S., et al., 2020. Electroencephalographic cross-frequency coupling as a sign of disease progression in patients with mild cognitive impairment: a pilot study. Frontiers in neuroscience, 14, p.790. (3区)
[4] Zhang, W., et al., 2023. Altered fronto-central theta-gamma coupling in major depressive disorder during auditory steady-state responses. Clinical Neurophysiology, 146, pp.65-76. (3区)
[5] Hülsemann, Mareike J., Ewald Naumann and Björn Rasch. 《Quantification of Phase-Amplitude Coupling in Neuronal Oscillations: Comparison of Phase-Locking Value, Mean Vector Length, Modulation Index, and Generalized-Linear-Modeling-Cross-Frequency-Coupling》. Frontiers in Neuroscience (2019). (3区)