日本神経回路学会 2009 年 9 月

動的スパイク相関の状態空間モデル

State-space Model of Dynamic Spike Correlations
島崎 秀昭 (PY)† , 甘利 俊一 † , Emery N. Brown‡ , Sonja Grün†
Hideaki Shimazaki(PY)† , Shun-ichi Amari† , Emery N. Brown‡ , and Sonja Grün†† 理化学研究所脳科学総合研究センター‡Massachusetts General Hospital, Depart. of Brain and Cognitive Sciences MITHarvard-MIT Division of Health Sciences and Technology, Cambridge MA USA
shimazaki@brain.riken.jp
Abstract— To trace active assembly activities ofneurons from multiple neural spike data, we developed
a method for estimating dynamic spike correlations bymeans of a state-space analysis. Discretized parallel
spike sequences are modeled as a non-stationary multivariate Bernoulli process using a log-linear function (amultinomial logit), which introduces states of higherorder interaction terms. We constructed a recursive
Bayesian filter/smoother to estimate time-varying loglinear parameters, and an optimization method ofhyper-parameters under the likelihood principle.

枠組みにより厳密に定義できる [Amari, IEEE Trans Inf

Keywords— State-space model, Log-linear model,
Information Geometry

N 個の同時スパイク時系列を多変量ベルヌーイ過程とみなし, 多変量ロジット連結関数 (対数線型モデル) を用いてモデル化する. そのために, n 回の試行から成るスパイク時系列を幅 ∆ の T 個の区間 (ビン) に区切る.
t,lXt,l = (X1t,l , X2t,l , . . . , XN) を l 試行目, t 番目のビンの二値データ (1, スパイク有り; 0, スパイクなし) とし, 生成モデルとして以下の対数線型モデルを考える.
xi (i = 1, . . . , N ) を二値変数, x = (x1 , x2 , . . . , xN ) を
N2 個のスパイクパターンを表すベクトルとする. 対数線型モデルでは, 同時分布を因子 Fi (x) = xi , Fij (x) =xi xj , Fijk (x) = xi xj xk (i < j < k), etc を用いて

1

序論
神経細胞集団の協調活動が神経系の情報処理に本質的な役割があるとする仮説が古くから提唱されている[Hebb, 1949] . この仮説に従えば, 複数細胞の協調的活動が同期的な発火活動という形で現れ, 観測されたスパ
イク時系列間の同期スパイクの構造に外界刺激や脳の内部状態 (注意や期待) の情報が符号化されている可能性がある. 実際, 多細胞同時計測技術の発展とともに, 神経細胞間のスパイク同期が動物の知覚・行動と密接な関係にあることを示唆するいくつかの実験的証拠が提出

Theory, 2001]. また状態空間モデルは一般に非定常システムの同定に用いられる. 我々は状態 (対数線型モデルの自然母数) の事後分布を対数二次形式で近似し, 非線形再帰濾波公式を導出, 濾波公式と固定幅平滑化公式を併せて時間変動する相関構造を推定する手法を開発した. これにより集団としての神経細胞活動の動的構造を明らかにし, 動物の認知・行動との関係について神経生理学的知見を得ることを目標としている.
2方法

2.1

対数線型モデルlog p (x | θ t ) = −ψ(θ t ) +されている [e.g. Riehle et al., Science, 1997].
神経細胞の同期的協調活動を調べるためには, 同時記と表す.

Xi∈{Ω1 ,...,Ωr }θit Fi (x) (1)tt
ここで, θ t = [θ1t , . . . , θ12, . . . , θ1···r, . . .]0 は録されたスパイク時系列間のスパイク相関構造を推定す

自然母数で, 時刻 t の同時分布を定める.

る必要がある. ただし, ３つ以上の神経細胞の協調活動正規化項.

を調べるためには, ２体間の相関だけでなく一般に多体部分集合 (Ω1 = {1, 2, . . . , N } , Ω2 = {12, 13, . . .},間の相関 (高次相関) も推定する必要があるかもしれなΩ3 = {123, 124, . . .}, etc). r ≤ N であり, r < N のとき式 (1) は最大エントロピーモデルと呼ばれる. 特に

い. また協調活動が動的であるとすれば, ダイナミックな相関構造を追跡できなければならない. 同時に, 覚醒・ψ(θ t ) はまた, Ωk は N 個から k 個を取り出した行動下で得られるデータを扱うためには, これらの推定

θi (i ∈ Ωr ) は i に含まれる二値変数同士の r 次の周辺相関を表す. 同期スパイク確率を表す ηit = E [Fi (x) |θ t ]は非定常なスパイク発火率の下で行える必要がある.

は期待値パラメータと呼ばれ, θ t と双対な座標系 η t =スパイク発火率と相関を同時に推定する技術を開発した.

tt[η1t , . . . , η12, . . . , η1···r, . . .]0 を構成する. 期待値パラメー¡¢Pntタの有効推定量は yi = n1 l=1 Fi Xt,l で与えられ,ttyt = [y1t , . . . , y12, . . . , y1···r, . . .]0 は r 次の最大エントロ

対数線形モデルを用いることで, 高次相関を情報幾何のピーモデルの十分統計量である.

この問題に答えるために, 我々は対数線形モデルを観測モデルとする状態空間モデルを用いて, 時間変動する

(a)Underlying Rates and Pairwise Correlations

Neuron 128Firing Rate
166 [Hz]37465000.1

0.2

0.3

0.4

0.5
Time [s](b)Simulated Spike Sequences

Positively Correlated

Neuron 2

Neuron 3

Neuron 4

Neuron 6

Neuron 7

Neuron 850 Trials

Neuron 1

Pairwise correlation

00.5

Time [s]Neuron 5(c)Negatively Correlated

Estimated Rates and Pairwise Correlations from 200 trials

図 1: 同時スパイク時系列データ (8 個) の動的相関構造の推定. (a) 二次の相関までを含む時間依存対数線型モデ
ル. 0.1 秒刻みのスナップショット. 二値変数をノードとし, 発火率をその濃淡で表す. 二次相関の強さはリンクの太さ及び濃淡で表す. 実線 (赤) 及び破線 (青) はそれぞれ正及び負の相関. (b) 生成モデルに従う同時スパイク時系列データ (50 試行). (c) 状態空間モデルによる動的発火率・相関構造の同時推定. 推定には 200 試行を用いた.

2.2

ただし, θ t|t−1 , Wt|t−1 は式 (4) より与えられる一期先状態空間モデル

全観測データを y1:T = {y1 , y2 , . . . , yT }, 背後の確率

予測平均及び共分散行列. 濾波分布は θ t に関する複変動を θ 1:T = {θ 1 , θ 2 , . . . , θ T } で表すと, 観測モデルは雑な式となるが, 対数凸関数の性質を保存する. そこ

p (y1:T |θ 1:T ) =TYexp [n (yt0 θ t − ψt )] .

で対数二次近似 (ラプラス近似) を行い, 事後平均 θ t|t(2)t=12で与えられる. 状態方程式として, 以下の式を仮定する

θ t = Fθ t−1 + ξ t ,ことである. このために, 一期先予測・濾波分布を逐次的に求めた後, 固定幅平滑化公式により平滑化分布を得た. 超パラメータは周辺尤度最大化を基準に選び, EMアルゴリズムに従って事後分布の獲得と最尤推定を交互に行うことで最適化した. ここでは, このうち最も重要な濾波分布の導出のみを記す.
濾波分布は Bayes の定理により, 尤度と一期先予測分布で表される:

p (yt |θ t , y1:t−1 ) p (θ t |y1:t−1 )p (yt |y1:t−1 )∝ exp[n (yt0 θ t − ψt )−−[ ∂θ∂t ∂θ0 log p (θ t |y1:t ) |θt =θt|t ]−1 として, 以下の非線形t濾波公式を得る:(3)(t = 2, . . . , T ). ここで ξt は, 平均 0・共分散 Q の多変量正規分布に従う乱数. 初期値は θ1 ∼ N (µ, Σ). 超パラメータをまとめて w = [F, Q, µ, Σ] とする.
我々の目的は平滑化事後分布 p (θ t |y1:T , w) を求める

p (θ t |y1:t ) =∂を ∂θ
log p (θ t |y1:t ) = 0 の解, 事後共分散行列 Wt|t をt¢¢0 −1 ¡1¡
θ t − θ t|t−1 ]. (4)θ t − θ t|t−1 Wt|t−12

θ t|t = θ t|t−1 + nWt|t−1 (yt − η t|t ),−1−1
Wt|t= Wt|t−1+ nGt .

(5)
(6)

ここで, キュムラント母関数 ψ(θ t ) の一階微分は双対座
標 η t を, 二階微分は Fisher 情報量行列 Gt を与えることに注意. 式 (5) は Newton-Raphson 法を用いて解ける.

3

結果
発火レート及び 2 次の相関が変動する 8 変数からなる対数線形モデル (図 1a) からスパイク時系列を生成した (図 1b). スパイクデータから背後の動的な発火率・2次相関構造を推定した (図 1c). 同様にして, 高次相関をモデルに組み込み, 推定することも出来る [1].
参考文献

[1] H. Shimazaki, S. Amari, E. N. Brown, and S.
Gruen “State-space Analysis on Time-varying Correlations in Parallel Spike Sequences” Proc. IEEEICASSP2009, 3501-3504.