神経ネットワークの振動パターンの解析
2008MI256 角田由紀
指導教員：大石泰章

1と書く．ここで M の (i, j) 成分 mij は j 番目の神経細胞

はじめに線虫やクリオネなどの生物の遊泳の動作は CPG (Cen-tral Pattern Generator) と呼ばれる神経ネットワークによって制御されている．CPG とは，神経細胞によって形成から i 番目の神経細胞までの結合関係を表す．mij が負な
らば結合は抑制性であり，正ならば結合は興奮性である．
以上をまとめると CPG の方程式は，

q = F (s)M Ψ(q)された非線形な発振器のことである．ある神経細胞が興奮するとその細胞は電気信号を出し，隣接する細胞との結合が抑制性であれば興奮状態を抑制し，結合が興奮性であれば興奮状態を励起する．このように CPG では神経細胞がとなる．

q が 1 次元の場合，式 (3) が周波数 ω の周期解を持つとしてフーリエ級数展開すると，作用し合うことで身体を動かすタイミングを生成し，その結果として身体の動きが生じる．この CPG を利用して，(3)q(t) =生物の複雑な動きをロボットに導入できないかという研究∞∑

ak sin(ωkt) + bk cos(ωkt)

k=0が進んでいる [2]．Iwasaki [1] は同じ動特性を持つ神経細を得る．ω が F (s)M の帯域幅と同程度であると考え，胞の集まりによって CPG をモデル化し，その CPG の解k ≥ 2 の成分を無視する．さらに b0 = 0 とすると，析と設計を行う系統的方法を開発した．しかし Iwasaki は近似を多く利用しているため，近似をすることが妥当であq(t) ∼= a1 sin(ωt) + b1 cos(ωt)= α sin(ωt + ϕ)るのか検証をする必要がある．
本研究では Iwasaki の方法に基づいて CPG の解析問題を解き，その方法の検証を行う．具体的には，Iwasaki の方法によって振動パターン (周波数，振幅，位相) を推定すとなる．q が N 次元の場合も，各成分について同様に考え，qi (t) ∼= αi sin(ωt + ϕi )(i = 1，. . . ，n)(4)

る．また，推定された振動パターンと実際の振動パターンと す る ．こ こ で αi は 振 幅 ，ϕi は 位 相 で あ る ．ま た ，を比較し，近似の妥当性を検証する．

tanh(qi ) は記述関数 κi によって近似し，2tanh(qi ) ∼= κi (αi )qi

CPG の方程式と振動パターンの推定

神経細胞が 1 個の場合を考える．神経細胞の動特性を， とする．これらの近似から CPG の方程式は，電気信号の入力が u，出力が v のとき，v = tanh(q)，q = f (s)u(F (jω)M K(α) − I)p = 0(1)pi := αi ejϕi ，(5)K(α) := diag(κi (αi ))とモデル化する．このとき tanh(q) は静的な非線形関数でとなる．ここで p は q を ejωt で割ったものである．このあり，f (s) は伝達関数関係から近似的に条件を満たす周波数，振幅，位相が求められ振動パターンを推定することができる．

ω0f (s) =s + ω0

qi の振動の振幅 αi が互いに等しくなるように CPG はとする．次に神経細胞が N 個ある場合を考えると式 (1)ある対称性を持つと仮定する．このときすべての i に対しは次のようになる．て αi = α とすると，式 (5) は，

vi = tanh(qi )，qi = f (s)ui(M − λI)p = 0(i = 1, 2, . . . , N )以 下 N 個 の 信 号 を ベ ク ト ル 表 現 し て ，v==(v1 v2 . . . vN )T ，q = (q1 q2 . . . qN )T ，u(u1 u2 . . . uN )T と書く．また，神経細胞の入出力関係もまとめて，v = Ψ(q)，q = F (s)u(2)pi := ejϕi ，λ :=(6)1
κ(α)f (jω)となる．ここで λ は M の固有値，p はそれに対応する固有ベクトルである．固有ベクトル p から位相 ϕi を推定し，固有値 λ から次の式を解くことによって周波数 ω と振幅α が推定できる．のように書く．ただし，Ψ は各成分の tanh(qi ) を計算することを表し，F (s) は f (s)I を表す．一方，N 個の神経細胞間の相互結合を行列 M を使って，u = Mvλκ(α) = 1 + j

ω
ω0(7)

このようにして行列 M に対する振動パターン (周波数，振幅，位相) を推定する．

3推定の実際とその評価

3.1

推定結果

321

20
−1
−10012002−1 0−1 −111020q(t)行列 M を



01

1
0
20-1-2-3
101214161820時間

図 2 シミュレーション (ω0 = 3)とし，ω0 = 1 とする推定を行う．行列 M の固有値は2 ± 1.90j ，2，2 ± 1.17j である．ここでは最大固有値である λ = 2 + 1.90j について考える．最大固有値とは，最も2010パワー/周波数 (dB/Hz)

大きい実部を持つ固有値の中で，さらに最も大きい虚部を持つ固有値のことである．式 (6)，(7) より ω と ϕ を求めると，ω = 0.1514 (Hz)，ϕ = 0，36，72，108，144 (° ) と推定できる．

0-10-20-30-40

次に，式 (3) を逆ラプラス変換した微分方程式を解きシ-50ミュレーションを行う．tanh(q) が非線形関数のためルン

0246

周波数 (Hz)

ゲクッタ法を用いて微分方程式を解くと，図 1 のようなグ
図 3 スペクトル密度ラフが得られる．
3クトル密度で表したものである．図 3 を見ると最初のピー2ク時の値に対して 3 倍，5 倍の値で周波数のスペクトル密q(t)1度が高くなっている．このことから第 3 次，第 5 次の高調0波の成分が少なからず含まれてしまうため，k ≥ 2 の成分-1を無視することには疑問が残る．しかし，Iwasaki の方法-2による推定値とシミュレーション結果の位相の比較は表 1-310のようになり，依然として良く一致している．このことか

1520時間2530

ら，ω0 の値によらず q のフーリエ級数展開の k ≥ 2 の成

図 1 シミュレーションによる CPG の振動パターン分を無視する Iwasaki の方法が有効であることがわかる．
初期値の影響を考えるため，初期値を変更してシミュ

式 (6)，(7) によって推定された値とシミュレーションでレーションを行う．その結果，変更前と同じ周期と歪みを得られた値を比較する．表 1 は推定とシミュレーションで持ったグラフが得られた．よって，初期値の影響は受けな得られた周期と ϕi の値を比較したものである．ω0 = 3 といことがわかる．

は ω0 の値を変更したときの結果である．これを見ると周4

表 1 推定とシミュレーションの値の比較
周期推定値
シミュ

ω0 = 36.61
7.242.61

本研究では，Isasaki の方法に基づいて CPG の方程式を解析問題について解き，その妥当性の検証を行った．近似

ϕ10.00.00.0

ϕ236.036.136.0

ϕ372.072.0
72.0

ϕ4108.0
108.1
108.0

ϕ5144.0
144.0
144.0

期は異なるが，位相は推定値とシミュレーションで得られた値はほぼ合致することがわかる．

3.2

おわりに近似の妥当性

フーリエ級数展開で近似を行うとき k ≥ 2 の成分を無視しても良いかを確かめるために，ω0 の値を変更してシミュレーションを行う．ω0 = 3 に変更しシミュレーションを行うと図 2 のようなグラフが得られる．図 3 は図 2 をスペを用いた CPG の方程式で解析問題を解いてもシミュレーションに近い結果が得られ，周波数 ω0 を変えても同様であったことから，Iwasaki の方法は妥当であると言える．

参考文献[1] T. Iwasaki: Multivariable harmonic balance for central pattern generators. Automatica, vol. 44 (2008),
pp. 3061–3069.
[2] 小田島正・湯浅秀男・伊藤正美：「均質モジュールから成る自律分散型多脚歩行ロボットシステム−機能とアルゴリズム−」．日本ロボット学会誌，第 16 巻第 7 号(1998)，pp. 957–964.