第7章 神経回路シミュレータ

1. NCS の概要

近年,計算機技術は,計算速度の向上や周辺機器の充実など飛躍的な発展を続けています.これに伴い,脳・神経系における情報処理機構に関する新しい研究アプローチとして,生理実験で得られた結果を元に細胞の特性を数理的に記述,再構成し,シミュレーション解析により,その機能メカニズムを探る,再構成的アプローチが注目されています.こうした数理モデルは,通常,多元の非線形連立微分方程式として記述され,数値計算によって解が求められます.しかしながら,このようなシミュレーションプログラムを,汎用プログラミング言語を用いて記述する場合,モデルパラメータの変更に伴う再プログラミングや,ファイル入出力など,モデル記述とは別な非本質的な作業が必要になり研究の効率を著しく低下させることが指摘されています.

そこで,こうした研究を支援するソフトウェアシステムとして神経回路シミュレータ NCS (Neural Circuit Simulator) が開発されました.NCS では細胞の特性や結合状態を専用のモデル記述言語である NCS 言語により記述します.NCS 言語のモデル記述には数値積分,ファイル入出力,計算結果の保存などの非本質的な処理は含まれておらず,研究者はモデル構築作業に専念することができます.また,モデルパラメータ,外部入力刺激などのシミュレーション条件は,シミュレーションプログラムとは別のファイルに格納されており,モデルを書き換えることなく,様々な条件下でのシミュレーションを効率的に行うことができます.このように NCS は神経回路モデルのシミュレーションにおいて,多くのプログラミングの知識を必要とせず,条件を変更しながらシミュレーションを効率良く行えるように設計されたシステムです.また,モデル記述の仕様等は汎用性を考え設計されていますので,神経回路モデルに限らず,連続系モデルであれば NCS を用いてモデリングシミュレーションが可能です.

以上のような NCS の特徴をまとめると,次のようになります.

  • 生理学的知見に基づいた大規模神経回路モデルシミュレーションが可能.

  • 生理学的知見に忠実なモデルはもちろん,連続系モデル一般を全て同様に扱うことができる.

  • 数理モデル構築以外のプログラミングが不要.

  • 再コンパイルをすることなく,条件を変更しながらのシミュレーションが可能.

1.1. 基本仕様

図 7.1. 「NCS のシステム構成」 に NCS システム構成を示します.NCS は 3 つの構成要素,すなわち NCS プリプロセッサ,NCS ライブラリ,条件設定関数群から構成されています.NCS 言語で記述された神経回路モデルは,NCS プリプロセッサを通して C 言語のシミュレーションプログラムとシミュレーション条件ファイル群に変換されます.変換された C 言語のシミュレーションプログラムは,C コンパイラによりコンパイルされたのち,NCS ライブラリとリンクされ,シミュレーション実行ファイルになります.NCS ライブラリは,シミュレーション開始前の処理,数値積分ルーチンなどを含む,シミュレーションを実現するための基本プログラムの集合体です.シミュレーション条件ファイル群は,外部入力条件,モデルパラメータ,信号遅延情報などシミュレーションの実行に関する情報を持ったファイルの集まりであり,条件設定関数群により,各種の条件に設定することができます.これらの条件ファイル群と実行ファイルによりシミュレーションが行われます.

NCS のシステム構成

図 7.1. NCS のシステム構成

1.2. NCS と SATELLITE

NCS は当初 MS-DOS 上で開発され利用されてきましたが,UNIX ワークステーション上で構築された SATELLITE のモジュールの一つとして機能するように変更が加えられると共に機能の追加,改良がなされてきました.すなわち,NCS と SATELLITE の統合により,モデリングからシミュレーション結果の解析,表示,評価に至る一連の作業を円滑に進める環境が実現されました.現在 NCS は,SATELLITE の理念をさらに推し進めた SATELLITE 言語上で稼動しており,一層の利便が図られています.

2. NCS 言語

NCS では,モデルを記述するために,NCS 言語と呼ばれる専用の言語を用いています.NCS にはより効率よくモデル記述を行うため,いくつかの特殊記述が用意されています.

NCS 言語は,予約語,NCS ライブラリ関数,文および特殊記述で構成されています.例えば,図 7.2. 「RL 直列回路」 に示す RL 直列回路は 図 7.2. 「RL 直列回路」 に示す微分方程式で記述されますが,この回路の NCS による記述は,リスト 1 のようになります.

NCS では,リスト 1 に示すように "module" 文から始まり "end" で終わる,あるシステムに関する一連のモデル記述を 1 つの単位とします.NCS は通常シングルモジュールからなるモデルのみ扱いますが,その拡張として複数のモジュールからなるモデルを扱うことも可能です.このとき,モジュール記述の他に,モジュール同士の結合状態を記述するネットワーク記述が必要になります.ネットワーク記述は,NCS が神経回路を対象にしていることから,神経回路を構成しやすいよう特化した記述になっています.ネットワーク記述をしようした複数のモジュールからなるモデルを NCS 言語で記述する方法などは 項4.2. 「NCS 言語による結合モデル」 で詳しく説明します.

以下には NCS 言語で使用される,予約語,NCS ライブラリ関数,文,特殊記述の詳細について述べていきます.

RL 直列回路

図 7.2. RL 直列回路

表 7.1. RL 直列回路モデル

電流

: 電流 [A]

: 電圧源 [V]

: インダクタンス 0.1[H]

インダクタンス電圧

: インダクタンス電圧 [V]

抵抗電圧

: インダクタンス電圧 [V]

: 抵抗 10[Ω]

リスト 1 : NCS言語による RL 直列回路モデルの記述

  1   /* RL circuit module */
  2   module:      circuit;
  3   exinput:     E;
  4   output:      i;
  5   observable:  Vr, Vl;
  6   constant:    L = 0.1;
  7   parameter:   R = 10., i0 = 0.;
  8   function:
  9                di = (E - Vr) / L;
  10               i = integral(i0,di);     /*  電流  */
  11               Vl = di * L;             /*  インダクタンス電圧  */
  12               Vr = R * i;              /*  抵抗電圧  */
  13               end;

2.1. 予約語

次の 2 つが予約語としてシステムにより定義されています.

  1. TIME : シミュレーション時刻を保持しています.全モジュールにおいて適用できます.

  2. CN : コンポーネント番号を保持しています.全モジュールにおいて適用できます.

2.2. ライブラリ関数

次の 5 種類の関数が用意されています.

  • pulse

    記述)

    y = pulse(a,b,c,d,e);

    引数)

    1. a : 入力開始時刻

    2. b : 入力初期値

    3. c : 高さ

    4. d : 時間幅

    5. e : 周期

    パルス関数

    図 7.3. パルス関数

  • ramp

    記述)

    y = ramp(a,b,c);

    引数)

    1. a : 入力開始時刻

    2. b : 入力初期値

    3. c : 勾配

    ランプ関数

    図 7.4. ランプ関数

  • integral

    記述)

    y = integral(a,p);

    式)

  • sigmoid

    記述)

    y = sigmoid(-x);

    式)

  • rcasb

    記述)

    y = rcasb(v,a,b,c,d,e,f,g);

    式)

NCS ライブラリ関数の他に,記述には C 言語の数学ライブラリ関数が使用できます.表 7.2. 「NCS 言語で使える C 言語数学ライブラリ関数の例」 に,その代表的な例の記述と式を示します.

また,ユーザが定義した関数も使用できます.

表 7.2. NCS 言語で使える C 言語数学ライブラリ関数の例

関数名 記述
exp
pow
sin
cos
tan

2.3. モジュール記述

あるシステムに関する一連のモデル記述をモジュールとして記述します.電気回路モデルの記述や膜モデルによるイオン電流レベルでの記述が可能です.モジュールはいくつかの文を用いて記述します.ここでは,モジュール記述に使用する文とその構成について説明します.

各文の基本的な構成は次のようになります.

宣言子  :  文の内容  ;

どの文の記述を行うか宣言する宣言子で始まり,コロン ":" で区切って文の内容を記述します.各文はセミコロン ";" で終わらなければいけません.文の中には以下の,"module","output","function" は必ず記述しなければいけません.その他の文は状況に応じて使用します.

  1. module [必須]

    モジュール名を定義します.先頭が英文字で始まる 8 文字以内の英数字綴りをモジュール名として宣言します.

    記述例)

    /*  "circuit"というモジュール名を宣言  */
    module  :  circuit;
  2. exinput [任意]

    セルモジュールにおいて,モデルの外部入力変数名を定義します.外部変数名は 1 つです.

    記述例)

    /*  "E"を外部入力変数として定義します.  */
    exinput  :  E;
  3. output [必須]

    モジュールの出力変数名を定義します.出力変数の数は 1 つです.

    記述例)

    /*  "i"を出力変数として定義  */
    output  :  i;
  4. observable [任意]

    後述する function 文で使用する変数のうち,計算値を観測したい変数名を宣言します.宣言した変数はシミュレーション実行時に,nout 関数で出力指定を行うことによって値を観測できます.

    記述例)

    /*  "Vr","Vl"を観測する変数として定義  */
    observable  :  Vr, Vl;
  5. constant [任意]

    function 文で使用する定数の,定数名とその値を宣言します.ここで定義された定数名は必ずしも function 文で使用される必要はありません.ただし,その分メモリを取っていますので,不必要な定数は宣言しないほうが良いでしょう.

    記述例)

    /*  "L"を定数として使用し,0.1を代入することを宣言  */
    constant  :  L = 0.1;
  6. parameter [任意]

    パラメーター変数の定義を行います.ここで定義されたパラメーターはシミュレーション条件ファイルに登録され,シミュレーション実行時に npara 関数によって値を変更することができます.ここで定義されたパラメーターは必ずしも function 文で使用される必要はありません.

    記述例)

    /*  "R","i0"をパラメーターとして使用し,
        それぞれ10,0を代入することを宣言  */
    parameter  :  R = 10, i0 = 0;
  7. function [必須]

    数式によりモジュールの特性を記述します.変数は使用される前に必ず値が定義されていなければなりません.数式はセミコロンで区切られます.

    function 文内の記述に柔軟性を持たせるため,if 文が用意されています.

    書式)

      if ( 条件式 ) { 文1 } [ else { 文2 } ]

    表 7.3. 条件式の関係演算子

    演算 記号
    等値 =
    非等値 <>
    比較 <,>,<=,>=
    論理和 .or.
    論理積 .and.
    否定 .not.

    条件式が真なら文 1 が,偽なら文 2 が実行されます.条件式に用いる関係演算子を 表 7.3. 「条件式の関係演算子」 に示します.

    記述例)

    /*  "Vl","Vr"の式を記述  */
    function  :
    Vl = di * L;
    Vr = R * i;
  8. コメント

    "/*","*/" で囲まれた部分はコメントになります.文中のどこでも,また複数行に渡る使用も可能です.

文の並びは module 文を最初,function 文を最後に記述する以外は順不動で,module,output,function 文以外の文は省略することが可能です.なお,モジュール記述の終わりに必ず,"end;" をつけることを忘れないでください.

2.4. モデル記述例

ここでは,NCS 言語で記述する方法を,実際のモデルを使って説明します.

●RL 直列回路モデル

ここでは,前節の 図 7.2. 「RL 直列回路」 の RL 直列回路のモデル記述について説明します.RL 直列回路モデルのインダクタンス電流は以下のように表されます.

式 (7.1),(7.2) を NCS 言語によって記述すると,

Vr = R * i;
Vl = di * L;

回路を流れる電流 は微分方程式 (7.3) で記述されていますので,次のようになります.

di = (E - Vr) / L;
i = integral(i0,di);

次に,こうした特性を持つセルタイプモジュールを記述します.モジュール名は "circuit" とします.

module  :   circuit;

外部入力に電圧源 を,そしてインダクタンス電圧 を出力するものとします.

exinput:   E;
output:    Vl;

またここで抵抗電圧 ,電流 を観測することにします.

observable:   Vr, i;

インダクタンス を定数として宣言します.

constant:   L = 0.1;

抵抗 ,初期電流 をパラメーターとして宣言します.

parameter:   R = 10., i0 = 0;

function 文には,前述した記述を用います.

こうして,リスト 1 の NCS 言語による記述が完成します.

●振り子の運動モデル

ここでは,力学系の簡単なシミュレーションとして,図 7.5. 「振り子の運動」 に示す振り子の運動について説明します.一端を固定した長さ [m] の糸の他端に重さ [kg] のおもりをつるし,鉛直面内で振動させます.

振り子の運動

図 7.5. 振り子の運動

最下点から円弧にそって測った長さを [m]( ) とすれば,その運動方程式は

と表せます.これより NCS のモデル記述は

module:      PENDULUM;
output:      theta;
observable:  d2theta,dtheta,theta;
parameter:   m = 1.0,g = 9.8,l = 1.0,theta0 = 3.141592 / 180.0 * 45.0;
function:
             d2theta = -g / l * sin(theta);
             dtheta = integral(0.0,d2theta);
             theta = integral(theta0,dtheta);
end;

となります.運動方程式は 2 階微分の形になっていますので,これを 1 階の連立微分方程式に書き換え, を求める記述とします.

比較のため解析的に運動方程式の解を求めます.問題を簡単にするために [rad] を 1 より充分小さいと仮定し,平衡点のまわりで線形近似します.このとき運動方程式は

となります.これを解くと

となります.

モデルファイル名を "pendulum.mdl" とし,シミュレーションするための SATELLITE スクリプトをリスト 2 に示します.シミュレーションの実行ついては,項3. 「NCS の使用法」 で詳しく説明します.

の初期値を 10 度として振り子を運動させた (減衰なし) 場合のシミュレーション結果を 図 7.6. 「振り子の運動モデルのシミュレーション結果」 に示します.実践が解析解,四角形は 10 刻み毎のモデル出力になります.

リスト 2 : 振り子の運動モデルのシミュレーション実行用バッチスクリプトの例

1.  npp("pendulum")     #モデルファイルの登録とプリプロセッサの起動
2.  nlink();                                  # 実行ファイルの作成
3.  ntime(10.0,0.01,0.01,0.01);               # 時間条件の設定
4.  time = (0~(10.0 / 0.01)) * 0.01;          # 時間軸データ
5.  series theta;
6.  nout(theta,"pendulum",0,1);               # 外部出力条件の設定
7.  theta0 = PI / 180 * 10.0;                 # 初期角度10度に設定
8.  npara("PENDULUM","theta0",theta0);        # パラメータの変更
9.  ncal();                                   # シミュレーションの実行
10. phi = PI / 180 * 90;                      # 解析解の計算
11. m = 1.0;
12. l = 1.0;
13. gra = 9.8;
14. true_theta = theta0 * sin(sqrt(gra / l) * time + phi);
15. theta = theta / (PI / 180.0);             # [rad]から[deg]に変換
16. true_theta = true_theta / (PI / 180.0);
17. wopen(1,"A4",0,1);                        # グラフ表示
18. title(1,"time[s]","angle[deg]");
19. scale("N","F","N","F",0,10,-30,30);
20. axis(1,1,"XY","XY",5,0,3,2,15,0);
21. frame();
22. graph(theta,time,0,1,10,0,2);
23. graph(true_theta,time,0,0,0,0,0);
振り子の運動モデルのシミュレーション結果

図 7.6. 振り子の運動モデルのシミュレーション結果

●Hodgkin-Huxley モデル

ここではイカの巨大神経軸索の電気的特性をモデル化した Hodkin-Huxley (以下 H-H モデル) の記述について詳しく説明します.

H-H モデルのナトリウム電流は次の式で表されます.

: ナトリウムの膜コンダクタンス,120[ ]

: ナトリウムの反転電位,115[mV]

式 (7.9),(7.10) を NCS 言語によって記述すると,

am = 0.1 * (25. - V) / (exp((25. - V) / 10.) - 1);
bm = 4. * exp(-V / 18.);

ただし, は V = 25[mV] 時の値を別に与える必要がありますので,次のように書き換えます.

if(V != 25.){
am = 0.1 * (25. - V) / (exp((25. - V) / 10.) - 1);
}else{
am = 0.1 * 10.;
}

ナトリウムチャンネルの開閉確率 は微分方程式で記述されていますので,次のようになります.

dmNa = am * (1. - mNa) - bm * mNa;
mNa = integral(mNa0,dmNa);

同様にして の記述は式 (7.11),(7.12),(7.13) より次のようになります.

ah = 0.07 * exp(-V / 20.);
bh = 1. / (exp((30. - V) / 10.) + 1.);
dhNa = ah * (1. - hNa) - bh * hNa;
hNa = integral(hNa0,dhNa);

最終的にナトリウム電流 は,式 (7.7) より次のようになります.

INa = GNa * pow(mNa,3.0) * hNa * (V - VNa);

カリウム電流は次式で表されます.

: カリウムの膜コンダクタンス,36[ ]

: カリウムの反転電位,-121[mV]

カリウムチャンネルの開閉確率はナトリウムの場合と同様にして,式 (7.15),(7.16),(7.17) より,次のように記述されます.

if(V != 10.){
     an = 0.01 * (10. - V) / (exp((10. - V) / 10.) - 1.);
}else{
     an = 0.01 * 10;
}
bn = 0.125 * exp(-V / 80.);
dnK = an * (1. - nK) - bn * nK;
nK = integral(nK0,dnK);

カリウム電流は式 (7.14) より,次のように記述されます.

IK = GK * pow(nK,4.0) * (V - VK);

漏れ電流は次式で表されます.

: 漏れ電流の膜コンダクタンス,0.3[ ]

: 漏れ電流の反転電位,10.6[mV]

上式より,漏れ電流の記述は次のようになります.

Il = Gl * (V - Vl);

膜を横切る電流 は,次式で表されます.

注入電流を とすると,全電流は次式で表されます.

これより,式 (7.19) は,次式のように変形できます.

以上より,NCS 言語での記述は次のようになります.

Iall = Iex - INa - IK - Il;
dV = Iall / Cm;
V = integral(V0,dV);

次に,こうした特性を持つセルタイプモジュールを記述します.モジュール名は,"hh" とします.

module:   HH;

外部入力に注入電流 ,膜電位 を出力するものとします.

exinput:   Iex;
output:    V;

ここで,ナトリウム電流 ,カリウム電流 ,漏れ電流 を観測することにします.

observable:   INa, IK, Il;

各イオン電流の反転電位 を,定数として宣言します.

constant:   VNa = 115.0, VK = -12.0, Vl = 10.6;

膜容量 ,積分の初期値,各イオンの膜コンダクタンス をパラメーターとして宣言します.

parameter:   Cm = 1.0, mNa0 = 0.05293, GNa = 120., GK = 36.,
             Gl = 0.3, hNa0 = 0.5961, nK0 = 0.3177, V0 = 0.;

function 文には,前述した記述を用います.こうして,リスト 3 の NCS 言語による記述が完成します.

パルス状の電流を注入した時のモデルシミュレーションを実行し,結果をグラフに表示するスクリプトの例をリスト 4 に示します.実行方法の詳細は 項3. 「NCS の使用法」 で説明します.結果を 図 7.7. 「H-H モデルのシミュレーション結果」 に示します.

リスト 3 : NCS 言語による Hodgkin-Huxley モデルの記述

1   /*   HH module   */
2   module:      hhmodel;
3   exinput:     Iex;
4   output:      V;
5   observable:  INa, IK, Il, Ig;
6   constant:    VNa = 115.0, VK = -12.0, Vl = 10.6;
7   parameter:   Cm = 1.0, mNa0 = 0.05293, GNa = 120.0, GK = 36.0,
                 Gl = 0.3, hNa0 = 0.5961, nK0 = 0.3177, V0 = 0.0;
8   function:
9                if(V != 25.0){
10                    am = 0.1*(25.0-V)/(exp((25.0-V)/10.0)-1.0);}
11               else{
12                    am = 0.1*10.0;}
13               bm = 4.0*exp(-V/18.0);
14               dmNa = am*(1.0-mNa)-bm*mNa;
15               mNa = integral(mNa0,dmNa);
16               ah = 0.07*exp(-V/20.0);
17               bh = 1.0/(exp((30.0-V)/10.0)+1.0);
18               dhNa = ah*(1.0-hNa)-bh*hNa;
19               hNa = integral(hNa0,dhNa);
20               INa = GNa*pow(mNa,3.0)*hNa*(V-VNa);
21               if(V != 10.0){
22                    an = 0.01*(10.0-V)/(exp((10.0-V)/10.0)-1.0);}
23               else{
24                    an = 0.01*10.0;}
25               bn = 0.125*exp(-V/80.0);
26               dnK = an*(1.0-nK)-bn*nK;
27               nK = integral(nK0,dnK);
28               IK = GK*pow(nK,4.0)*(V-VK);
29               Il = Gl*(V-Vl);
30               Iall = Iex-INa-IK-Il;
31               dV = Iall/Cm;
32               V = integral(V0,dV);
33               end;
H-H モデルのシミュレーション結果

図 7.7. H-H モデルのシミュレーション結果

リスト 4 : H-H モデルのシミュレーション実行用バッチスクリプトの例

1   nassign("hhmodel");                       # モデルファイルの登録
2   npp();                                    # プリプロセッサの起動
3   nlink();                                  # 実行ファイルの作成

4   ntime(10,0.001,0.01,1);                   # 時間条件の設定
5   nstim("hhmodel",0,"P",1,0,100,3,999);     # 外部入力条件の設定

6   series V,Iin,INa0;
7   nout(Iin,"hhmodel",0,2);                  # 外部出力条件の設定
8   nout(V,"hhmodel",0,1);
9   nout(INa0,"hhmodel",0,3,"INa");
10  ninteg("R");                              # 積分方式の変更
11  ncal();                                   # シミュレーションの実行

12  wopen(1,"A4",0,0);                        # グラフ表示
13  time=(0~1000)/100;
14  graph(V,time,0,0,0,1,3.5);
15  title(1,"time[ms]","membrane potential[mv]");
16  axis(1,1,"XY","XY",5,0,3,0,0,0);

3. NCS の使用法

本システムは SATELLITE 上で稼動しており,会話型処理,またはバッチファイルを用いた一括処理ができます.コマンドについての詳細は,SATELLITE 言語リファレンスマニュアルを参照して下さい.

NCS でシミュレーションを行う場合の手順は,次のようになっています.

  1. モデルファイルの作成

    NCS 言語を用いてモデルを記述したファイルを作成します.NCS 言語を用いた記述方法の詳細は 項2. 「NCS 言語」 を参照してください.

  2. モデルファイルの登録

    シミュレーションを行うモデルの,モデル記述ファイルを登録します.NCS は以後,このモデルファイルに対して処理を行います.

  3. 実行ファイル,シミュレーション条件ファイルの作成

    登録されたモデルファイルに対してプリプロセッサを起動,そしてリンクを行うことにより,実行ファイルとシミュレーション条件ファイル群を作成します.

  4. シミュレーション条件の設定

    シミュレーションの際に必要なシミュレーション時間,外部入力,出力変数などの設定を行います.

  5. シミュレーションの実行

    シミュレーションを実行します.

  6. シミュレーション結果の表示・解析

    シミュレーションが終わって得られた結果のグラフを表示します.

以下にはこの手順の詳細を,リスト1のモデル記述を使ったシミュレーションを例に用いて述べていきます.

3.1. モデルファイルの作成

NCS 言語の仕様に従って,モデルファイルを記述します.モデルファイルのファイル名は拡張子に ".mdl" を付けてください.

注) ver2.95

次項で説明する nassign 関数でモデルファイル名が登録してあるとき,ne 関数を実行すると,そのモデルファイルを対象としてエディタが実行されます.この時立ち上がるエディタのデフォルトは vi エディタですが,UNIX の環境変数 editor を設定することによって,好きなエディタに変更できます.

3.2. モデルファイルの登録

後述する npp 関数によるプリプロセッサの起動や,nlink 関数による実行ファイルの作成は,SATELLITE のワークエリア内に登録されたモデルファイルに対して実行されます.従って,NCS を使ってシミュレーションを行うためには,使用するモデルを登録する必要があります.このモデルファイルの登録は nassign 関数か,次項に述べる npp 関数によってなされます.npp 関数は引数にモデルファイル名をとり,そのファイル名のモデルファイルを登録します.このモデルファイル名には拡張子 ".mdl" を付ける必要はありません.

[]satellite[]~/home/demo:[1]% nassign("circuit")

上記の例では,"circuit.mdl" をモデルファイルとして登録しています.

3.3. 実行ファイル,シミュレーション条件ファイルの作成

npp 関数実行により,プリプロッセサが起動し,NCS 言語で記述されたモデルファイルが C 言語ソースファイルとシミュレーション条件ファイル群に変換されます.

[]satellite[]~/home/demo:[2]% npp()

とすることにより,プリプロセッサが起動します.

npp 関数はモデルファイル名を引数にとることができます.nassign 関数でモデルファイルを登録しなくても,npp 関数でモデルファイル名を指定することで登録が可能です.この場合も nassign 関数同様,モデルファイル名に拡張子 ".mdl" を付ける必要はありません.

[]satellite[]~/home/demo:[3]% npp("circuit")

とすることにより,"circuit.mdl" をモデルファイルとして登録して,プリプロセッサを起動することになります.

npp 関数実行後,nlink 関数により,C 言語のソースファイルをコンパイル後,NCS のライブラリとリンクされ,実行ファイルが作成されます.

[]satellite[]~/home/demo:[4]% nlink()

3.4. シミュレーション条件の設定

シミュレーションを行うにあたって,各種シミュレーション条件を設定します.このシミュレーション条件は npp 関数が実行された時点でモデルファイルより設定されているものと,設定されていないものが存在します.シミュレーションを行うにあったて,設定の必要なシミュレーション条件がありますので注意してください.

以下に述べるシミュレーション条件の設定関数は,npp 関数実行によって生成されたシミュレーション条件ファイル群に対して行われます.従って,シミュレーション条件設定関数実行以前に,npp 関数が実行されている必要があります.npp 関数が実行されていない状態で,シミュレーション条件設定関数を実行した場合,

sl : Error [<NCS:nout> No.1] : Improper Model File Name
near at line <n>

というエラーが表示されます.また,各種シミュレーション条件を設定後,再度 npp 関数を実行すると,シミュレーション条件は初期化されます.

●時間条件

時間条件は ntime 関数によって設定します.ntime 関数の書式は次のようになっています.

書式) ntime( last, cal, str, itv )

引数)

  1. last : 計算時間

  2. cal : 計算刻み

  3. str : 計算結果の抽出時間幅

  4. itv : バッファに書き込む時間幅

計算時間はシミュレーションが終わるまでの時間です.計算刻みは,数値計算の際の時間刻み幅で,これの大小によって誤差や,収束度,シミュレーションの実行時間が左右されます.適当と思われる値に設定してください.ストア刻みは,計算結果を抽出する時間幅で,この値を小さくすればするほど,細かい時間幅の実行結果が得られます.ただし,データ容量の面から言えば,実行結果の容量が大きくなります.使用しているコンピューター環境を考えて適切な値に設定してください.書き込み刻みは,実行結果を後で説明する nout 関数で指定されたバッファに書き込む時間幅です.書き込み刻みを小さくしすぎると,シミュレーション実行時間が非常に長くなりますので注意してください.

[]satellite[]~/home/demo:[5]% ntime(0.1,0.005,0.005,0.05)

としたときには,計算時間は 0.1[s],計算刻みは 0.005[s],計算結果の抽出時間幅を 0.005[s],バッファに書き込む時間幅は 0.05[s] と設定したことになります.

設定した時間条件は nsclist 関数の引数に "T" を用いて表示することができます.

[]satellite[]~/home/demo:[6]% nsclist("T")

TIMER

  Last  Time =    0.1
  Calc. Step =  0.005
  Store Step =  0.005
  BufferStep =   0.05

nsclist 関数を用いることによって,時間条件が正しく設定されているか確認してください.

●外部入力条件

外部入力条件は nstim 関数によって設定します.nstim 関数はモジュール名,コンポーネント番号,入力波形を指定する書式になっています.入力波形には,パルス関数,ランプ関数とファイル入力,バッファ入力が用意されています.nstim 関数の基本的な書式は次のようになっています.

書式) nstim( mdul, com, type, p1 [, p2, p3, p4, p5] )

引数)

  1. mdul : モジュール名

  2. com : コンポーネント番号

  3. type : 入力波形の指定 ("P", "R", "B")

  4. : パラメーター列

入力波形
パルス関数 入力開始時刻 入力初期値 高さ 時間幅 周期
ランプ関数 入力開始時刻 入力初期値 傾き    
バッファ入力 バッファ名        
パルス関数

図 7.8. パルス関数

ランプ関数

図 7.9. ランプ関数

この関数で指定するモジュールでは,exinput 文による外部入力の記述がなされている必要があります.

[]satellite[]~/home/demo:[7]% nstim("circuit",0,"P",0,1,0,10,999)

とすることにより,モジュール名 "circuit" のモジュールに対してパルス関数を入力する設定を行います.シングルモジュールからなるモデルを扱う場合は,コンポーネント番号を必ず 0 に設定します.

設定した外部入力条件は nsclist 関数の引数に "S" を用いて表示することができます.

[]satellite[]~/home/demo:[8]% nsclist("S")

EXTERNAL  INPUTS
   Data No.      Component    Function Data No.
   1          circuit(0)

No.1  Function     start_tm init_out height width   period
      <  PULSE>    [   0]   [   1]   [   0] [   10] [   999]

nsclist 関数を用いることによって,外部入力条件が正しく設定されているか確認してください.

●出力条件

出力条件は nout 関数によって設定します.nout 関数は出力値を格納するバッファ名,モジュール名,コンポーネント番号,出力情報の属性を引数に持ちます.ただし,出力情報の属性を内部変数値に設定した場合のみ,内部変数名を引数に指定する必要があります.nout 関数の基本的な書式は次のようになっています.

書式) nout( buff, mdl, com, type [, val] )

引数)

  1. buff : 出力値を格納するバッファ名

  2. mdl : モジュール名

  3. com : コンポーネント番号

  4. type : 出力情報の属性 (1 : 出力値,2 : 入力値,3 : observable に指定した変数値)

  5. val : type = 3 を指定したとき,observable に指定した変数名を指定する.

ここで,出力値とはモジュール名 mdl のモジュールの記述において output 文で指定されている変数の値が入ります.また入力値は,exinput 文で指定されている変数の値が入ります.

なお,出力値を格納するバッファは nout 関数実行以前に series 型のバッファとして宣言されている必要があります.例えば,nout 関数実行以前に,

[]satellite[]~/home/demo:[9]% series Vin, i, vr, vl

とすれば,nout 関数で,バッファ Vin,i,vr,vl を出力値の格納用バッファとして扱うことができます.バッファが宣言されていない場合には,nout 関数実行時にエラーが表示されます.注意してください.

[]satellite[]~/home/demo:[10]% nout(Vin,"circuit",0,2)
[]satellite[]~/home/demo:[11]% nout(vl,"circuit",0,1)

によって,モジュール名 "circuit" のモジュールの入力値がバッファ Vin に,出力値が vl にそれぞれ出力されます.

また,その他の内部変数を出力するためには次のように設定します.

[]satellite[]~/home/demo:[12]% nout(vr,"circuit",0,3,"Vr")
[]satellite[]~/home/demo:[13]% nout(i,"circuit",0,3,"i")

これによって,モジュール名 "circuit" のモジュールの内部変数 "Vr" の値がバッファ vr に,"i" の値がバッファ i にそれぞれ出力されます.

設定した外部出力条件は nsclist 関数の引数に "O" を用いて表示することができます.

[]satellite[]~/home/demo:[14]% nsclist("O")

OUTPUT
   Variable    OUTPUT    VARIABLE
   Vin [   0]  EX.INPUT  OF circuit(0)
   vl  [   0]  OUTPUT    OF circuit(0)
   vr  [   0]  Vr        OF circuit(0)
   i   [   0]  i         OF circuit(0)
Max of Number of Output = 255

●パラメーター値の変更

パラメーター値は npara 関数によって変更することができます.npara 関数の書式は次のようになっています.

書式) npara( mdl, var, num )

引数)

  1. mdl : モジュール名

  2. var : パラメーター名

  3. num : パラメーターに設定する値

npara 関数の実行にあたって,モデルファイルの記述で,モジュール名 mdl の変数 var が parameter 文で宣言されている必要があります.

[]satellite[]~/home/demo:[15]% npara("circuit","R",20)

npara 関数で変更した値は,npp 関数等によってシミュレーション条件ファイル群が初期化されるまで有効になります.

現在のパラメーター値は,nlist 関数で表示することができます.

[]satellite[]~/home/demo:[16]% nlist("circuit")
Parameter    List
Module name :   circuit
   R        =   [    20]
   i0       =   [     0]

nlist 関数はモジュール名を引数とし,そのモジュールのパラメーター値を表示します.

積分方式の選択

通常,数値積分法には自動刻み積分法が用いられています.NCS はこの他に,ルンゲクッタ - ギル法とオイラー法を備えており,ninteg 関数によって,これらの積分法を用いることができます.ninteg の書式は次のようになっています.

書式) ninteg( type )

引数)

  1. type : 積分方式 ("F" : 自動刻み積分法,"R" : ルンゲクッタ‐ギル法,"E" : オイラー法)

3.5. シミュレーションの実行

各種シミュレーション条件を設定後,ncal 関数によりシミュレーションが実行されます.

[]satellite[]~/home/demo:[17]% ncal()

と ncal 関数を実行すると,

## NCS Ver.6.8.3 SIMULATION PROGRAM on osf4.0 ##
>> THE CALCULATION HAS FINISHED ................... !! << 0.0% done.

と表示されます."0.0%" は現在の実行の割合をパーセンテージで示しており,この数字が "100%" に達すると,

## NCS Ver.6.8.3 SIMULATION PROGRAM on osf4.0 ##
>> THE CALCULATION HAS FINISHED ................... !! << 100.0% done.

と表示が変わり,シェルは関数待ち状態になって,シミュレーションが終了したことを示します.

3.6. バッチファイルの使用

ここまで述べてきた,シミュレーションを実行するまでを 1 つのバッチファイルにまとめることができます.そして,inline 関数を使用することにより,バッチファイルを用いたシミュレーションの実行が行えます.

リスト 5 : シミュレーションを行うためのバッチスクリプトの例

1   nassign("circuit");                     # モデルファイルの起動
2   npp();                                  # プリプロセッサの起動
3   nlink();                                # 実行ファイルの作成
4   ntime(0.1,0.005,0.005,0.05);            # 時間条件の設定
5   nstim("circuit",0,"P",0,1,0,10,999);    # 外部入力条件の設定
6   series Vin,i,vr,vl;
7   nout(Vin,"circuit",0,2);                # 外部出力条件の設定
8   nout(i,"circuit",0,1);
9   nout(vr,"circuit",0,3,"vr");
10  nout(vl,"circuit",0,3,"vl");
11  ninteg("R");                            # 積分方式の変更
12  ncal();                                 # シミュレーションの実行
[]satellite[]~/home/demo:[18]% inline("バッチファイル名")

リスト 5 にモデルファイル circuit.mdl を使ったシミュレーションを行うための,バッチファイル例を示します.このファイルを "circuit.sl" としますと,シミュレーションの実行は次のようになります.

[]satellite[]~/home/demo:[19]% inline("circuit.sl")

3.7. シミュレーション結果の表示・解析

計算結果は,nout 関数で指定したバッファに格納され,SATELLITE の他のシステムモジュールによって,結果のグラフ表示,解析などが行えます.例えば,circuit.sl を用いたシミュレーション結果は,GPM を使用することにより,次のようにしてグラフに表示することができます.

[]satellite[]~/home/demo:[20]% wopen(1,"A4",0,0)
[]satellite[]~/home/demo:[21]% time = (0~20) / 200
[]satellite[]~/home/demo:[22]% graph(i,time,0,0,0,0,0)
[]satellite[]~/home/demo:[23]% title(1,"time[s]","current[a]")
[]satellite[]~/home/demo:[24]% axis(1,1,"XY","XY",5,0,0,0,0,0)

グラフ表示の詳細は SATELLITE システムの GPM モジュールマニュアルを参照してください.

RL 直列回路モデルのシミュレーション結果

図 7.10. RL 直列回路モデルのシミュレーション結果

4. NCS による神経回路網のモデリングシミュレーション

4.1. モジュール化の概念

ここでは NCS におけるタイプ,モジュール,コンポーネントの概念を述べます.これらの関係を 図 7.11. 「NCS における神経回路素子とモデル構造」図 7.12. 「実際の神経回路との対応」 に示します.

NCS における神経回路素子とモデル構造

図 7.11. NCS における神経回路素子とモデル構造

実際の神経回路との対応

図 7.12. 実際の神経回路との対応

神経回路は,多数の神経細胞が化学シナプス,電気シナプス (ギャップジャンクション) を介して結合することにより構成されています.すなわち,神経回路は,神経細胞,化学シナプス,電気シナプスの 3 種類の組み合わせにより成り立っています.神経回路の各構成要素はそれぞれ特定の機能をもっており,神経細胞は信号の処理,化学シナプスは化学伝達物質による情報伝達,電気シナプスは電位差による情報の伝達という機能を果たしています.NCS ではモデル記述を行う上で,これら 3 つの要素を基本の型 (タイプ) と考え,神経回路を神経細胞型 (cell type : セルタイプ),化学シナプス型 (synapse type : シナプスタイプ),電気シナプス型 (gap type : ギャップタイプ) としています.つまり,

神経回路の構成要素 = { cell type, synapse type, gap type }

と表記できます.

各タイプはさらに細分化でき,それらをモジュールと呼びます.モジュールは,各タイプの中で同一の特性を持つ素子全体を表現するものです.例えば網膜の場合,セルタイプモジュールには光を受容する錐体視細胞 (CONE),視覚物質を修飾,制御する水平細胞 (HC) など特性の異なる神経細胞があります.そこで,NCS では各タイプを異なる特性を有する素子ごとに分割し,分割された一つ一つをモジュールと定義します.セルタイプは

cell type = { HC module, CONE module, … }

と表され,他のタイプも同様です.

モジュールは同じ特性を持つ素子一つ一つの集合であり,HC モジュールは次のように表現できます.

HC module = { HC[0], HC[1], … }

HC[0],HC[1],…はコンポーネントと呼ばれ,これは神経回路網を構成する素子の実態に対応します.コンポーネントはそれの属するモジュール名に数字を添え,番号付けすることにより表され,一つの素子が特定されます.

4.2. NCS 言語による結合モデル

NCS では大規模なモデルを記述しやすく,かつ柔軟性のあるモデル記述を実現するため,前節のモジュール化の概念のもと,結合記述を用いることができます.例えば,表 7.4. 「Hodgkin-Huxley モデル」 に示した Hodgkin-Huxley モデル (以下 H-H モデル) が 図 7.13. 「神経回路モデルの例」 のように抵抗で直列に継がっている神経回路モデルの記述は,リスト 6 のようになります.

  • 2~12 行目

    ネットワーク記述部.図 7.13. 「神経回路モデルの例」 の結合状態を定義しています.

    神経回路モデルの例

    図 7.13. 神経回路モデルの例

  • 14~47 行目

    H-H モデルの記述,ギャップ性電流 を入力として,膜電位 が出力されます.function 文は,表 7.4. 「Hodgkin-Huxley モデル」 の連立微分方程式の記述です.

    表 7.4. Hodgkin-Huxley モデル

    膜電位

    : 膜電位[ ]

    : 膜容量[ ]

    : 全膜電流[ ]

    ナトリウム電流

    : ナトリウム電流 [ ]

    : ナトリウムの膜コンダクタンス 120[ ]

    : ナトリウムの反転電位 115[ ]

    カリウム電流

    : カリウム電流[ ]

    : カリウムの膜コンダクタンス 36[ ]

    : カリウムの反転電位 -12[ ]

    リーク電流

    : リーク電流の膜コンダクタンス 0.3[ ]

    : リーク電流の反転電位 10.6[ ]

  • 49~56 行目

    電気シナプスの記述.膜電位に応じて電流を出力する記述になっています.

このように NCS 言語による記述は,微分方程式などで表されるモジュール記述 (リスト 6 の 14~62 行目),モジュール記述の結合状態を表したネットワーク記述 (リスト 6 の 2~12 行目) の 2 つに大別できます.

ここでは,こうした結合記述を用いた結合モデルの記述および使用方法を説明します.

リスト6 : NCS 言語による結合記述を用いた Hodgkin-Huxley モデルの記述

1   /*    Hodgkin-Huxley's cell model    */
2   type:          NETWORK;
3   module:        SQUID;
4   cell:          HH[11];
5   gap:           G[10];
6   connection:
7                  HH[0] < ( G[0] < HH[1] );
8                  for( n = 1; n <= 9; n++ ){
9                       HH[n] < (G[n-1] < HH[n-1] + G[n] < HH[n+1]);
10                 }
11                 HH[10] < ( G[9] < HH[9] );
12                 end;
13  /*    HH module    */
14  type:          CELL;
15  module:        HH;
16  exinput:       Iex;
17  input:         Ig;
18  output:        V;
19  observable:    INa, IK, Il, Ig;
20  constant:      VNa = 115.0, VK = -12.0, Vl = 10.6;
21  parameter:     Cm = 1.0, mNa0 = 0.05293, GNa =120.0, GK = 36.0,
22                 Gl = 0.3, hNa0 = 0.5961, nK0 = 0.3177, V0 = 0.;
23  function:
24                 if(V != 25.){               /* sodium current */
25                      am = 0.1 * (25.-v)/(exp((25. - V)/10.) - 1.);}
26                 else{
27                      am = 0.1 * 10.;}
28                 bm = 4. * exp(-V / 18.);
29                 dmNa = am * (1. - mNa) - bm * mNa;
30                 mNa = integral(mNa0,dmNa);
31                 ah = 0.07 * exp(-V / 20.);
32                 bh = 1. / (exp((30. - V) / 10.) + 1.);
33                 dhNa = ah * (1. - hNa) - bh * hNa;
34                 hNa = integral(hNa0,dhNa);
35                 INa = GNa * pow(mNa,3.0) * hNa * (V - VNa);
36                 if(V != 10.){               /* potassium current */
37                       an = 0.01*(10.-V)/(exp((10. - V)/10.) - 1.);}
38                 else{
39                       an = 0.01 * 10.; }
40                 bn = 0.125 * exp(-V / 80.);
41                 dnK = an * (1. - nK) - bn * nK;
42                 nK = integral(nK0,dnK);
43                 IK = GK * pow(nK,4.0) * (V - VK);
44                 Il = Gl * (V - Vl);         /* leakage current */
45                 Iall = Iex - INa - IK - Il + Ig;
46                 dV = Iall / Cm;
47                 V = integral(V0,dV);
48                 end;
49  /*    G module    */
50  type:          GAP;
51  module:        G;
52  input:         VOP(0.1,0);
53  output:        Ig;
54  parameter:     Gl = 5.0;
55  function:
56                 Ig = Gl * (VOP - POSOUT);
57  end;

4.3. 結合モデルの予約語

結合モデル記述の際は,項2.1. 「予約語」 の予約語に加えて,さらに次の 3 つの予約語がシステムにより定義されています.

  1. PRECN

    化学・電気シナプスへの入力値を出力しているセルモジュールを,そのシナプスのシナプス前セルモジュールと呼びます.PRECN はシナプス前セルモジュールのコンポーネント番号を保持しています.化学シナプス,電気シナプスモジュールにおいて適用できます.

  2. POSTCN

    化学・電気シナプスからの出力値を入力しているセルモジュールを,そのシナプスのシナプス後セルモジュールと呼びます.POSTCN はシナプス後セルモジュールのコンポーネント番号を保持しています.化学シナプス,電気シナプスモジュールにおいて適用できます.

  3. POSOUT

    シナプス後セルモジュールからの出力値を保持しています.化学シナプス,電気シナプスモジュールにおいて適用できます.

4.4. 結合モデルにおけるモジュール記述

結合モデルで使用する各モジュールは次項に述べる文により記述します.これらの文の組み合わせはタイプによって決められており,ユーザはそれにしたがって記述を行います.以下には文の詳細とタイプ別記述法について順に説明していきます.

●文

結合記述では,項2.3. 「モジュール記述」 に加えて,さらにいくつかの文を使います.文の詳細を順に説明していきます.

  1. type

    結合記述を用いる場合はモジュールのタイプを宣言しなければなりません.タイプには以下の3 種類があります.

    • CELL : 細胞型 (セルタイプ)

    • SYNAPSE : 化学シナプス型 (シナプスタイプ)

    • GAP : 電気シナプス型 (ギャップタイプ)

    また,結合記述の場合には次のように宣言を行います.

    type : NETWORK

    記述例)

    /*    結合記述の宣言    */
    type : NETWORK;
    /*    セルタイプの宣言    */
    type : CELL;
  2. cell

    モデルで使用するセルモジュールの宣言を行います.モジュール名とコンポーネント数を記述します.複数のセルモジュールを宣言するときは "," で継げます.

    記述例 )

    /*    "HH"というモジュール名のセルモジュールを
          コンポーネント数11個分使用することを宣言    */
    cell :    HH[11];
  3. synapse

    結合記述において,モデルで使用する化学シナプスモジュールの宣言を行います.モジュール名とコンポーネント数で記述します.複数のセルモジュールを宣言するときは "," で継げます.

    記述例)

    /*    "SYN" というモジュール名の化学シナプスモジュールを
          コンポーネント数5個分使用することを宣言    */
    synapse :    SYN[5];
  4. gap

    結合記述において,モデルで使用する電気シナプスモジュールの宣言を行います.モジュール名とコンポーネント数で記述します.複数のセルモジュールを宣言するときは "," で継げます.

    記述例)

    /*    "G" というモジュール名の電気スナプスモジュールを
          コンポーネント数10個分使用することを宣言    */
    gap :    G[10];
  5. input

    モジュールの入力変数名を定義します.セルモジュールの入力変数の数と記述順序は後述する connection 文の結合関係式における記述と一致している必要があります.化学シナプスモジュール,電気シナプスモジュールの入力変数の数は 1 つです.複数のセルモジュールを宣言するときは "," で継げます.また,化学シナプス,電気シナプスの input 文には遅延情報を付加することができます.遅延情報は次の形で記述します.

    input : 変数名 (遅延時間,初期値) ;

    時刻 0 から遅延時間までの間は初期値が出力されます.

    記述例)

    /*    "Ig" を入力変数として定義    */
    input :    Ig;
    /*    "VOP" を入力変数として定義
          ただし,遅延時間0.1で初期値0    */
    input :    VOP(0.1,0);
  6. output

    モジュールの出力変数名を定義します.出力変数の数は 1 つです.

    記述例)

    /*    "V" を出力変数として定義    */
    output :    V;
  7. connection

    コンポーネント同士の結合関係を記述します.これを結合関係式と呼びます.結合関係式は "<" をはさんで,右辺の出力が左辺に入力される構成になります.また,1 つの入力は "()" で囲むことで記述します.したがって,

    mdl[i] < (input1)(input2);

    という書式の場合,input1 がモジュール mdl のコンポーネント i の 1 番目の入力に,input2 がモジュール mdl のコンポーネント i の 2 番目の入力変数に代入されることになります.input 文で宣言されている変数の数と,結合関係式で表される入力の数は等しくなければいけません.つまり,この場合モジュール名 mdl の記述で,

    input : V1, V2;

    のように 2 つの入力変数が宣言されている必要があります.この時,V1 には input1 からの出力が,V2 には input2 からの出力が代入されます.

    ところで,入力をとらない,つまり input 文が存在しないモジュールのコンポーネントへの入力は存在しないことになります.そこで,そのようなセルモジュールのコンポーネントの結合関係式は,次のような書式になります.

    mdl[i] < ();

    また,セルモジュールの入力は,化学または電気シナプスモジュールのコンポーネントを介している必要があります.したがって,セルモジュール mdl のコンポーネント i の 1 番目の入力に,セルモジュール mdl のコンポーネント j からの出力が化学 (または電気) シナプスモジュールのコンポーネント k を通って入る場合の書式は次のようになります.

    mdl[i] < (syn[k] < mdl[j])(input2);

    さらに,各モジュールのコンポーネントへの入力には,モジュールのコンポーネントからの出力を加算したものを指定することも可能です.この場合,"+" で継いで列挙します.例えば,セルモジュール mdl のコンポーネント n の出力がシナプスモジュール syn のコンポーネント m を介したものと,セルモジュール mdl のコンポーネント l の出力がシナプスモジュール syn のコンポーネント o を介したものが結合して,セルモジュール mdl のコンポーネント i の 2 番目の入力となっている場合の記述は次のようになります.

    mdl[i] < (input1)(syn[m] < mdl[n] + syn[o] < mdl[l]);

    結合関係式において右辺に現れるセルモジュールのコンポーネントは,別の結合関係式の最左辺に表れている必要があります.また,化学および電気シナプスモジュールのコンポーネントは単体では用いられず,入力側と出力側に必ずセルモジュールのコンポーネントを持つことに注意してください.

    記述例)

    connection:
    /*    セルモジュール "HH" のコンポーネント1からの出力が,
          シナプスモジュール "G" のコンポーネント0を通って,
          セルモジュール "HH" のコンポーネント0に入れる.    */
    
    HH[0] < (G[0] < HH[1]);
    /*   セルモジュール "HH" のコンポーネント4からの出力が,
         シナプスモジュール "G" のコンポーネント4を通ったものと,
         セルモジュール "HH" のコンポーネント6からの出力が
         シナプスモジュール "G" のコンポーネント5を通ったものと結
         合され,セルモジュール "HH" のコンポーネント5に入力される.  */
    HH[5] < (G[4] < HH[4] + G[5] < HH[6]);
    /*   セルモジュール "HH" のコンポーネント9からの出力が,
         シナプスモジュール "G" のコンポーネント9を通って,
         セルモジュール "HH" のコンポーネント10の1番目の入力になり,
         セルモジュール "CA" のコンポーネント0からの出力が,
         シナプスモジュール "SYN" のコンポーネント0を通って,
         セルモジュール "HH" のコンポーネント10の2番目の入力になる.  */
    HH[10] < (G[9] < HH[9])(SYN[0] < CA[0]);

●セルタイプモジュールの記述

このタイプはシングルモジュールで用いていたものと同じタイプです.セルタイプモジュールで用いられている文は次のようになっています.

  1. type 文

  2. module 文

  3. exinput 文

  4. input 文

  5. output 文

  6. observable 文

  7. constant 文

  8. parameter 文

  9. function 文

3 から 8 は順不同です.また各文の詳細は 項2.3. 「モジュール記述」 および 項4.4. 「結合モデルにおけるモジュール記述」 を参照してください.なお,モジュール記述の終わりには必ず,"end;" をつけることを忘れないでください.

記述例)

/*    celll type module    */
type:        CELL;
module:      CL;
exinput:     Iex;
input:       Ig;
output:      V;
observable:  IK, Il, Ig, nK;
constant:    VK = -12.0, Vl = 10.6;
parameter:   Cm = 1.0, GK = 36.0, Gl = 0.3;
finction:
if(V! = 10.){
     an = 0.01 * (10. - V) / (exp((10. - V) / 10.) - 1);
}else{
     an = 0.01 * 10.;
}
bn = 0.125 * exp(-V / 80.);
dnK = an * (1. - nK) - bn * nK;
nK = integral(0.3177,dnK);
IK = GK * pow(nK,4.) * (V - VK);  /*  カリウム電流  */
Il = Gl * (V - Vl);               /*  リーク電流  */
Iall = Iex - IK - Il + Ig;
dV = Iall / Cm;
V = integral(0.,dV);              /*  膜電位  */
end;

●化学シナプスタイプモジュールの記述

化学シナプスタイプモジュールは入出力が 1 つのモジュールであり,セルタイプのモジュール同士を継ぐ役割を果たします.化学シナプスタイプと電気シナプスタイプの記述上の違いはありません.化学シナプスタイプモジュールで用いられている文と,その並びは次のようになっています.

  1. type 文

  2. module 文

  3. input 文

  4. output 文

  5. observable 文

  6. constant 文

  7. parameter 文

  8. function 文

3 から 7 は順不同です.また各文の詳細は 項2.3. 「モジュール記述」 を参照してください.なお,モジュール記述の終わりには必ず,"end;" をつけることを忘れないでください.

記述例)

/*    synapse type module    */
type:       SYNAPSE;
module:     GABA;
input:      P0(0.1,0);
output:     Tr;
parameter:  FB = 0.01;
function:
Tr = P0/ (FB + P0);
end;

●電気シナプスタイプモジュールの記述

化学シナプスタイプモジュールは入出力が 1 つのモジュールであり,セルタイプのモジュール同士を継ぐ役割を果たします.化学シナプスタイプと電気シナプスタイプの記述上の違いはありません.電気シナプスタイプモジュールで用いられている文と,その並びは次のようになっています.

  1. type 文

  2. module 文

  3. input 文

  4. output 文

  5. observable 文

  6. constant 文

  7. parameter 文

  8. function 文

3 から 7 は順不同です.また各文の詳細は 項2.3. 「モジュール記述」 を参照してください.なお,モジュール記述の終わりには必ず,"end;" をつけることを忘れないでください.

記述例)

/*    gap type module    */
type:        GAP;
module:      G;
input:       VOP(0.1,0);
output:      Ig;
parameter:   GL = 5.0;
function:
Ig = GL * (VOP - POSOUT);
end;

●結合記述

モデルで使用するモジュールを宣言し,各コンポーネントの結合関係を数式的に記述します.結合記述は,他のモジュール記述より先に,つまりモデルファイルの最初に記述するようにして下さい.結合記述で用いられている文と,その並びは次のようになっています.

  1. type 文

  2. module 文

  3. cell 文

  4. synapse 文

  5. gap 文

  6. connection 文

この文の並びは,絶対的なものですので,文を省略するのは構いませんが,順番を入れ替えないように注意してください.また各文の詳細は 項2.3. 「モジュール記述」 を参照してください.尚,結合記述の終わりには必ず "end;" をつけるのを忘れないでください.

結合記述には for 文が用意されており,大幅に記述を簡略化することができます.for 文の書式は次のようになっています.

書式)

for(式1 ; 式2 ; 式3){
      文
}

式 1,2,3 には,代入式,関係式,演算式などの式を 1 つ記述します.式 1 は初期式であり,for ループが始める前に参照されます.そして式 2 が満たされている間,文が実行されます.式 3 は for ループの 1 回のループ毎の最後に行われる式です.この for 文において文を囲む "{","}" は不要です.式 2 には,表 7.3. 「条件式の関係演算子」 に示した関係演算子が使用できます.

記述例)

/*    Hodgkin-Huxley's cell model    */
type:       NETWORK;
module:     SQUID;
cell:       HH[11];
gap:        G[10];
connection:
         HH[0] < (G[0] < HH[0]);
         for( n = 1; n <= 9; n++ ) {
               HH[n] < (G[n-1] < HH[n-1] + G[n] < HH[n+1]);
         }
         HH[10] < (G[9] < HH[9]);
end;

4.5. 神経回路網のモデリング

ここでは,H-H モデルを 図 7.13. 「神経回路モデルの例」 のように接続した神経回路網の記述を例に,そのモデリングについて説明しています.

細胞の結合を考える場合,隣接する細胞からの電流も考える必要があります.

まず,セルタイプモジュールを記述していきます.モジュールのタイプはセルタイプ,モジュール名は "HH" です.

type:         CELL;
module:       HH;

外部入力に注入電流 を,入力に隣接する細胞からの電流 ,そして膜電位 を出力するものとします.

exinput:      Iex;
input:        Ig
output:       V;

ナトリウム電流の反転順位 ,カリウム電流 ,漏れ電流 ,隣接する細胞からの電流 を観測することにします.

observable:   INa, IK, Il, Ig;

各イオン電流の反転電位 を,定数として宣言します.

constant:     VNa = 115.0, VK = -12.0, Vl = 10.6;

膜容量 ,積分の初期値,各イオンの膜コンダクタンス をパラメーターとして宣言します.

parameter:    Cm = 1.0, mNa0 = 0.05293, GNa = 120., GK = 36.,
              Gl = 0.3, hNa0 = 0.5961, nK0 = 0.3177, V0 = 0.;

function 文には,前述した記述を用います.

細胞同士を結合する電気シナプスタイプモジュールの記述は次のようになります.電気シナプスタイプモジュールのモジュール名は "G" です.

type:         GAP;
module:       G;

電圧を入力とし,電流を出力とします.ただし,電圧は初期値が 0 で 0.1 秒の遅れを持って入力されます.

input:        VOP(0.1,0);
output:       Ig;

電気シナプスのコンダクタンスをパラメーターとして宣言します.

parameter:    GL = 5.0;

出力電流 図 7.13. 「神経回路モデルの例」 より次式で表されます.

ここで, は電気シナプス両端の電圧です.これより,function 文は次のように記述されます.

function:
              Ig = GL * (VOP - POSOUT);

最後に結合記述を考えます.このモデルの名前を "SQUID" とします.

type:         NETWORK;
module:       SQUID;

図 7.13. 「神経回路モデルの例」 より,"HH" を 11 個と "G" を 10 個が必要ですので宣言します.

cell:         HH[11];
gap:          G[10];

これらの結合関係を for 文を用いて記述します.

connection:
         HH[0] < (G[0] < HH[0]);
         for( n = 1; n <= 9; n++ ) {
               HH[n] < (G[n-1] < HH[n-1] + G[n] < HH[n+1]);
         }
         HH[10] < (G[9] < HH[9]);

結合記述部をまとめると,次のようになります.

記述例)

/*    Hodgkin-Huxley's cell model    */
type:         NETWORK;
module:       SQUID;
cell:         HH[11];
gap:          G[10];
connection:
         HH[0] < (G[0] < HH[1]);
         for( n = 1; n <= 9; n++ ) {
              HH[n] < (G[n-1] < HH[n-1] + G[n] < HH[n+1]);
         }
         HH[10] < (G[9] < HH[9]);
end;

結合記述は,モデル記述ファイルの先頭に必ず書かれていなければなりません.

3 つの 3 種類のモジュールと結合記述の記述例を全てまとめると,リスト 6 の NCS 言語による神経回路の結合モデル記述が完成します.

4.6. 結合モデルの使用法

作成したモデルファイルを使用し,シミュレーションを行う方法について述べます.ただし,基本的な手順は 項3. 「NCS の使用法」 に記したものと同様ですので,ここでは簡単に示すことにします.関数の書式の詳細などは 項3.4. 「シミュレーション条件の設定」 を参照してください.

さて,はじめにモデルファイルを登録し,プリプロセッサを起動します.

[]satellite[]~/home/demo:[25]% nassign("hh-netowrk")
[]satellite[]~/home/demo:[26]% npp()

または,nassign 関数を用いずに

[]satellite[]~/home/demo:[27]% npp("hh-network")

とすることにより,"hh-network.mdl" をモデルファイルとして登録して,プリプロセッサを起動することになります.

npp 関数実行後,nlink 関数により,C 言語ソースファイルをコンパイル後,NCS のライブラリとリンクされ,実行ファイルが作成されます.

[]satellite[]~/home/demo:[28]% nlink("-02")

npp 関数を実行し,シミュレーション条件ファイル群を生成したら,次に,各種シミュレーション条件を設定します.

時間条件は ntime 関数によって設定します.

[]satellite[]~/home/demo:[29]% ntime(10,0.001,0.01,1)

としたときには,計算時間を 10[ms],計算刻みを 0.001[ms],ストア刻みを 0.01[ms],書き込み刻みを 1[ms] に設定したことになります.

外部入力条件は nstim 関数によって設定します.コンポーネント番号の指定をすれば,設定の範囲内で任意のコンポーネントに入力を行うことができます.

[]satellite[]~/home/demo:[30]% nstim("HH",10,"P",1,0,100,3,999)

とすることにより,モジュール名 "HH" のコンポーネント 10 に対してパルス関数を入力する設定を行います.

出力条件は nout 関数によって設定します.コンポーネント番号の指定により,任意のコンポーネントの出力を得ることができます.

[]satellite[]~/home/demo:[31]% nout(Iin,"HH",10,2 )

によって,モジュール名 "HH" のモジュールのコンポーネント 10 の入力値がバッファ Iin に出力されます.出力値を格納するバッファは nout 関数実行以前に series 型のバッファとして宣言されている必要があります.

また,複数のコンポーネントからの出力をまとめて配列に格納したい時などは,あらかじめ series 型の配列を宣言しておき,nout 関数の第一引数と第二引数との間に新たに配列番号を指定する引数を設けることで実現できます.例えば,

[]satellite[]~/home/demo:[32]% series V[3]
[]satellite[]~/home/demo:[33]% nout(V,0,"HH",0,1)
[]satellite[]~/home/demo:[34]% nout(V,2,"HH",5,1)
[]satellite[]~/home/demo:[35]% nout(V,3,"HH",10,1)

とすることによって,モジュール名 "HH" のモジュールのコンポーネント 0,5,10 の出力値が配列 V[0],V[1],V[2] に出力されます.

●遅延条件の変更

遅延条件は,ndelay 関数によって変更することができます.ndelay 関数の書式は次のようになっています.

書式) ndelay( mdl, var, dt [, init] )

引数)

  1. mdl : モジュール名

  2. var : 内部変数名

  3. dt : ディレイ時間

  4. init : 出力の初期値

モジュール名 "G" の内部変数 "VOP" の遅延条件を変更する場合には次のようになります.

[]satellite[]~/home/demo:[36]% ndelay("G","VOP",0.25)

設定した外部入力条件は nsclist 関数の引数に "D" を用いて表示することができます.

[]satellite[]~/home/demo:[37]% nsclist("D")

DELAY

DELAY INFORMATION
   Data No.,  Input Name,  Delay Time,  Initial Output
     1        G(  VOP)       0.25            0

nsclist 関数を用いることによって,遅延条件が正しく設定されているか確認してください.

電気シナプスや化学シナプスに遅延条件を付加した場合,積分方式にはデフォルトである自動刻み積分法を用いることはできません.したがって,ninteg 関数により,ルンゲ - クッタ - ギル法またはオイラー法を指定する必要があります.

[]satellite[]~/home/demo:[38]% ninteg("R")

●シミュレーションの実行

シミュレーションを実行する前に,パラメータ値が正しく設定されているか nlist 関数で確認しましょう.引数に "ALL" を指定すると,全てのモジュールのパラメーター値を表示します.

[]satellite[]~/home/demo:[39]% nlist("ALL")

Parameter List

Module name :   HH
   Cm       =   [1]
   mNa0     =   [0.05293]
   GNa      =   [120]
   GK       =   [36]
   Gl       =   [0.3]
   hNa0     =   [0.5961]
   nK0      =   [0.3177]
   V0       =   [0]

Module name :   G
   Gl       =   [5]

リスト 7 : 神経回路網モデルのシミュレーション実行用のバッチファイル

1   nassign("hh-network");              # モデルファイルの登録
2   npp();                              # プリプロセッサの起動
3   nlink();                            # 実行ファイルの作成
4   ntime(10.0,0.001,0.01,1);           # 時間条件の設定
5   nstim("HH",10,"P",1,0,100,3,999);   # 外部入力条件の設定
6   series Iin, V[3];
7   nout(Iin,"HH",0,2);                 # 外部出力条件の設定
8   nout(V,0,"HH",0,1);
9   nout(V,1,"HH",5,1);
10  nout(V,2,"HH",10,1);
11  ndelay("G","VOP",0.25);             # 遅延条件の変更
12  ninteg("R");                        # 積分方式の変更
13  ncal();                             # シミュレーションの実行

さて,単一セルタイプモデルと同様に,各種シミュレーション条件を設定後,ncal 関数によりシミュレーションが実行されます.

[]satellite[]~/home/demo:[40]% ncal()

実行後,シェルはコマンド待ち状態になって,シミュレーションが終了したことを示します.

リスト 7 にファイル hh-network.mdl を使ったシミュレーションを行うためのバッチファイル例を示します.このファイルを "hh-network.sl" とするとシミュレーション実行は次のようになります.

[]satellite[]~/home/demo:[41]% inline("hh-network.sl")

シミュレーション結果は,システムモジュール GPM を使用することにより,次のようにしてグラフに表示することができます.

[]satellite[]~/home/demo:[42]% wopen(1,"A4",0,0)
[]satellite[]~/home/demo:[43]% time=(0~1000)/100
[]satellite[]~/home/demo:[44]% graph(V[0],time,0,0,0,0,0)
[]satellite[]~/home/demo:[45]% scale("N","D","N","D")
[]satellite[]~/home/demo:[46]% graph(V[1],time,0,0,0,0,0)
[]satellite[]~/home/demo:[47]% scale("N","D","N","D")
[]satellite[]~/home/demo:[48]% graph(V[2],time,0,0,0,0,0)
[]satellite[]~/home/demo:[49]% title(1,"time [ms]",
                                        "membrane potential [mV]")
[]satellite[]~/home/demo:[50]% axis(1,1,"XY","XY",5,0,0,0,0,0)
神経回路網モデルのシミュレーション結果

図 7.14. 神経回路網モデルのシミュレーション結果

グラフ表示の詳細は,SATELLITE モジュール GPM のマニュアルを参照してください.

Last updated: 2005/11/12