Go to Qiita Advent Calendar 2024 Top
LoginSignup
7
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@gkhtnk(Gaku Hatanaka)

（神経科学で学ぶ）matlab 初心者向け講座 STA編 2/2

Last updated at Posted at 2018-09-01

はじめに

この講座ではmatlab初心者向けに解析で使ういくつかの機能について、神経科学に関係する内容に触れながらその使い方を解説します。

この投稿は関連する前半での実装をもとにしています。前半がお済みでない方はまずはそちらから先にお読みください。

モデルニューロンの応答の計算

前半でGaborフィルタの実装が終わったので、モデルニューロンの応答を計算していきましょう。

モデルの定義

まず初めに、単純型細胞のモデルとして、Linear-Nonlinear(LN) Cascadeタイプのモデルを想定します。

時刻$t$におけるスパイクの発火の有無を以下のように定義します。今回は話は簡単にするために、モデルニューロンは時刻$t_i$の刺激$\mathbf{s} (t_i)$に対して発火するか、しないかの２状態のみを取るものとします。

\begin{aligned}
p(t) &= \mathbf{k}^T \mathbf{s}(t) & \text{Linear process} \\
q(t) &= f\left( p(t) \right) & \text{Nonlinear process} \\
r(t) &\sim P \left( q(t) \mid \Theta \right) & \text{Spike generation}
\end{aligned}

数式を目の前にすると身構えてしまう人もいるかもしれませんが、今回はかなり単純な定義を採用します。1

単純型細胞モデルの実装

$p(t) = \mathbf{k}^T \mathbf{s}(t)$は内積です。直感的な説明をすると、ここではフィルタと刺激がどの程度似ているのかを見積もります。今回は$\mathbf{k}$としてGaborフィルタを採用します。

2次元のフィルタ行列K1(size: [16 x 16])と3次元の刺激行列NS(size: [16 x 16 x 2048])は1次元目（行、y軸）と2次元目(列、x軸)の長さがともに等しいので「基本的な演算で互換性のある配列サイズ」となっています。そのため「暗黙的な拡張」が適用されて、計算結果は今回の場合NSと同じサイズ(size: [16 x 16 x 2048])の3次元行列になります。

matlabでは要約統計量の計算は基本的に若い次元（1次元目、２次元目、、、）から大きさが1より大きい次元に適用されます。K1.*NSの計算結果に対して2回sumを計算していますが、それぞれが1次元目、2次元目について合計を計算することになります。丁寧な記述を心がけるのであれば次元を指定する引数に値を渡しましょう。squeeze関数は大きさが1の次元を削除します。このため、計算結果P1は1次元の配列(size: [2048 x 1])になります。

% prepare a set of noise stimuli
NS = randn(16, 16, 2048);

% make 2D Gabor filter as a model (kernel) of neuron
K1 = make2DGabor(16, 0.25, 0, pi/6, 2, 1/4, 0);

% projection of the stimuli NS by kernel K1 to 1D space
P1 = squeeze(sum(sum(K1.*NS)));

$q(t) = f\left( p(t) \right)$は非線形処理です。モデル自体を内積だけの式 $ r(t) = \mathbf{k}^T \mathbf{s}(t) $ だけで定義してもよいのですが、非線形関数$f(\cdot)$を入れておくと、処理が都合よくなることがあるので使っておきます。今回は$f(\cdot)$として半波整流max(0, x)を仮定します。半波整流はスパイクが負の数を表現できないということを表現します。

% rectify the projected stimuli
Q1 = max(0, P1);

$r(t) \sim P \left( q(t) \mid \Theta \right)$は細胞の発火が適当なパラメータ$\Theta$のもとで内部変数$q(t)$に依存して確率的に定まるということを意味します。真面目に考える際は、細胞の発火が何らかの確率分布に従うという条件を置くこともあります。今回は$q(t)$がある閾値$\theta \in \Theta$をこえた場合に必ず発火するという仮定を置きます。この部分は非線形処理$f(\cdot)$の中で記述されることもあります。2

% generate a spike responce by threshold
R1 = 2 < Q1;

ここまでで、各処理段階での結果を得ることができました。3

ここまでの実装では、各処理をそれぞれ丁寧に実装していました。ただ、複数のモデルニューロンを扱う際などは、毎回同じ処理を記述することは間違いのもとですし、コードの見通しが悪くなってしまいます。だからと言ってわざわざ関数ファイルを作ってまで実装する必要性もそこまで感じられません。そんな中途半端な処理を実装する方法として、「無名関数」を使った実装と、スクリプト内関数を使った実装を以下で紹介します。

(補足) 無名関数を使った処理の実装

それぞれの処理を「無名関数」を用いて実装してみましょう。

無名関数はわざわざ関数ファイルを作るほどではない程度の処理を実装する際に使用します。funcname = @(input) implementationの形で記述します。4

% project: projection of the stimuli NS by kernel K1 to 1D space
project = @(Kernel, Stimuli) squeeze(sum(sum(Kernel.*Stimuli)));

% rectify: nonlinear process
% ramp function a.k.a ReLU(Rectified Linear Unit)
rectify = @(projected) max(0, projected);

% isfired: spike generation
isfired = @(rectified, threshold) threshold < rectified;

P1 = project(K1, NS);
Q1 = rectify(P1);
R1 = isfired(Q1, 2);

(補足) サブ関数による実装例

functionキーワードを用いた関数は1つのファイル内に複数定義することができます。5

参考までに同じ実装をサブ関数での実装例を紹介しておきます。関数の記述は普通の関数の場合と同じです。サブ関数はスクリプトファイル内のみから実行可能で、スクリプトや関数の実装に続けて記述します。ファイル内では必ず末尾に実装します。

... % script implementation

P1 = project(K1, NS);
Q1 = rectify(P1);
R1 = isfired(Q1, 2);

... % script implementation

function P = project(K, S)
% project: projection of the stimuli by kernel to 1D space
% K: Kernel
% S: Stimuli
P = squeeze(sum(sum(K.*S)));
end

function Q = rectify(P)
% rectify: nonlinear process
% ramp function a.k.a ReLU(Rectified Linear Unit)
% P: projected stimuli
Q = max(0, P);
end

function R = isfired(Q, T)
% isfired: spike generation
% R: rectified (and projected) stimuli
% T: threshold
R = T < R;
end

時系列データの可視化

それぞれの処理段階での結果の時系列を表示してみましょう。

グラフが見やすいように先ほど計算したデータの一部のみを表示するように設定します。
結果が格納される配列に対して、インデックス配列を適用することで、指定した部分のみを取り出すことができます。

plot関数を用いることで時系列データを表示します。
入力引数でグラフの色やスタイルを変更することができます。

% define time axis
t = 1:256;

figure();

subplot(3, 1, 1);
plot(t, P1(t));

subplot(3, 1, 2);
plot(t, Q1(t), '.');

subplot(3, 1, 3);
plot(t, R1(t), 'r:*');

Fig4_1.PNG

画像ファイルの保存

解析などの実行結果を画像ファイルとして出力しましょう。今回はSTAによる受容野の推定結果を題材に開設します。

Spike Triggered Averageの実行

STAによる受容野の推定結果を計算します。

STAの基本的な考え方は、Spikeを誘発した刺激の取り出し、その平均像を可視化するといったところです。ノイズ刺激の時間軸（3次元目）に対してニューロンが発火したかどうかを表すR1を論理インデックスとして取り出してきます。平均の計算はmean関数を利用し、応答を誘発したノイズ刺激データの3次元目に対して適用することで各ピクセルごとに平均像を計算することになります。

subplotを使って画像を並べることでモデルのフィルタK1と推定結果E1を比較してみましょう。発火数によりますが、受容野のおおよその形状は再現できていると思います。

% Estimate model receptive field by STA
E1 = mean(NS(:, :, R1), 3);

subplot(1, 2, 1);
imagesc(K1);
axis image off;
title('Model');

subplot(1, 2, 2);
imagesc(E1);
axis image off;
title('Estimation');

Fig5_1.PNG

matlabには画像を保存する関数がいくつかあります。必要な状況に合わせて使い分けましょう。ここではその中からprint関数とimwrite関数の２つを紹介します。6

print関数によるFigureの保存

print関数は表示しているFigureの描画を画像ファイルとして保存します。

画像ファイルの形式や解像度などをオプションで細かく設定できるのでドキュメントをお読みください。

print('resultSTA_print', '-dpng');

resultSTA_print.png

imwrite関数による行列データの画像化

imwrite関数は、行列データ自体を画像として保存します。

関数が様々なデータ形式、出力ファイル形式に対応している都合上、仕様に沿った記述を行わないと思った通りの出力が得られません。とりあえず何も考えずに推定結果E1をimwriteの引数に入れてみましょう。7

imwrite(E1, 'resultSTA_imwrite.png');

resultSTA_imwrite.png
グレースケールかつ、明るかった部分（ON領域）しか見えなくなってしまいました。
これはimwrite関数が行列データを[0 ~ 1]のデータであると認識しているためです。

mat2gray関数を利用してデータの最小値を0、最大値を1とするように正規化します。
引数

imwrite(mat2gray(E1), 'resultSTA_imwrite_normalized.png');

resultSTA_imwrite_normalized.png

カラーマップを反映させたいときはgray2ind関数とind2rgb関数を利用します。
gray2ind関数はグレースケール画像をn段階のインデックス画像に変換します。
ind2rgb関数はインデックス画像にカラーマップを適用し、RGB画像に変換します。

imwrite(ind2rgb(gray2ind(mat2gray(E1), 256), parula(256)), 'resultSTA_imwrite_normalized_parula.png');

imwrite関数自体もインデックス画像を入力に受けることができるため、以下のようにしても同様の結果を得ることができます。

imwrite(gray2ind(mat2gray(E1), 256)), parula(256)), 'resultSTA_imwrite_normalized_parula.png');

resultSTA_imwrite_normalized_parula.png

動画ファイルでの保存

実行結果を動画ファイルとして保存しましょう。

時間変化の可視化

ディスプレイ上でデータを可視化する以上、可視化は2次元に制限されます。8

時間とともに表示を変化させることで、表示する軸を1つ拡張することができます。
ここでは時間とともに鮮明となるSTAの推定結果を可視化してみます。tiの値をforで順に増やしていきます。9

titleを使用することで、各座標軸ごとにグラフのタイトルをつけることができます。さらにここでは時間の変化とともに、タイトルの文字列を更新しています。num2strは数値を文字列に変換します。大かっこ[ ]を使用することで文字列を結合することができます。ここでは['Stimulus: ', num2str(ti)]とすることで、Stimulus: 512のように何番目のノイズ画像を処理したのかを表示します。

各ループの最後にpauseで指定した時間(秒)だけ処理を停止します。こうすることで、描画がコマ送りに更新されます。10

実行を途中で停止したいときはCtrl + Cで処理を中断できます。

% for-loop to update STA ewsult with time
figure();

for ti = 1:2048
  P1 = project(K1, NS(:, :, 1:ti));
  Q1 = rectify(P1);
  R1 = isfired(Q1, 2);
  E1 = mean(NS(:, :, R1), 3);
  
  subplot(1, 2, 1);
  imagesc(NS(:, :, ti));
  axis image off;
  title(['Stimulus: ', num2str(ti)]);
  
  subplot(1, 2, 2);
  imagesc(E1);
  axis image off;
  title('Estimation');
  
  pause(1/30);
end

image.png

getframeによるFigure画像の取得

getframeを用いることで、現在描画されているFigureをデータとして確保することができます。
とりあえず以下を実行して、各時刻tiにおける刺激と、STAの結果を取得します。

% for-loop to update STA ewsult with time
figure();
Frame(2048) = struct('cdata',[],'colormap',[]);

for ti = 1:2048
  P1 = project(K1, NS(:, :, 1:ti));
  Q1 = rectify(P1);
  R1 = isfired(Q1, 2);
  E1 = mean(NS(:, :, R1), 3);
  
  subplot(1, 2, 1);
  imagesc(NS(:, :, ti));
  axis image off;
  title(['Stimulus: ', num2str(ti)]);
  
  subplot(1, 2, 2);
  imagesc(E1);
  axis image off;
  title('Estimation');
  
  drawnow();
  Frame(ti) = getframe(gcf);
end

VideoWriterによる動画の書き込み

VideoWriterを利用することで、動画ファイルを作製することができます。

これまでと処理の書き方が少々異なりますが、以下を実行することでMP4形式の動画ファイルとして書き込むことができます。11

% Create VideoWrite Object
vidObj = VideoWriter('resultSTA.mp4', 'MPEG-4');

% Set properties bofore file opend
vidObj.FrameRate = 30;
vidObj.Quality = 100;

% Writeing video file, do not forget open and close
vidObj.open();
vidObj.writeVideo(Frame);
vidObj.close();

残念ながらQiitaではmp4ファイルを貼りつけることができないようです。

imwriteによるgif動画の書き込み

また、imwriteを使用することにより、gif動画ファイルを作成することもできます。12

filename = 'resultSTA.gif';
DelayTime = 1/30;

for ti = 1:2048
  [A, map] = rgb2ind(frame2im(Frame(ti)), 256);
  
  if ti == 1
    % For first frame
    imwrite(A, map, filename, 'gif', 'LoopCount', Inf, 'DelayTime', DelayTime);
  else
    imwrite(A, map, filename, 'gif', 'WriteMode','append', 'DelayTime', DelayTime);
  end
end

resultSTA_4.gif

お疲れ様です。
とりあえず今回は以上です。13

前半へ戻る

  1. 今回の目的はSTAの実装ではなく、matlabの練習ですのでかなり粗い記述となっています。この辺りの定義は論文によってかなり違いがありますので調べてみてください。前半でも紹介しました、Meyer AF et al. 2016 Front Syst Neurosci、辺りが参考になるかと思います。 ↩

  2. ここでは$\Theta$を固定したパラメータとして扱いましたが、これは必ずしも必須ではありません。実際的には、刺激に対するスパイクの発生がポアソン分布などの何らかの確率分布に従うという仮説が置かれます。また、今回のモデルニューロンでは半波整流したのちに閾値処理をおこなっています。この場合、結果からみると半波整流は不要ですがモデルを変更するときのことを考えて分けて考えます。 ↩

  3. それぞれの処理結果が変数P1,Q1,R1へと代入されています。このような変数名の付け方はよくありません、変数の名前は後から読んでもわかりやすいものにしましょう。適当なアルファベットや単語を略したものを付けがちですが、できるだけ避けましょう。コードが長くなった時に同じ名前で変数を再定義し、値を上書きしてしまうという悲劇を招きます。 ↩

  4. 「無名関数」は「関数ハンドル」という変数として扱われます。関数に対して関数を引数として渡したいときなどに役に立ちます。 ↩

  5. スクリプト内サブ関数はR2016b以降からサポートされています。それ以前のバージョンをお使いの方は，この実装は使用できません。 ↩

  6. そのほかにもsaveas関数や、savefig関数などを利用することで、図を保存することができます。 ↩

  7. ここでは出力された画像データを拡大して貼り付けています。上記スクリプトを実行して得られるデータは縦横16 pixなので、画像ファイルとしてはかなり小さなものになります。 ↩

  8. 3次元データを可視化する関数もありますが、それらはあくまで2次元に投影した図になります。 ↩

  9. tiが順に増えていく際に、論理インデックス配列となるR1の大きさはNSの3次元目の要素数より少ないですが、「より小さな配列の論理インデックス」はインデックス配列の欠損要素がゼロに設定されるので、エラーが出ずに実行できます。 ↩

  10. ここで停止時間として、1/30を使用していますが、実際にはデータの処理と描画にかかる時間があるため、必ずしも描画が30　Hzでおこなわれるということを意味しません。また、各時間における描画をギリギリまで早くするにはdrawnowを使用します。 ↩

  11. VideoWriterはクラスという形式で、定義されています。今回はクラスについて扱わないので詳細は割愛します。 ↩

  12. ここで表示しているgif動画はアップロードするファイル容量の関係で、4フレームおきに描画した結果になっているため、ここでのスクリプトの4倍のスピードで描画されています。 ↩

  13. 今回の記事を書くにあたり、いくつかの内容を削りました。また気が向いたら関連する内容を追加しようと思います。 ↩

7
5
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
7
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?