Simulink による反応拡散系 CNN のシミュレーション
宮崎大学工学部教育研究支援技術センター
小島丈英
はじめに
Simulink は、MATLAB のツールボックスの一つとして力学システムをブロック図としてモデル化することでシミュレーションを行い、その結果を分析するソフトとして世界中で広く使われている。MATLAB を用いてシミュレーションする際は微分方程式を算出した後、その解法プログラムとしての「ODE45」を用いてプログラムを作成した上でシミュレーションを実行しなければならず、MATLAB のプログラムを作成する際にMATLAB に対するある程度の知識が必要となってくる。しかし、Simulink を用いる場合は、微分方程式を算出したら、マウスを用いてその微分方程式に基づいたブロック図を作成し、シミュレーションのパラメータを入力することでシミュレーションを行うことができ、ブロック図を作ることに慣れさえすれば非常に簡単にシミュレーションを行うことができる。
本稿では、電気電子工学科回路システム研究室に導入されている simulink を用いて、32 個のセルが一次元空間上に配列された反応拡散 CNN(Cellular Neural Network)[1]を構築し、traveling wave, trigger wave のシミュレーションを行う。最初に、simulink の立ち上げ、電気回路への簡単な応用例について説明した後、反応拡散CNN のシミュレーションについて述べる。
キーワード：MATLAB、Simulink、RL 回路、ＣＮＮ、反応拡散系 CNN[1]1.Simulink の起動方法
まず、簡単に Simulink の起動方法について説明する。起動方法は MATLAB 上で「Simulink」と入力するか、をクリックする。すると、図1 のような画面が出てくるので図１左上の○部分をクリックしてブロック図を新規作成するか、ツールバーから｢ファイル→新規作成→モデル｣とクリックしていくと作ることができる。
次に Simulink で使用するライブラリについては、口頭発表にて説明する。
2.電気回路への簡単な応用例
Simulink を用いたシミュレーションの簡単な例として図２の回路[2]を考える。この回路図について節点方程式を立てて、整理すると次式のようになる。
1⎧ di L ( t )[E − R 0 {i L ( t ) + i R ( t ) }]=⎪⎪dt
L⎨
⎪ i R (t ) = L di L ( t )⎪⎩R1dt上式を用いて図 3 のように、ブロック図を描いた後、パラメータとシミュレーション時間を入力して実行すると、図４のように抵抗とコイルに流れる電流の時間波形を得ることができる。
なお、図５のように Matlab のワークスペースへ出力するブロックも Simulink にはあり、そのブロックを用いることで Matlab 上の Figure コマンドを用いてグラフの編集をすることができる。

図１ Simulink のスタート画面

図２ 簡単な RL 回路[2]図４(上図) 抵抗 R1 に流れる電流の時間波形(下図) L に流れる電流の時間波形

図３

RL 回路のブロック図 例１

このように、数式で表すことができるものはSimulink を用いて、シミュレーションすることができる。Simulink に最初から入っているデモには、数学や物理を用いた一般のアプリケーションのほかにも自動車の制御システムや航空宇宙アプリケーションとして、航空機などに用いられる制御システムなどのシミュレーション例がある。
3.Simulink を用いた反応拡散系 CNN のシミュレーション
本研究では、Simulink を用いて反応拡散系 CNNに関する論文[1]の結果をシミュレーションし、その妥当性を検討する。
3.1 CNN の概要
まず、本論に入る前に、CNN の概要について述べる。CNN は、脳科学や神経回路網などで用いられる Cellular Neural Network の意味と結合振動
力 学 系 な ど で 用 い ら れ る Cellular Nonlinear
Network の意味の二つに分けられている。今回の応用例では、CNN を Cellular Neural Network の意味で用いる。

図５

RL 回路のブロック図

例２

A. 孤立セル
CNN の一つのセルは、一般に入力 u ij 、しきい値

z ij 、状態変数 xij 、出力 y ij の 4 変数を用いて表される。
B. CNN セルの近傍圏
M 行 N 列の 2 次元配列で構成される CNN を考える。任意の i 行 j 列に位置する CNN セル C ij は近傍セルにのみ局所的に結合され、半径 r の近傍内で直接影響する範囲 S ij (r ) は次式で定義される。

Sij = {Ckl : max ( k −i , l − j ) ≤ r ,1≤ k ≤ M ,1≤ l ≤ N}

C. 標準 CNN の状態方程式(1)x&ij = − xij + z ij + ∑ a kl y kl + ∑ bkl u kl
k ,l∈sijk .l∈S i , j

D. 標準 CNN の出力方程式(2)

yij = f (xi, j )ただし、 f (xi, j ) は非線形関数である。
E. CNN セルの局所的結合式(1)中の A( y ij ) =∑ a y は、近傍セルの出k ,l∈sijkl∑ b u は、

k ,l∈sijklkl

近傍セルの入力(外力)同士の線形結合である。
F. CNN セルの境界条件
CNN セルにおける境界条件の考え方は、3 つあり、外側を 0 と考える固定境界条件、境界を境にして対称な値をとるノイマン境界条件、反対側境界にある変数同士を同じ値と考える周期的境界条件がある。このうち今回のシミュレーションで採用したのは、周期的境界条件である。
3.2 反応拡散系 CNN
パターン生成や波動の生成は、反応拡散方程式として一般に知られている非線形偏微分方程式∂u = f (u ) + D∇ 2 u∂t(3)によって記述される。ただし、
u(x,t)：時間 t における位置ｘのシステム状態を表す m 次元ベクトルf(u)：状態変数ベクトル u の非線形関数
D：m×m の対角行列

∂ 2ul
∂ 2ul
∇ ul =+∂x 2
∂y 2

図６kl

力同士の線形結合であり、 B(u ij ) =2u ≤ b1
⎧m1u + b1 (m0 − m1 )⎪f (u ) = ⎨m0 ub1 ≤ u ≤ b2
⎪m u + b ( m − m )b2 ≤ u202⎩ 2(4 )とする。
複雑な時空間的ふるまいをモデル化するために使われる反応拡散系微分方程式の多くは式(3)で近似できる。
これらの現象は連続体理論を空間的に離散化することにより近似される。その場合、式(3)は、式(6)の非線形性の区分線形近似

表 1 種々の反応拡散系に対するパラメータ値[1]Travelling Trigger Spiral Turing’s Turing’s
wave
wave
wave Patterns Patterns(blobs) (stripes)
m0-1-0.8-0.5-0.1-0.1
m111111m211111b1-1-1-1-1-1b2115101εβ1/35
0.7

1/350.71/35
0.71/3501/100γG

1.231.26211111Du0.3

0.2

0.1

0.1

0.1

Dv0005050

次式で表される。
du (i , j ) = f ( u ( i , j )) + D [ u ( i + 1 , j )dt+ u ( i , j + 1) + u ( i − 1, j )+ u ( i , j − 1 ) − 4 u ( i , j )],(5 )ただし、u(i , j)は、位置(i , j)にあるセルにおける状態変数ベクトルである。
状態変数 u, v を有する 2 次反応拡散 CNN[1]du= − f ( u ) − gv + D u ∇ 2 udt(6 )(dv= ε u + β − γv + D v∇ 2vdt)

図７

孤立セルの null isoclines を描くブロック図%function y = fcn(u)function y = nlfunc1(u)% This block supports an embeddable subset of
the MATLAB language.
% See the help menu for details.
m0=-0.9;
m1=1;
m2=1;
b1=25e-9;b2=45e-9;
if u <= b1
y= m1*u+b1*(m0-m1);
elseif (u > b1 && u<=b2)
y = m0*u;
elsey = m2*u+b2*(m0-m2);end;

図 8 孤立セルの null isoclines

ddt(7)
τdvn = un + βn −γ vn + voff + Dv (vn+1 + vn−1 − 2vn)
ε dt

τ un = − f (un) − gvn +uoff + Du (un+1 +un−1 − 2un)

なお、上式では、一次元の反応拡散 CNN を考えている。
3.4 シミュレーション式(7)を一次元空間上で、32 個の点に離散化しシミュレーションを行う。
A 孤立セルのシミュレーション
まず、はじめに孤立セル (Du = Dv = 0 ) に対して
図 7 のブロック図を用いて、微分方程式(7)の nullisoclines(零等傾線)を求める。null isoclines は、式(7)のそれぞれの式の右辺を零にする (u n , v n ) 曲

図 8.ブロック nlfunc1 の内部構造を考える。非線形項の f(u)は図６に表される区分線形の非線形関数である。
適切なパラメータ設定をすることで、(6)式は、Turing のパターンを形成したり、trigger wave やtraveling wave の伝播をおこなったり、Spiral waveを形成したりする。参考のために、それぞれのパラメータ値を表 1 に示す[1]。
3.3 変数の変換
電子回路の変数として都合が良いように、変数変換を行うと、式(6)は式(7)に書き直される。

線である。
図 ８ は 幾 つ か の パ ラ メ ー タ に 対 す る null
isoclines を示している。
図 8(a)では、興奮しやすい動作を実現するように、セルパラメータが設定されている。この場合、セルは 1 個の安定平衡点 Q1 をもっている。状態が平衡点からずれると、セルは活性化し、弛張期間をへて安定な平衡点 Q1 に戻ってくる。図８(a)では 、セルパラ メータは m1=1, m2=1, m0=-0.9,b1=25×10-9, b2=45×10-9, βn=4.5×10-9, γ=1.2, g=1.1,voff=0, uoff=0, Du=0, Dv=0 である。
図 8(b)では 2 安定動作を実現するように設定されている。

図 9. null isoclines を描くブロック図
図 11. 図 8(a)の場合の孤立セルの解軌道

図 10. 孤立セルの解軌道を描くブロック図
つまり、各孤立セルは 2 安定平衡点 Q1 ,Q3 を持っている。なお、 Q2 は、不安定な平衡点である。
平衡点からずれているとき、一つの平衡点から他の平衡点までトリガーをかけることができる。
図 8(b) で は 、 セ ル パ ラ メ ー タ は m1=1, m2=1,m0=-0.9, b1=2.5×10-9, b2=4.5×10-9, βn=0.7×10-9, γ=5,g=1.2, voff=4×10-9, uoff=-1.35×10-9, Du=0, Dv=0 である。
図８(c)のシミュレーションにおいて、セルのパラメータは、同図のように一つの安定平衡点を持つように設定されるが、特定の結合状態の下では不安定となる。図 8(c)では、セルパラメータは、
m1=1, m2=1, m0=-0.1, b1=1.5×10-9, b2=12.5×10-9,βn=0, γ=6, g=1, voff=0.5 × 10-9, uoff=0.25 × 10-9,Du=0, Dv=0 である。図 11 は、図 8aの場合に対する式(7)の軌道を (u n , v n ) の位相平面上に描いたものである。この計算は、Simulinkでは、図 10 のブロック図を用いて行われる。
図 11 において、図 8(a)のnull isoclinesの図面上に解軌道を描いたことにより、状態変数uとvの初期値をどこにおいても、同じ平衡点Q1に落図 12. 図 8(b)の場合の孤立セルの解軌道

図 13. 図 8(c)の場合の孤立セルの解軌道ち着くことから、図８(a)の null isoclines 曲線の安定な平衡点が一つしかないことが分かる。
この結果は、線形理論と一致する。
次に、図 12 においては、解軌道が微分方程式の初期値を変えることで 2 つの平衡点に落ち着くことがわかり、図 8(b)の null isoclines 曲線は一

図 16. 各セルの状態変数 u の時間波形

図 14. traveling wave を生成する一次元反応拡散CNN(セル数 32)の Simulink メインブロック

図 17. 各セルの状態変数 v の時間波形

図 15. 図 14 で用いたサブシステム(セル)の定義つの安定平衡点 Q1 ,Q3 と一つの不安定の平衡点Q2 を持つことが分かる。
図 13 の場合は、図 8(c)の null isoclines 曲線に１つしか平衡点が存在しないため図 11 の時と同様でどこに微分方程式の初期値をおいてもその平衡点の場所に解軌道が落ち着くことが分かる。
また、興味深いのは図 13 の場合他の二つとは違って平衡点近くで spiral wave のように巻きつきながら平衡点に落ちつくことも分かった(渦型安定平衡点)。

図 18. 250μs ごとに観測した各セルの状態変数 uの遷移図(横軸はセルの番号)：traveling wave
B．一次元反応拡散 CNN のシミュレーション
この節では、反応拡散系 CNN セルを 32 個一次元で結合させた時の traveling wave と trigger waveをシミュレーションする。
このシミュレーションのために用いるブロック
図はそれぞれ 2 つずつある。
まず、図 14 と図 19 は、traveling wave や trigger図 21. 各セルの状態変数 u の時間波形

図 19. trigger wave を生成する一次元反応拡散CNN(セル数 32)の Simulink メインブロック図 22. 各セルの状態変数 v の時間波形

図 20. 図 19 で用いる subsystem(セル)の定義
wave を生成する一次元反応拡散 CNN(セル数 32)の Simulink メインブロックである。
次に、図 15 と図 20 は、メインブロック図で用いたサブシステム(セル)の定義である。
まず、32 個のセルの一次元配列で伝播する
traveling wave の simulink によるシミュレーシ
ョン結果を図 16～図 18 に示す。セルパラメータは 、 τ=42μs, ε=1/35, m1=1, m2=1, m0=-0.9,
b1=25×10-9, b2=45×10-9, βn=4.5×10-9, γ=1.2, g=1.1,voff=0, uoff=0 とし、結合パラメータは、Du=0.5,Dv=0 を用いた。
図 16 と図 17 は、32 個のセルの状態変数 u と v

図 23. 250μs ごとに観測した状態変数 u の遷移図(横軸はセルの番号)：trigger waveの時間波形を示す。図 18 は約 250μ 秒ごとに観測した状態変数 u の遷移図である。縦軸が u の値で横軸はセルの通し番号である。
図 16,図 17 より、図 8(a)の状態で traveling waveが生成しているので、32 個のどのセルも、時間がたつにつれて、状態変数 u と v の値が同じ平衡状態に落ち着いていることが分かる。
また、図 18 により、状態変数 u の値が 250μ秒ごとに移り変わっていき最後に 20nA あたりで平衡点に達している。
次に、32 個のセルの一元配列で伝播させた
trigger waveのSimulinkによるシミュレーション結果を図 21～図 23 で示す。ここでもsimulinkを用いて(7)式を空間的に離散化した式によりシミュレーションを行った。セルパラメータは、
τ=42μs, ε=1/35, m1=1, m2=1, m0=-0.9, b1=2.5×10-9,b2=4.5×10-9, βn=4.7×10-9, γ=5, g=1.2, voff=0,uoff=-1.35×10-9とし、結合パラメータは、Du=0.5,Dv=0 とした。図 21 と図 22 は、32 個のセルの状
態変数uとvの時間波形を示す。
さらに図 23 は約 250μ秒ごとに観測した状態変数 u の遷移図を示している。上図より trigger waveを伝播させた場合でも状態変数 u と v の値は、各セルにおいて初期値を変えても必ず時間が推移するにつれて同じ平衡点の値へと落ち着いていることが分かる。
また、図 23 により、状態変数 u の値が 250μ秒ごとに移り変わっていき最後に 20nA あたりで平衡点に達している。
こ の 結 果 は 、 論 文 [1] で 述 べ ら れ て い る
MATLAB や H-SPICE によるシミュレーション結果と一致していることから、Simulink で得られた上記のシミュレーション結果は妥当なものと考える。
4.まとめ
今回は、Simulinkにより反応拡散系CNN[1]をモデル化し、セルを 32 個一次元配列上に結合させた時のtraveling waveとtrigger waveを再現した。まず、孤立セルの解軌道は、安定平衡点が一つの場合はどこに初期値をおいてもその平衡点に落ち着くことを確かめた。二つの安定平衡点がある状態では、初期値を変えることでQ1とQ3のどちらかに落ち着くことも確認した。次に 32 個結合させた 一 次 元 反 応 拡 散 系 CNN は trigger wave やtraveling waveのどちらの場合を生成させた時でも初期値から安定点となる所まで急激な変化とともに落ち着いていることが分かった。このようにsimulinkを用いることでMATLABのプログラムの知識が無くても微分方程式の算出が理解できブロック図を作成することができシミュレーシ
ョン結果を得ることができるので今後ますますsimulinkを利用されていくのではないかと考えられる。
今後の課題としては、二次元反応拡散 CNN においても Simulink を用いてモデル化しシミュレーションを行いたい。
参考文献[1] T.S-Gotarredona and A.R-Vazquez, On theDesign of Second Order Dynamics
Reaction-Diffusion CNNs, Journal of VLSI
Signal Processing 23, 1999, 351-371

青山貴伸, 森口肇, 蔵本一峰, 吉田郷弘
使える！MATLAB, 講談社, 106-123[3] 高谷邦夫, Matlab の総合応用, 森北出版,117-132[2]