NISS2002演習準備

	MATLABによるシミュレーションとデータ解析の基礎

１．MATLAB入門

MATLABでは，行ベクトルは

>>x = [ 1 2]

のように空白，あるいはカンマでで区切って，列ベクトルは

>>y = [ 3; 4]

のうようにセミコロンで区切って表します．行列は，

>>A = [ 1 2; 3 4]

のように表します．その要素は，例えば1行目の第2列なら

>>A(1,2)

のように表します．

MATLABでは，+, - などの演算子，sin(), exp() などの関数は，行列の要素ごとに
作用するようになっているので，膨大なデータ処理もforループなど使わずに記述
することができます．例外は，*, /, ^で，

>>A*B

はAとBの行列としての積，

>>A^2

は行列AA，

>>A^-1

はAの逆行列を表します．

>>A/B

はA=xBの解，つまりA*B^-1を表します．同様にAx=Bの解は，A\Bで求まります．

要素ごとのかけ算，わり算やべき乗をしたいときは，

>>A.*B
>>A./B
>>A.^3

のように"."つきの演算子を使います．


２．MATLABで微分方程式を解こう

ではさっそく，MATLABでシミュレーションを走らせてみましょう．MATLABには，
ode...()という常微分方程式の初期値問題を解くための関数がいくつか用意されて
います．

状態空間のあるところで急に感度が高くなるような系を"stiff"な系と呼びますが，
スパイクニューロンモデルもその類に入るので，ode23s(), ode15s()など，"s"つ
きの関数を使います．

(1) １次微分方程式

まず最初に，簡単な１次の微分方程式 dx/dt = ax を解いてみましょう．

>>type first.m

function xdot = first( t, x)
% First: simple first-order dynamics

a = -1;
xdot = a*x;

>>ode23s( 'first', [0 10], 1)

おなじみの指数関数が表示されたでしょうか？ここで'first'は微分方程式の右辺
が書いてある関数名，[0 10]は開始と終了時刻，1はxの初期値です．係数aや初期
値を適当に変えて試してみてください．

(2) 多次元微分方程式

多変数の微分方程式の場合，状態変数は列ベクトルで表すしきたりになっています．

>>type second.m

function xdot = second( t, x)
% Second: simple second-order dynamics

A = [ 0 -1; 1 0];
xdot = A*x;

>>ode23s( 'second', [0 10], [1;0])

列ベクトル[1;0]は状態変数xの初期値です．青がx(1), 緑がx(2)の波形です．

こうやってグラフを描くだけではなくて，その結果をいろいろ使いたいという場合，

>>[t,x] = ode23s( 'second', [0 10], [1;0])

というように，時刻と変数の時系列を行列に代入します．例えば，

>>plot( x(:,1), x(:,2))

とすると，２次元の状態空間での軌跡を見ることができます．ここでx(:,1)は，行
列xの第1列のすべての要素，という意味です．

ode...()はシミュレーションの時間刻みを自動調節してくれるので，時刻のベクト
ルtの中身は必ずしも一定間隔ではありません．もし等間隔のデータが欲しい，と
いう場合は，

>>[t,x] = ode23s( 'second', 0:0.1:10, [1;0])

とすると，時刻0.1間隔に補間し直したデータを出力してくれます．

MATLABでは，時系列は縦方向に並べる慣例になっていて，行列xの中で各時刻の状
態変数は行ベクトルとして並んでいます．

さっき状態変数は列ベクトルで表すと言ったばかりなのに一貫性がないじゃないか，
と不審に思う方もいるかと思いますが，その通り，一貫性がないのです．気をつけ
ましょう．

(3) 入力波形の扱い

システムの入力に対する応答を見たい場合，入力波形を関数の中で読み込む必要が
あります．ここでは，global変数を使った方法を示します．

>>type alpha.m

function [xdot,xinit,option] = alpha( t, x, flag)
% Alpha: second order response

global stim
switch flag
case ''
	s = interp1(stim(:,1),stim(:,2),t);	% linear interpolation
	a = 1;
	xdot = [ s - a*x(1);
		x(1) - a*x(2)];
case 'init'
	xdot = [stim(1,1) stim(end,1)];	% time span
	xinit = [0;0];		% initial condition
	option = odeset('MaxStep',0.1);	% not to skip the stimulus
end

一般に関数の中の変数はその場限りのものとして扱われますが，ここでは(時刻,入
力値)のデータを入れた'stim'という行列を参照するため，"global stim"と宣言し
ています．

flagが空の場合は，interp1でstimの値を線形補完した入力値を計算し，それに応
じたxdotの値を返します．flagが'init'の場合は，xdotの代わりに開始と終了時刻
を，さらにxの初期値と計算上のオプションを返すようになっています．

例えば，時刻t=1でのインパルスを

>>global stim
>>stim = [0 0; 0.9 0; 1 1; 1.1 0; 10 0]

と定義し，

>>ode15s( 'alpha')

とすると，ode15s()がまずalpha('','','init')を呼んで初期設定値を取り込んで
から走り始めます．

x(2)の波形は俗に"alpha function"と呼ばれる応答で，シナプス電流のモデルなど
によく使われます．


３．Hodgkin-Huxleyニューロンモデル

ではいよいよ，ニューロンモデルのシミュレーションに入りましょう．まずは古典
中の古典，Hodgkin and Huxley (1952)のイカ軸索モデルから．

>>type hh.m

function [xdot,xinit,option] = hh( t, x, flag)
% HH: Hodgkin-Huxley (1952) model of squid axon

global	stim
switch flag
case ''
	% state variables
	v=x(1);		% membrane potential (mV)
	m=x(2);		% sodium current activation
	h=x(3);		% sodium current inactivation
	n=x(4);		% potassium current activation
	% input current (uA/cm^2)
	Ie = interp1(stim(:,1),stim(:,2),t);	% linear interpolation

	% membrane capacitance (uF/cm^2)
	C=1;
	% maximum conductances (uS/cm^2)
	gna = 120;		% sodium
	gk = 36;		% potassium
	gl = 0.3;		% leak
	% reversal potentials (mV)
	Ena = 50;		% sodium
	Ek = -77;		% potassium
	El = -54.4;		% leak
	% activation and inactivation rates
	am = 0.1*(v+40)/(1-exp(-(v+40)/10));
	bm = 4*exp(-(v+65)/18);
	ah = 0.07*exp(-(v+65)/20);
	bh = 1/(1+exp(-(v+35)/10));
	an = 0.01*(v+55)/(1-exp(-(v+55)/10));
	bn = 0.125*exp(-(v+65)/80);
	% xdot
	xdot = [ (Ie - gna*m^3*h*(v-Ena) - gk*n^4*(v-Ek) - gl*(v-El))/C;
		am*(1-m)-bm*m;
		ah*(1-h)-bh*h;
		an*(1-n)-bn*n];
case 'init'			% Initialization
	xdot = stim(1,1):0.1:stim(end,1);	% time span
	xinit = [ -65; 0.1; 0.3; 0.3];	% initial v,m,h,n
	option = odeset('MaxStep',10);	% option
end

変数xやパラメタの値を，それらしい名前の変数に代入しています．これは必ずし
も必要というわけではないですが，プログラムを読みやすく，後でバグを発見しや
すくするためです．

では刺激電流Ieを，20msまでゼロ，その後100msで100microA/cm^2まで増やしてい
くというシミュレーションをしてみましょう．

>>global stim
>>stim = [0 0; 20 0; 100 100]
>>ode15s( 'hh')

標準のグラフだと膜電位と他の変数のスケールが揃っていなくて見にくいですが，

>>[t,x] = ode15s( 'hh');
>>subplot(2,1,1); plot(t,x(:,1)); subplot(2,1,2); plot(t,x(:,2:4));

とすると，２つに分けて見ることができます．

では，刺激の波形や細胞のパラメタを変えたときにどうなるか，いろいろ試してみ
てください．

（続く）

****