Pythonで実装するEEG解析演習#1 ERPを計算し、トポグラフィカルマップで可視化する

はじめに

 本記事は、脳波解析の初学者向けの記事です。本記事では、Mike X Cohen. ANALYZING NEURAL TIME SERIES DATA. Massachusetts: The MIT Press, 2014. に掲載されている脳波解析の演習問題をPythonで実装し、解説します。Pythonに関する基礎知識（計算やグラフ作成など）や書籍が手元にあると、本記事の理解が深まります。書籍が手元にない方でも、記事の内容を理解できるようにEEG解析に関する概要や問題文を記載しております。演習問題で利用するEEGのデータセットは、書籍で利用されているデータセットを使用します（誰でもダウンロード可）。なお、記事掲載にあたり、事前に著者のMike X Cohen様から許諾を得ております。

 第一回目は、書籍の第9章 Overview of Time-Domain EEG Analysesの内容をもとにERPに関する演習問題を実装します。

ERPとは

 ERP: Event Related Potentials（事象関連電位）は、刺激やイベントに関連した脳の電気的応答をEEG信号から抽出・可視化したものです。時間軸でイベントや刺激に同期させて、EEG信号を平均することでERPを計算できます。計算結果で得られたERPの特徴から、P300やN400といったERP成分を観察でき、この成分の振る舞いから脳の認知処理などを推測します。PやNは、それぞれPositiveとNegativeを表しており、信号の振れの極性（陽性・陰性）を表しています。また、PやNに続く数字は、その成分が観察できる時点を示しています。

ERPの主なメリットです。

1. 計算が簡単で高速（試行ごとの信号を足して平均するだけ）
2. 時間的な精度が高い（ミリ秒単位）
3. 被験者の数が少なくても、データを迅速に品質チェックできるので実験設計の正しさの確認が容易

また、ERPはデメリットもあります。

1. 計算結果を神経生理学的メカニズムに結びつけるのが難しい
   これは、ERPが空間的精度が低く、脳のどの部分の活動がERPと対応しているかを特定することが難しいためです。また、試行ごとに位相の揃わない成分は、イベントや刺激に関連していても平均化の過程で相殺されてしまう可能性があります。
2. S/N比の高いERPを取得するためには、数十〜数百の試行が必要
   これにより、被験者の疲労や集中力の低下が、実験結果に影響する可能性があります。

トポグラフィカルマップ

 EEG解析におけるトポグラフィカルマップは、電極ごとの電位分布を頭部表面にマッピングしたものです。下図は、電極ごとのERPをトポグラフィカルマップで示したものです。トポグラフィカルマップは、EEG解析における基本的な可視化手段の一つで、脳活動の空間的なパターンを視覚的に理解するために有効な方法です。

演習問題

 各電極におけるERPを計算し、トポグラフィカルマップに表示する5つの時点を選択せよ。例えば 0~400msまで100msステップずつなど。1つの図にこれらの時点におけるトポグラフィカルマップを作成すること。
また、S/N比を向上させるため、各時点の20ms前から20ms後までの脳活動の平均値を表示すること。例えば、時点が200msの場合、180ms~220msまでの脳活動の平均値を表すこと。
各トポグラフィカルマップのタイトルには、中心となる時点を記載すること。

以下、原文 (Mike X Cohen. ANALYZING NEURAL TIME SERIES DATA. Massachusetts: The MIT Press, 2014, p106-107)

Compute the ERP at each electrode. Select five time points at which to show topographical plots(e.g., 0 to 400 ms in 100-ms steps). In one figure, make a series of topographical plots at these time points. To increase the signal-to-noise ratio, make each plot show the average of activity from 20ms before until 20 ms after each time point. For example, the topographical plot from 200 ms should show average activity from 180 ms until 220 ms.Indicate the center time point in a title on each subplot.

実行環境

 本記事では、PythonでEEG解析をするためにMNEを使用します。MNEは、PythonでEEGやMEGといった生体信号を処理・分析・可視化するためのオープンソースパッケージです。

Python 3.12.0
matplotlib==3.8.2
mne==1.11.0
numpy==1.26.3
scipy==1.15.3

実装

 ERPの計算は、numpyのmean関数を使います。今回使用するEEGのデータ形式は(チャンネル, 時間, 試行)です。チャンネルごとにERPを計算するので、mean関数の引数にaxis=(1,2)を入れています。
トポグラフィカルマップの作成は、mneのviz.plot_topomap関数を使用します。viz.plot_topomap関数内でカラーバーを表示する方法が見当たらなかったので、matplotlibのcolorbar関数を使います。

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat # matファイルの読み込み
import mne

## EEGデータの読み込み
EEG = loadmat('sampleEEGdata.mat')['EEG'][0, 0]
eeg = EEG['data'][:, :, :] # eeg.shape→(64, 640, 99)
sfreq = EEG['srate'][0] # SamplingRate: 256

## 作成する5つの時点を選択(0ms, 100ms, 200ms, 300ms, 400ms)
time_points = np.arange(0, 401, 100)
time_window = 20 # 20ms

## mne用のチャンネル座標情報を作る
info = mne.create_info(
    ch_names = [ch[0][0] for ch in EEG['chanlocs'][0]],
    sfreq = sfreq,
    ch_types = "eeg"
)

montage = mne.channels.make_standard_montage("standard_1020")
info.set_montage(montage)

## 描画、ERP計算
fig, axes= plt.subplots(1, 5, figsize=(18, 6), constrained_layout=True)
## 各時点の前後20msのインデックスを取得
for i, tp in enumerate(time_points):
    start_time = tp - time_window
    end_time = tp + time_window
    start_index = np.argmin(np.abs(EEG['times'][0] - start_time))
    end_index = np.argmin(np.abs(EEG['times'][0] - end_time))
    
    ## 各時点のERP計算
    erp = np.mean(eeg[:, start_index : end_index + 1, :], axis=(1,2))

    ## トポグラフィカルマップ作成
    im, _ = mne.viz.plot_topomap(
        erp,
        info,
        axes = axes[i],
        cmap = 'RdBu_r', # カラーマップを指定
        show = False,
        contours = 6, # 等高線の数を指定
        vlim = (-10, 10) # スケール調整
        )
    axes[i].set_title(str(tp) + "ms")

## カラーバーの表示
cbar = fig.colorbar(im, ax=axes[-1], shrink=0.7)
cbar.set_label('Amplitude (uV)')

plt.show()

実行結果

 0〜400ミリ秒までの5つの時点のERPをトポグラフィカルマップに可視化できました。