はじめに

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

今回はSpike Triggered Average(STA)による受容野の測定を取り上げます。その過程で簡単な視覚野の神経細胞のモデルの作成、仮想視覚刺激呈示によるスパイク列の生成などを行うことで、matlabの基本的な機能に触れていきたいと思います。
長くなったので前半と後半に分けています。

matlabを全く触ったことの無い人はまず、「MATLAB入門」から始めましょう。matlabが用意している関数にはできるだけ、mathworks社のサイトへリンクを貼っています「関数 - カテゴリ別」。公式のドキュメントが整っているという点がmatlabの長所の1つです。使い方の分からない関数はdocまたはhelpを使うことで調べることができます。

% Reference page in Help browser
doc doc;

% Help for functions in Command Window
help help

実装例は整理した後で、githubにでもあげる予定です。²

Spike Triggered Average(STA)とは

視覚野の神経細胞の受容野を測定する手法の1つです。

この講座では、matlabの基本的な使い方の解説に主眼を置いています。受容野やSTAについて厳密な議論はせず、できるけ単純なモデルを使った説明を努めます。ざっくりとした概要を理解したい方は下記の阪大・大澤先生のページを参考にしてください
「受容野の定量的測定法」。³

2次元データの可視化

モデルの実装に入る前にmatlabの練習を1つ、おこないましょう。

matlabのデータの1番基本になるものはその名の通り行列です。行列形式で与えられたデータを可視化する方法を試してみましょう。使用するデータは何でも良いのですが、ここではrand関数で一様乱数を生成します。
figureで新規ウィンドウを生成し、imagescでデータを表示します。⁴

% Create 4x4 matrix by random
M = rand(64);

% Create new figure
figure();

% Display matrix M
imagesc(M);

座標軸の調節にはaxisを使います。ここではアスペクト比の調節image, 軸の不可視化offオプションを用います。データを描画する際に不必要な情報はできるだけ削りましょう。

% Adjust the aspect ratio and hide the axis
axis image off;

デフォルトのカラーマップはparula が指定されています。⁵

カラーマップはcolormapで変更することができます。例えば黒から白へ変化するグレースケールはgrayで使用することができます。

% Change the colormap to gray
colormap gray

colormapの色と値の対応はcolorbarを使うことで確認できます。

% Display colorbar
colorbar;

スクリプトファイルの作成

Gaborフィルタの計算を題材にスクリプトファイルを作成してみましょう。

コマンドラインからeditコマンドを用いて以下を実行します。ここではscript2DGaborという名前のファイルを作成していますが、各自で好きな名前を付けても構いません。このような形式で保存されたファイルはスクリプトファイルと呼ばれます「プログラミングとスクリプト」。一連の手続きをスクリプトファイルに記述し、保存しておくことで再び同じ処理を実行することができます。

% Create a script file
edit script2DGabor

今回は簡単のため、時間方向は無視して、空間方向（XY平面）のみで構成される受容野を考えます。正弦波とガウシアンの掛け合わせで表現されるGaborフィルタを仮定します。

格子点行列の生成

まずは古典的にmeshgridを使って格子点行列を生成します。⁶

linspaceを使って[-1 ~ 1]のベクトルを入力として与えます。出力X, Yはそれぞれの軸の格子点の値に対応します。subplotは1つのウィンドウに複数の描画を並べる際に活躍します。

% Create X, Y grid by [-1 ~ 1] N=64 vector
[X, Y] = meshgrid(linspace(-1, 1, 64));

figure();

subplot(2, 3, 1);
imagesc(X);
axis image off;

subplot(2, 3, 2);
imagesc(Y);
axis image off;

鋭い人なら気が付いたかもしれませんが、Yの増加する方向が上下逆転していると思います。matlabのimage表示関数はデータを上から下へと描画するように設定されています。⁷

座標軸の回転

方位を定義して座標軸を回転します。

Gaborフィルタの重要なパラメータの1つに方位ORがあります。方位の定義（0度の定義, 増加する回転方向）は研究グループによってバラバラです。先ほど生成した格子点に対して回転ベクトルをかけて座標軸を回転させます。ここでは縦じまを0度とし、方位の値が大きくなるにつれ、時計周りに回転するように定義します。

円周率の値としてpiを使用します。matlabではいくつかの定数が定義されています。これらの値は上書きできます。同じ名前の変数は定義しないようにしましょう「定数とテスト行列」。

ここで注意しなければならないのはmatlabでの演算の扱いです。matlabは線形代数の規則に従う演算（行列演算: A*B）と要素単位で行う演算（配列演算: A.*B）を区別します「配列と行列の演算」。matlabでは、いわゆる「暗黙的な拡張」によりエラーが発生せずに、意図しない結果が得られることがあります「基本的な演算で互換性のある配列サイズ」。今回、ＸＹ軸の格子点データX, Yに対して回転ベクトルの係数（スカラー）cos(OR), sin(OR)を掛ける際、配列演算であることを強調するために*ではなくあえて.*を用います。⁸

% Define OR: orientation (30 degree)
OR = pi/6;

% Rotated axis
A = cos(OR).*X + sin(OR).*Y;

subplot(2, 3, 3);
imagesc(A);
axis image off;

2次元の正弦波の生成

座標軸の回転が定義出来たら三角関数cosを適用して正弦波を作成します。

matlabの関数のほとんどは行列の入力に対応しており、たいていの場合、行列の各要素に演算を適用した結果を返します。パラメータとして空間周波数SFを定義しましょう。

% Define SF: spatial-frequency (2 cycle per image)
SF = 2;

% Calculate Sinusoidal wave
Sinus = cos(2*pi*SF.*A);

subplot(2, 3, 4);
imagesc(Sinus);
axis image off;

2次元のガウシアンの生成

受容野は空間的に限局します。ガウシアンを用いることでその様子を近似します。

まず関数hypotを使って中心からの距離を定義する行列Rを定めます。正規分布関数の定義では、係数 $ \frac{1}{\sqrt{2\pi\sigma^2}} $ を掛けますが、ここでは省きます。

% Define SG: spatial envelope (1/4 of image)
SG = 1/2;

% Radius from center
R = hypot(X, Y);

% Calculate Gaussian envelope
Gauss = exp( -R.^2./(2*SG^2) );

subplot(2, 3, 5);
imagesc(Gauss);
axis image off;

Gaborフィルタの生成

2つの行列を要素ごとに掛け合わせ、Gaborフィルタを生成します。

先ほども述べた通り、matlabでの行列演算と配列演算の違いに注意しましょう。特に今回の場合、同じサイズの正方行列を利用するため、どちらの書き方でもエラーは出ませんが、異なる計算結果が出力されます。

% Calculate Gabor-filter by multiplication of Sinusoidal wave and Gaussian
Gabor = Sinus.*Gauss;

subplot(2, 3, 6);
imagesc(Gabor);
axis image off;

スクリプトファイルの実行

一度ワークスペースを空にした後スクリプトファイルを実行してみましょう。

ワークスペースはclearを使うことでリセットできます。
また、ワークスペース上の変数を表示するにはwhosを利用します。

% Remove all the items from workspace
clear all;

% excute scriptfile
script2DGabor;

% List variables in workspace, with sizes and types
whos;

関数ファイルの作成

関数ファイルの作成することで、処理を引数ともに実行できます。

初期視覚野神経細胞のモデルとしてGaborフィルタを定義したわけですが、パラメータとして設定した方位ORや空間周波数SFなどは個々の細胞ごとに異なります。これらのパラメータを入力引数として持ち、Gaborフィルタを出力するような関数を考えます。

先ほど作成したGaborフィルタを計算するスクリプトファイルをもとに順に関数を作成していきます。

functionキーワードを使い、ファイル名を関数名と一致させることで、関数ファイルを作成することができます。ここでは関数名、ファイル名共にmake2DGaborを使用します。

コマンドライン上で以下を実行し、新規ファイルを作成します。

%Create function file
edit make2DGabor

ファイルの拡張子からはスクリプトファイルと関数ファイルの違いはわかりません。
関数ファイルでは以下のように１番初めの実行可能行でfunctionキーワードを用いて関数の宣言を行う必要があります。

make2DGabor.m

function Gabor = make2DGabor ...

入力引数を整理

はじめに入力引数を整理しましょう。

先ほどのスクリプトファイルをよく読んでどの変数を入力引数にすべきかを考えます。ここでは追加として、出力するフィルタのサイズSZ, 受容野の中心位置を定めるパラメータCX, CY, フィルタの位相PHを追加しておきます。

make2DGabor.m

function Gabor = make2DGabor(SZ, CX, CY, OR, SF, SG, PH) ...

ヘルプテキストの作成

functionのすぐ下にコメントを挿入することで、ヘルプテキストを作成することができます。

ヘルプテキストはhelpを使用すると参照できます。

help make2DGabor

関数の動作、入力引数・出力引数の説明などを記述しておきます。

make2DGabor.m

function Gabor = make2DGabor(SZ, CX, CY, OR, SF, SG, PH)
% make2DGabor make 2D Gabor filter
%   Gabor = make2DGabor(SZ, CX, CY, OR, SF, SG, PH)
%   SZ    : Size of output filter
%   CX    : Center position X
%   CY    : Center position Y
%   OR    : Orientation
%   SF    : Spatial-Frequency
%   SG    : Spatial envelope
%   PH    : Phase
%   
%   To visualize with imagesc function, use below option
%   axis image off
...

関数の実装

先ほど作成したスクリプトファイルの内容を元に関数の処理を実装していきます。

出力引数と同名の変数を定義し、そこに値を代入することで、関数の出力を得ることができます。(ここでは変数Gaborが出力引数になります)

今回はファイルの先頭から処理を記述していきましたが、関数を実装していく中で、新しい入力引数が必要になることは多々あります。先に実装の見通しを立てておくことは役立ちますが、必ずしもこの順番で関数を作っていく必要はありません。

make2DGabor.m

function Gabor = make2DGabor(SZ, CX, CY, OR, SF, SG, PH)
% make2DGabor make 2D Gabor filter
%   Gabor = make2DGabor(SZ, CX, CY, OR, SF, SG, PH)
%   SZ    : Size of output filter
%   CX    : Center position X
%   CY    : Center position Y
%   OR    : Orientation
%   SF    : Spatial-Frequency
%   SG    : Spatial envelope
%   PH    : Phase
%   
%   To visualize with imagesc function, use below option
%   axis image off

% Create X, Y grid by [-1 ~ 1] axis vector
[X, Y] = meshgrid(linspace(-1, 1, SZ));

% Rotated axis
A = cos(OR).*(X-CX) + sin(OR).*(Y-CY);

% Calculate Sinusoidal wave
Sinus = cos(2*pi*SF.*A - PH);

% Radius from center
R = hypot(X-CX, Y-CY);

% Calculate Gaussian envelope
Gauss = exp( -R.^2./(2*SG^2) );

% Calculate Gabor-filter by multiplication of Sinusoidal wave and Gaussian envelope
Gabor = Sinus.*Gauss;

様々なパラメータを持つGaborフィルタの可視化

関数を使用することで、異なるパラメータを持つGaborフィルタを簡単に生成することができます。

自分の実装した関数の挙動が正しいか、それぞれの引数に適当な値を入れて確かめましょう。

figure();
colormap parula(1024);

subplot(2, 3, 1);
imagesc(make2DGabor(64, 0, 0, 0, 2, 1/4, 0));
axis image off;
caxis([-1, 1]);

subplot(2, 3, 2);
imagesc(make2DGabor(128, 0.25, 0, 0, 2, 1/4, 0));
axis image off;
caxis([-1, 1]);

subplot(2, 3, 3);
imagesc(make2DGabor(128, 0.25, 2/3, 0, 2, 1/4, 0));
axis image off;
caxis([-1, 1]);

subplot(2, 3, 4);
imagesc(make2DGabor(128, 0.25, 2/3, pi/4, 2, 1/4, 0));
axis image off;
caxis([-1, 1]);

subplot(2, 3, 5);
imagesc(make2DGabor(128, 0.25, 2/3, pi/4, 4, 1/4, 0));
axis image off;
caxis([-1, 1]);

subplot(2, 3, 6);
imagesc(make2DGabor(128, 0.25, 2/3, pi/4, 4, 1/8, 0));
axis image off;
caxis([-1, 1]);

後半へ続く

神経科学と言いましたが、今回の話が関係するのはそのごくごく小さなシステム神経科学という分野でのお話になります。もともと研究室向けに書いていた記事になりますのでご了承ください。 ↩
大した量にはならないと思います。コード部分を順番にコピペして実行できるようになると思います。エラー等が発生しましたら教えてください。 ↩
詳しく勉強したい方は以下のレビューがフリーでアクセスできて周辺も話題も含めて俯瞰できるように思います。 ↩

Meyer AF, Williamson RS, Linden JF, Sahani M. (2016)
Models of Neuronal Stimulus-Response Functions: Elaboration, Estimation, and Evaluation.
Front Syst Neurosci 10:109.
https://dx.doi.org/10.3389/fnsys.2016.00109
PMID: 28127278
matlabには行列を表示する関数としてimshow、image、imagescの３つが用意されています。imagescは自動的に行列の値の最大・最小値をもとに色の配置を決めます。とりあえず行列の中身を見たいときはこれがお勧めです。描画設定をいじれば細かい見た目は変更できるので好みのものを使ってください。 ↩
このカラーマップの利点は以下の記事を読んでください。
* A New Colormap for MATLAB – Part 1 – Introduction
* A New Colormap for MATLAB – Part 2 – Troubles with Rainbows
* A New Colormap for MATLAB – Part 3 – Some Reactions
* A New Colormap for MATLAB – Part 4 – The Name ↩
暗黙的な拡張を利用すれば、格子点行列を作る必要はなくなります。可読性が悪くなるためここでは素直にmeshgridを利用します。 ↩
これは画像データがコンピュータ上でどのように扱われているのかに関係します。表示の向きを変えるには以下のような方法があります。
* axis xyオプションを使う「axis y 軸の方向」
* AxesオブジェクトのYDirプロパティを'normal'に設定する「Axes のプロパティ」
* flipudで行列の上下を逆転する
* 気にせず他の場所で都合を合わせる。 ↩

非常に細かいところですが、ここで実装しているのは点の回転移動ではなく座標軸の回転です。
回転後の座標軸として用いている変数Aは以下の$X'$に相当します。
これを取り違えると回転の方向が逆転します。 ↩

座標軸の回転

\begin{pmatrix}
X' \\
Y'
\end{pmatrix}
= 
\begin{pmatrix}
cos{\theta} & sin{\theta} \\
-sin{\theta} & cos{\theta}
\end{pmatrix}
\begin{pmatrix}
X \\
Y
\end{pmatrix}

点の回転移動

\begin{pmatrix}
x' \\
y'
\end{pmatrix}
= 
\begin{pmatrix}
cos{\theta} & -sin{\theta} \\
sin{\theta} & cos{\theta}
\end{pmatrix}
\begin{pmatrix}
x \\
y
\end{pmatrix}

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

はじめに

Spike Triggered Average(STA)とは

2次元データの可視化

スクリプトファイルの作成

格子点行列の生成

座標軸の回転

2次元の正弦波の生成

2次元のガウシアンの生成

Gaborフィルタの生成

スクリプトファイルの実行

関数ファイルの作成

入力引数を整理

ヘルプテキストの作成

関数の実装

様々なパラメータを持つGaborフィルタの可視化