目次
Math
Analysis
Experiment
Simulations

05脳波の位相同期解析

周波数の解析はフーリエやウェーブレットを使って今までにも説明をしてきましたが，位相に関する話はほぼありませんでした．ここでは，脳波の最後の次元，位相を使った解析手法に触れていきます．

▶︎
all
running...

はじめに

位相とは，の説明は大丈夫だと思いますが，「脳波の位相」が何を表しているのかをまず考えていきます．

脳波の位相，はより正確には周波数分解した脳波のそれぞれの周波数についての位相，です．取り出すには複素数である必要があるため，まずウェーブレットやヒルベルト変換を用いて複素信号にします．

そしたら普通にangle(x)のようにして瞬時位相の情報を取り出します．

次に脳波位相の考え方です．

物理的には，それぞれの周波数を持った振動子として捉えられる神経細胞(群)ですが，こいつらは積分発火型の振動子です．つまり発火する直前くらいのタイミングで受け取った信号に対する反応性が高く，逆に発火した後の不応期とかだと低くなります(シナプス前細胞Ａからの信号だけじゃ発火しないし，発火するのはＡからの信号をもらったずっと後だからあまり影響がないみたいな)．この特性ゆえ，発火のタイミング，つまり位相が近い振動子同士ほど情報伝達がちゃんと起きるし，可塑性によって強化されてもいきます．

言い換えると，同期的に活動している神経群同士は互いに影響を及ぼしやすく，それによって同期も強まるという性質があると考えられます．

この考え方は，Fries., 2005なんかで言われていることですね．さらに，Okun., et al. 2008などで言われるような興奮性/抑制性の均衡のことなんかを考慮すると，尤もらしい説明に聞こえます．

この前提ゆえに，異なる領域の脳波で同期が見られたりした場合，その領域間でnetworkがある，情報のやりとりが行われている，などと考えるのが位相同期解析の基本的な方針です．

あるいは，同じ情報を処理しているのだったら同調したように活動しているはずだよね，という考え方も取られます．が，個人的には上述のような説明の方が妥当というか説得力があるような気がします．

ともあれ，さっそく解析を見ていきましょう．

試行間位相同期

位相同期は，何と何の位相が同期しているのかで種類が分けられます．まずは基本の，特定電極の試行間での位相同期を見ます．ERPの位相版って感じです．

毎回同じような反応を示す，つまり毎試行位相がそろっているようなものであれば良い結果，そうでなければ微妙な結果が出る解析です．名前は厄介なことに何通りもあって，全部同じ事なのにみんな好き勝手呼んでいます．

Phase-locking Index, Phase-locking factor, Inter trial coherence, Inter trial phase coherenceなどがあります．やめてほしい．そういいつつ筆者も気分で言い分けてる気がします．ごめんなさい．式は

P(t)=∣∑k=1Kexp(iϕkf(t))K∣P(t) = |\frac{\sum_{k=1}^K exp(i \phi_k^f(t))}{K}|

のようになります．ERP同様，kkは試行数，ttはタイムポイントで，ffは周波数，iiは虚数単位，それからϕ\phiは瞬時位相になります．瞬時位相はヒルベルト変換などで解析信号をだし，そいつの位相時系列のうちタイムポイントttに対応する値です．

簡単ですね．ERPの時は波形x(t)x(t)そのものの加算平均をしていたけど，今回は波形のうちの瞬時位相を加算平均しているわけです．

結果は0--1の値を取り，試行間で位相がそろっているほど1に近付き，そうでないと0に近付きます．複素数というか極座標表示なので，加算平均しないで一つのトライアルのみで絶対値を取ると当然1です．

これが加算平均しても1に近い値にあるのなら同期しているよね，とする解析です．ただ，値自体は試行数によっても変動するので大小の比較とかは難しく，差があるかどうかは数値でなく統計的議論をする必要があります．

上の図が概念図です．それぞれの試行での位相が極座標で灰色で示されています．これらの加算平均したもの，位相同期度が黒の太線で描かれています．左の図は，試行間である程度位相が揃っているので，黒線が長い，つまり同期度が高くなっている一方，右の図はばらばらになっているため，同期度が低くなっています．

特性として，計算の際に振幅は考慮しないために位相のみに影響され，振幅が大きかったり小さかったりしても影響を受けないという点があります．expの前に振幅はついていませんよね．

だからこそ綺麗に同期していれば1になります．これゆえ，α波などの大きな波に埋もれた小さなγ波の位相が同期しているような状態でもしっかりと検出する事が出来ます．

ただ解釈は(個人的には)ちょっと難しい気がして，ERP同様に「だからなに？」になりかねないというか，位相が試行間で揃うこと自体が重要だとするようなうまい説明が求められるような気もします．

論文としては位相同期を見るのはここ最近でもまだトレンド的なものではありますし，学士の卒論とかでここまでやってれば結構いい感じなんじゃないですかね．

部位間位相同期

こっちの方が個人的にはしっくりきます．部位間での位相同期です．基本的には，試行間位相同期の発展的な式なのですぐ理解できるかと思います．名前はPhase locking valueとか言ったりしますが，これも論文によってみんな名前がてきとうすぎて，部位間なのか試行間なのか分かりにくいです．

筆者のラボはそれゆえ，Phase Synchlony Index(PSI)と称しています．最近筆者も慣れてきましたが，これならSynchronyなので部位間同期の方だな，と分かりやすいので確かに便利です．

式はこう

I(t)=∣∑k=1Kexp(iϕk1f(t)−ϕk2f(t))K∣I(t) = |\frac{\sum_{k=1}^K exp(i \phi_{k1}^f(t) -\phi_{k2}^f(t) )}{K}|

ここで，ϕ\phiの下につくkの数字は電極1,2を意味します．二つの電極(or 領域)の周波数fの瞬時位相を算出し，そいつらの位相差を加算平均します．やっぱり位相差がそろっていれば1，ばらばらなら0の値に近付いていきます．

ここで少し面白いのは，位相差の加算平均なので位相差自体がどれほどの大きさであるかは問題ない点です．神経系も物理・化学的なワイヤを使って情報通信をしているようなものなので当然伝達のラグはあります．それ以外にも積分発火だから閾値に達するまでのラグもあるし．ワイヤの伝達時間は距離にもよりますし，何が言いたいかというと離れた領域の細胞が100%ラグなく同期するなんてことはあり得ません．(余談ですが非線形力学的にも同期している振動子は若干の位相ずれが必ず存在します)．

なのでどんな位相差であっても問題なく，とにかく位相差が一定かどうかを検出する手法になります．

Phase-lag Index

Cross Frequency Coupling

Cross Frequency Couplingは，異なる周波数成分同士の同期を考える手法です．

Phase Phase Coupling

二つの振動子を仮定して，高周波と低周波の位相同期を考える手法です．言うまでもなく，周波数が違うので単純に位相を比較しても高周波がどんどん引き離していく事になってしまいます．なので比較のためにはちょっと工夫が必要で，いくつかやり方があるようです．一般にCross Frequency Couplingといったらこいつを指すことが多いかもです．

まず一つ目のやり方は基本的に，高周波が低周波の倍数であるときのみに成り立つ計算です．低周波の位相を高周波との比の分進めてあげる事で，疑似的に低周波の周波数を引き上げます．その時の位相差がどうなっているかを見るのが手法です．

PPC(t)=∣∑k=1Kexp(i(mϕk1flow(t)−ϕk2fhigh(t)))K∣where.flow:fhigh=1:mPPC(t) = |\frac{\sum_{k=1}^K exp(i( m \phi_{k1}^{f_{low}}(t) -\phi_{k2}^{f_{high}}(t) ))}{K}| \\ where. f_{low}:f_{high} = 1:m \nonumber

もう一つのやり方は，低周波高周波それぞれをバンドパスによって抜き出した上で，更に高周波成分に対し同様の低周波でバンドパスをかけ，現れた低周波を比較する手法です．

PPC(t)=∣∑k=1Kexp(i(ϕlow(t)−ϕAhigh(t)))K∣PPC(t) = |\frac{\sum_{k=1}^K exp(i(\phi_{low}(t) - \phi_{A_{high}}(t)))}{K}|

Phase Amplitude Coupling

低周波の位相と，高周波の振幅位相を比較する手法です．低周波の位相に依存して高周波の活動がmodulateされているとか，そういった考え方が出来ます．一般に低周脳波は遠距離通信的な役割を持っているのではないかなどと言われているため，たとえばトップダウンの信号を低周波がもってきて，その情報量(位相)に依存してローカルな高周波活動，たとえば知覚処理の動態が変化しているのではないか，なんて議論に使われていたりします．

PAC(t)=∣∑k=1KAhigh(t)exp(iϕlow(t))K∣PAC(t) = |\frac{\sum_{k=1}^K A_{high}(t) exp(i\phi_{low}(t))}{K}|\\

ここで，AhighA_{high}は高周波成分の瞬時振幅，つまりエンベロープを取るやつです．よって，バンドパスフィルタによって取り出した高周波成分をヒルベルト変換し，絶対値を取る事で求められます．ϕlow\phi_{low}は低周波成分の瞬時位相です．こちらも低周波をバンドパスで抜き出し，そいつの位相を取り出す事で求められます．

プログラム 
%% Cross Frequency Coupling
%PAC
clear;
clc;
srate=2000;
t = (-pi:1/srate:pi)';
n = length(t);
f_low=5;
f_high = 120;
x1 = ((2.0 + sin(2*pi*f_low*t)) .* sin(2*pi*f_high*t))';
x2 = ((2.0 + sin(2*pi*f_low*t)) .* sin(2*pi*f_high*t)+randn(length(t),1))';

% high-gamma(80-150)Hz amplitude time series
amp_high1 = abs(hilbert(eegfilt(x1, srate, 80, 150)));
amp_high2 = abs(hilbert(eegfilt(x2, srate, 80, 150)));

% theta (4-8)Hz  phase time series
phase_low1 = angle(hilbert(eegfilt(x1, srate, 4, 6)));
phase_low2 = angle(hilbert(eegfilt(x2, srate, 4, 6)));

% phase time series of envelope of high_gamma
phase_env1 = angle(hilbert(eegfilt(amp_high1, srate, 4, 6)));
phase_env2 = angle(hilbert(eegfilt(amp_high2, srate, 4, 6)));

% comupute PAC
pac1 = abs(mean(amp_high1.*exp(1i*phase_low1)));
pac2 = abs(mean(amp_high2.*exp(1i*phase_low2)));

figure;
subplot(2,2,1)
plot(x1)
hold on
plot(amp_high1)
xlim([0,1000])
legend('raw data','high-amp(120Hz)')
title('Phase Amplitude Coupling', 'FontSize', 12, 'FontName', 'Arial')
hold off

subplot(2,2,2)
plot(x2)
hold on
plot(amp_high2)
xlim([0,1000])
legend('raw data','high-amp(120Hz)')
title('Phase Amplitude Coupling with noise',  'FontSize', 12, 'FontName', 'Arial')
hold off

subplot(2,2,3)
plot(phase_low1);
hold on
plot(phase_env1);
xlim([0 1000])
title('Phase of low and high-amplitude envelope', 'FontSize', 12, 'FontName', 'Arial');
legend('Low-freq', 'high-greq-Envelope');
hold off

subplot(2,2,4)
plot(phase_low2);
hold on
plot(phase_env2);
xlim([0 1000])
title('Phase of low and high-amplitude envelope',  'FontSize', 12, 'FontName', 'Arial');
legend('Low-freq', 'high-freq-Envelope');
hold off

Modulation Index

Modulation Indexは，Phase Amplitude Couplingがどれだけ起きているかを統計的に議論するためのメソッドです．今度まとめます．

部位間同期によるnetwork解析

部位間の同期をPSIなどを使って計算した場合，それぞれの周波数帯域での全脳のコネクティビティマップのようなものを算出することが可能です．

これを用いて，特定の部位間の同期だけでなく脳全体でどのようなネットワークが形成されていたか，またそのネットワークがどのように変化するかを考えることも可能です．

static Functional Connectivity

dynamic Functional Connectivity

【目次へ】