神経ネットワークの自己組織化臨界現象を調べる上で、global lability of synchronization(同期のグローバルな不安定性)法という手法を使うことになったので、その計算方法をメモしておく。
大まかにいうと、ある神経ネットワークの活動において時刻$t$で同期しているニューロンのペアの個数$M(t)$と時刻$t+1$で同期しているペアの個数$M(t+1)$がどれだけ違うか、ということを評価する。値が大きければ大きいほど全体の同期/非同期プロセスを表し、小さいほど局所的なものを表す。
位相同期の評価
時系列データとして、二つの信号$s_i(t)$と$s_j(t)$の二種類が与えられたとする。
この二つの信号がある時刻$t$において
- 位相差$\theta_{i, j}(t)$の絶対値が$\pi/4$未満である
- 時間窓$v$の範囲での位相差$\theta_{i, j}(t)$に対する同期パラメータ$\gamma_{i,j}(t)$が$\sqrt{1/2}$より大きいこと
この二つの条件を満たす時、二つの信号は同期しているとする。
ここで位相差$\theta_{i, j}(t)$は以下のように計算を行う。
$$ \Delta \theta_{i,j}(t)=\arctan[{\tilde s_i(t)s_j(t)- \tilde s_j(t)s_i(t)\over s_i(t)s_j(t)- \tilde s_j(t)\tilde s_i(t)}] $$
ここで$\tilde s(t)$は信号$s(t)$のヒルベルト変換を表す。
$$ \tilde s(t) = {1\over \pi}\int^\infty_{-\infty} {s(\tau)\over t-\tau} d\tau$$
ヒルベルト変換って何？という方は以下↓
参考: https://www.onosokki.co.jp/HP-WK/eMM_back/emm180.pdf
ざっと言うと実部のみの信号$s(t)$が本来含む虚部の信号$\tilde s(t)$を返す操作。さらに大雑把に言えば位相を90°ずらす操作に当たるらしい。
同期パラメータ$\gamma_{i,j}(t)$は以下のように求める。
$$ \gamma_{i,j}(t) = \sqrt{[{1\over v}\sum^{t+v}_t \cos \Delta \theta _{i,j}(t)]^2+[{1\over v}\sum^{t+v}_t \sin\Delta \theta _{i,j}(t)]} $$
ここで$v$は時間窓を表す。
ネットワークの中の全てのニューロン信号について、これらの条件を満たす一致するペアの個数を$M(t)$とする。以下のような式で表すことができる。
$$ M(t)=\sum[|\Delta \theta_{i,j}|<\pi/4 ~and ~ \gamma_{i,j}(t)>\sqrt{1/2}] $$
不安定性の評価
得られたペアの個数$M(t)$と$M(t+1)$から、以下のようにして不安定性$l(t)$を求めることができる。
$$ l(t) = |M(t+1)-M(t)|^2 $$
実装
前提として、ヒルベルト変換は以下のように計算することができる。
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import hilbert
x = np.arange(-10,10, 0.1)
y = np.cos(x)
z = hilbert(y)
plt.plot(x,z.real, label="hil_real")
plt.plot(x,z.imag, label="hil_imag")
plt.legend()
plt.show()
元の入力であるcos波に対して虚部であるsin波が取得できていることがわかる。
これを用いて実装を行ったのが以下のクラス。
import itertools
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import hilbert
class DLS():
def __init__(self, signals, v=50):
self.signals = [hilbert(s) for s in signals]
self.v = v #時間窓を50ステップとする
def theta(self, s1, s2, t_i):
y = s1[t_i].imag*s2[t_i].real - s2[t_i].imag*s1[t_i].real
x = s1[t_i].real*s2[t_i].real - s1[t_i].imag*s2[t_i].imag
return np.arctan(y/x)
def gamma(self, i, j, t_i):
if t_i + self.v >= len(self.signals[i]):
return 0
T_i = np.arange(t_i, t_i+self.v)
cos_ave = np.cos(self.theta(self.signals[i], self.signals[j], T_i)).sum()/len(T_i)
sin_ave = np.sin(self.theta(self.signals[i], self.signals[j], T_i)).sum()/len(T_i)
return np.sqrt(cos_ave**2 + sin_ave**2)
def isPair(self, i, j, t_i):
theta_abs = abs(self.theta(self.signals[i], self.signals[j], t_i))
gamma = self.gamma(i,j,t_i)
if theta_abs < np.pi/4 and gamma > np.sqrt(1/2):
return 1
else:
return 0
def M(self, t_i):
count = 0
for i, j in itertools.combinations(range(len(self.signals)), 2):
count += self.isPair(i, j, t_i)
return count
def l(self, t_i):
return (self.M(t_i+1)-self.M(t_i))**2
試しに全く同じsin波を3つ入れて計算してみる。
v = 50
x = np.sin(np.arange(0,10, 0.1))
signals = [x, x, x]
dls = DLS(signals, v)
for i, x in enumerate(x[:v]):
print(dls.M(i))
# 一致する組み合わせ数3が出力される
print(dls.l(i))
# 完全に安定なので0が出力される
こうして計算を行うことができた。
参考文献
Kitzbichler, M. G., Smith, M. L., Christensen, S. R., & Bullmore, E. (2009). Broadband criticality of human brain network synchronization. PLoS computational biology, 5(3), e1000314.
Aguilar-Velázquez, D., & Guzmán-Vargas, L. (2019). Critical synchronization and 1/f noise in inhibitory/excitatory rich-club neural networks. Scientific reports, 9(1), 1-13.