特集「時空間統計解析：新たなる分野横断的展開」［原著論文］

統計数理（2012）
第 60 巻 第 1 号 149–157
c 2012 統計数理研究所


単一試行生体イメージングデータ解析のためのAR 型モデルに基づく時空間フィルタリング法
三分一 史和

1,2

3

4

5

・越久 仁敬 ・岡田 泰昌 ・川合 成治 ・
田村 義保

1,2

・石黒 真木夫6,7（受付 2011 年 7 月 1 日；改訂 9 月 8 日；採択 9 月 27 日）

要

旨

回帰分析や相互相関解析は脳の動的イメージングデータにおける神経賦活の検出に広く用いられている．これらの解析法には神経賦活の時間的変化を反映した参照関数を先験的に仮定する必要があり，言い換えれば，神経賦活のうち，その変動パターンが参照関数と相似性の高いものが検出されるということになる．また，参照関数が定義出来ない場合は，解析そのものが困難になってしまうという問題がある．我々の先行研究において，これらの問題点を回避するため，時空間フィルタリング法を提案した．この方法は繰り返し計測（多重試行）
された呼吸関連データにおいて，神経賦活に対応した動力学的状態変化が生じる時刻と場所の検出を可能とした．しかし，この方法は自励活動などの単一試行データには直接適用することができない．
そこで，本研究では，先行研究における時空間フィルタリング法に sliding time window を組み込むことにより改変し，単一試行データにおいても神経賦活の検出が可能であることを示した．
キーワード： 時空間フィルタリング，イノベーションアプローチ，脳機能マッピン
グ，光イメージング．
1.

はじめに

ここ 10–20 年間の脳神経活動の計測技術の進歩により，膜電位イメージングなどの光イメー
ジングや機能的磁気共鳴画像（fMRI）などのように，脳活動を時空間データとして計測することが可能となった．イメージングデータは特に高空間分解能で神経ダイナミクスの多くの空間情報を得られることが特徴である．これらのデータは 2 次元もしくは 3 次元空間に配置された各グリッドに時系列が対応した構造を持つ．そして，時空間データから統計的に有意な神経賦活を検出するために多くの方法が提案されてきた．その中でも，一般線形モデルの枠組みを用いた Statistical Parametric Mapping が fMRI データ解析において標準的に用いられている
1
2
3
4
5
6
7

統計数理研究所：〒190–8562 東京都立川市緑町 10–3
総合研究大学院大学 複合科学研究科統計科学専攻
兵庫医科大学 生理学講座：〒663–8501 兵庫県西宮市武庫川町 1–1
独立行政法人国立病院機構 村山医療センター：〒208–0011 東京都武蔵村山市学園 2–37–1
統計数理研究所 外来研究員：〒190–8562 東京都立川市緑町 10–3
統計数理研究所 名誉教授：〒190–8562 東京都立川市緑町 10–3
総合研究大学院大学 名誉教授

統計数理

150

第 60 巻

第 1 号 2012

（Friston et al., 1995）．
光イメージングデータは fMRI データに比べ空間の次元が 1 つ低いだけで類似の構造を持っているが，これまで標準的な解析方法は提案されていない．このような状況において，Oku etal.（2007, 2008）と Okada et al.（2007）は time-lagged 相互相関解析を多重計測において加算平均により事前処理した光イメージングデータに適用し，ラットやカエルの脳幹における呼吸リ
ズムの生成メカニズムやパターン形成過程を明らかにした．彼らは，横隔膜の吸息活動を反映する第 4 頚髄前根からの出力信号を参照関数として用いている．解析は相互相関係数の有意性を統計検定し，有意性が認められた部位をマッピングするという流れで，視覚的にも非常に分かりやすいものであるが，いくつかの問題点がある．相互相関解析あるいは回帰分析は参照関数と与えられた時系列との形状的な類似性を評価するものであり，参照関数と異なるパターンをもつ神経賦活は検出することが困難である．さらに，そもそも参照関数を定義することができないデータの場合には解析そのものが行えない．
このような問題点を回避するために，我々はイノベーションアプローチによる時空間フィル
ターと統計検定を組み合わせた方法を開発した（Miwakeichi et al., 2011）．この方法では，各々のピクセルに対応する時系列の神経賦活が生じていない時間区間において近傍のピクセルからの影響も考慮した Nearest-Neighbor AR（NNAR）モデルを同定し，残りの時間区間をその同定された NNAR モデルのパラメーターを用いてフィルタリングしてイノベーション（フィルター
出力）を推定する．この方法では参照関数を先験的に仮定する必要が無く，NNAR モデルを同定した区間とフィルタリングを施した区間のイノベーションの強度を比較することにより時系
列の動力学的性質の変化，つまり神経賦活が生じた時刻と部位を検出することが可能となる．
先行研究では，イノベーションの統計検定のためには，ある程度のサンプル数（試行回数）が必要であるため，適用の範囲は外部刺激により誘発される反応や，呼吸などの周期的な自励活動など繰り返し計測が可能なデータに限定される．そこで，本研究では単一試行のデータの解析のために先行研究で提案した方法の改良を行い，その適用性を議論する．
2.

方法

2.1 データ計測
新生ラットの脳幹において，呼吸リズムを生成する部位として para-facial respiratory group（pFRG）と pre-Bötzinger complex（preBötC）と呼ばれる 2 つの領域の存在が一般的に知られている（図 1）．本研究ではこれらの部位をイメージングデータの計測領域に設定し，吸息活動は第 4頚髄前根（C4VR）
からの信号でモニターした．C4VR からの信号はサンプリング周波数 1kHz で計測し，15Hz–3kHz の帯域でバンドパスフィルターと全波整流を施し，さらに減衰時定数 100ms
で積分処理を行った．記録用チャンバーは蛍光顕微鏡（MVX-10， Olympus Optical，Tokyo，Japan）に設置し，試料は電位感受性色素（di-2-ANEPEQ，Invitrogen，Carlsbad，CA，USA）
にて染色した．励起光の光源にはタングステン ハロゲンランプを用い（150W，λ = 480 − 550 nm），
long-pass barrier filter（λ > 590 nm）を通過した落射蛍光は CMOS センサーアレイ（MiCAM
によりサンプUltima L-camera，BrainVision；100 µm × 100 µm pixel size，100×100 pixel array）
リング周波数 50Hz で計測した．C4VR からのアナログ信号は A–D 変換の後，光イメージン
グデータと共に記録媒体に格納した．アナログ C4VR 信号は，オーディオモニターで呼吸活動と判定された最小振幅以上のものをウインドウディスクリミネータで選択し，それに対応する Transistor-Transistor Logic（TTL）パルスを生成させ，そのパルスは光計測システムを駆動するトリガーに用いた．記録したタイムフレーム 1024 で，記録はトリガー信号の 768 タイム
フレーム前から開始するよう設定した．そして，この計測を 30 回繰り返し，比較的アーチファ
クトの重畳が少ない 29 回のデータを解析に用いた．

単一試行生体イメージングデータ解析のための AR 型モデルに基づく時空間フィルタリング法

151

図 1. （a）新生ラットを腹側から見た様子とおよその脳幹の位置を示す模式図，（b）新生ラッ
トの脳幹と計測領域（正方形内）を示す模式図，（c）第 4 頚髄前根（C4VR）出力信号のオ
リジナルデータ（上），積分処理された C4VR 出力信号（下）．

各々のピクセルにおいて，タイムフレーム毎に最初の時刻における蛍光強度
（F0 ）
からの変化
量
（∆F）を計算し，試料内における細胞膜に沈着した色素量と照射光強度のばらつきを補正するため，各ピクセルにおける背景蛍光強度はその最大値で正規化した．そして，∆F と正規化された背景蛍光強度（F）の比率（∆F/F）は各々のタイムフレームにおいてピクセル毎に算出した．F の値には閾値を設定し，F< 0.25 であるピクセルにおいては，∆F/F= 0 とした．また，色素の退色に起因する背景の線形トレンドは一次方程式をフィッティングすることにより除去した．
2.2 時空間フィルタリング法
時系列 η(t), t = 1, . . . , S において，AR モデルは過去の p 点の値の線形結合として以下のよう
に表される
(2.1)

η(t) − β =p

α(i)η(t − i) + ε(t) .

i=1

ここで，ε(t) は予測誤差，β は定数，p はモデルの次元， α(i) は AR 係数である．時系列の線
形動力学的性質はパラメーターベクトル ϑ = {α(1), . . . , α(p), β, σε2 } により同定される，ここで，σε2 は ε(t) の分散である．イノベーションは時系列の残差から以下のように推定できる
 p


(2.2)
ε(t) = η(t) − η̂(t) = η(t) −α(i)η(t − i) + β .
i=1

そして，線形 AR プロセスで推定できない信号はイノベーションの中に残されることになり，
それは AR プロセスのフィルター出力と解釈できる．
この線形 AR プロセスに外場からの入力 γ(t), t = 1, . . . , S, を仮定すると，AR モデルは外生変数型 AR モデル
（ARX）へと一般化することができる(2.3)

η(t) − β =

p

i=1

α(i)η(t − i) +q
ϕ(j)γ(t − j) + ε(t) .

j=1

光イメージングデータの場合，観測点は 2 次元平面上のピクセルに対応し，あるピクセル座標
ν = (l, m) に対して近傍のピクセルからの影響のみを考慮すると，その制約された ARX モデル

統計数理

152

第 60 巻

第 1 号 2012は以下のようになる，η ν (t) − β ν =

(2.4)

p


αν (i)η ν (t − i) +

i=1

q
 
δ υ (j)η υ (t − j) + εν (t) ,υ∈N(ν) j=1

ここで， N(ν) はピクセル ν = (l, m) の近傍のピクセル座標の集合である．近傍のピクセルをν に対して辺で接しているものに限定すると，それらの座標はN(ν) = {(l + 1, m), (l, m − 1), (l, m + 1), (l − 1, m)}(2.5)となる．この制限された ARX モデルを Neighbor Autoregressive model（NNAR）と呼ぶことにし，ピクセル ν に対するイノベーションはε(t)ν = η ν (t) − η̂ ν (t)
 p

q
 ν
 
ν
ν
υ
υ
ν
α (i)η (t − i) +
δ (j)η (t − j) + β= η (t) −

(2.6)

i=1

υ∈N(ν) j=1と計算される．
この NNAR モデルを区間 t1 ≤ t ≤ t2 において同定し，その同定されたモデルで任意の区間

t1 ≤ t ≤ t2 にフィルタリングを施すことができ，各々の区間に対応するイノベーションを ε1,ν (t)と ε2,ν (t) とする．フィルタリングを施した区間においてシステムの動力学的特性が変化すると，それは同定された NNAR モデルでは予測ができない信号として出現し，その時刻 t において ε2,ν (t ) の振幅レベルが増加することになる．その増加がモデルを同定した区間のイノベー
ションの振幅レベルと比較して有意なものかどうかを評価する必要があるが，時刻 t に対応するイノベーションは 1 つしかサンプルを持たないので，直接統計検定を行うことはできない．
そこで，w + 1（w は偶整数）の幅をもつ sliding time window t − (w/2) ≤ t ≤ t + (w/2) を導入
t +(w/2)
2,ν し，window 内のイノベーションの平均値 (w + 1)−1 n=t(t ) とモデルを同定した区 −(w/2) ε

t2
−1
1,ν
間全体 t1 ≤ t ≤ t2 のイノベーションの平均値 (t2 − t1 + 1)(t) の差の統計検定を行
t=t1 εう．本研究では t 検定を用いた．そして，時刻 t をシフトすることにより，t 値の時間情報を得ることが可能となり，さらに，t 値をマップ化（t-map）
することにより動力学的特性が変化した部位を検出することができる．
3.

結果

各ピクセルに対応する時系列において，呼吸関連の神経賦活が開始される前の区間において，背景の定常振動を NNAR モデルで同定した．本研究では計算コストを考慮し，対象とするピクセルに辺で接する 4 つの近傍のピクセルを NNAR モデルに用い，パラメーターは最小二乗法で推定した．次に，この同定したモデルを用いて時系列の残りの区間にフィルタリングを施しイノベーションを推定した．時間軸の原点（0s）は C4VR からの信号において観測される吸息活動の立ち上がり時刻とし，NNAR モデルは −4.22s から −2.24s の区間で同定し（100time frames），その後の −2.22s から 5.24s の区間（374 time frames）にフィルタリングを施した．
NNAR モデルの定義式（2.4）におけるモデルの次元 p と q は先行研究と同様 p = q = 7 に設定した（Miwakeichi et al., 2011）．Sliding time window の幅と神経賦活が生じる時刻の不確定さは
トレードオフの関係にあり，幅を小さくすると時間精度は上がるが， t 検定の信頼性は低下する．そこで，パラメトリック検定において一般的に必要とされるサンプル数が 30 であるので，
window 幅の下限を w + 1 = 31 とし，そして，幅を 61，91 を変えたときの解析結果を比較した
ところ，検出された部位についてはほぼ同じであったため，本研究では w + 1 = 31 を window 幅

単一試行生体イメージングデータ解析のための AR 型モデルに基づく時空間フィルタリング法

図 2.

153

帯域毎の全ピクセルの平均パワースペクトルの空間強度分布．（a）オリジナルデータ，（b）イノベーション時系列，（c）全ての全ピクセルの平均パワースペクトル（オリジナル
データ（上），イノベーション（下））．円内は呼吸関連の神経組織の部位を示している．
なお，パワースペクトルは最大値と最小値を [1 0] の区間に合わせるように規格化した．

に採用した．そして，NNAR モデルの推定とフィルタリングの手順を 29 回分の多重試行デー
タ全てに適用した．
まず最初に，背景信号の基本的な振動特性を調べるため，FFT によりオリジナルデータとイノベーション時系列のパワースペトルを計算した．図 2（c）はフィルタリングを施した区間
における全てのピクセルのオリジナルデータとイノベーション時系列のパワースペクトルの平均値を示しており，オリジナルデータにおいては 0.1–0.5 Hz，3.0–6.0 Hz，12.0–14.0 Hz に対応，
する主な 3 つのピークが認められる．各々のピークに対応するパワーの分布マップは図 2（a）（b）に示されている．イノベーション時系列に対応するパワースペクトルの平均値ではこれらの周波数帯域におけるパワーが低減されており，分布マップにおいては，特段に強度の高い部
位は認められず，全体としてフラットになっている．これは，NNAR モデルにより背景振動が効率よく同定されたということを示している．
0.1–0.5 Hz に対応するピークはオリジナルデータとイノベーション時系列の両方に認められ，分布マップにおいては，高い強度をもつピクセルのグループ（円内）
はほぼ呼吸関連部位と対応している．呼吸関連活動はフィルタリングを施した区間で生じるので，システムの動力学的特性変化としてイノベーション時系列に予測不能な信号として出現したといえる．
NNAR モデルで同定した区間のイノベーションとフィルタリングを施した区間のイノベー

統計数理

154

図 3.

第 60 巻

第 1 号 2012

t-map の時間変化を示す代表的なタイムフレーム．
（a）オリジナルデータ，（b）イノベーション時系列．
（c）
C4VR からの出力時系列とタイムフレームに対応する時刻点．
（t-mapの閾値は p<0.05 で設定し False Discovery Rate（FDR）により補正した．また，5 ピ
クセルサイズ未満のクラスターは除去してある．点 A∼ B は t 値の時間変動を調べた
ピクセルの位置に対応．）

ションの強度レベルの平均値は t 検定で評価し，全てのピクセルにおいてこの検定を繰り返した．そして，t 値の空間分布である t-map を作成し，有意にシステムの動力学的特性変化が生じる時刻と部位を調べた．図 3（b）は代表的な 5 枚のタイムフレームにおける t-map を示している（t 値の閾値を p < 0.05 の水準で設定し，False Discovery Rate で補正したところ対応する t 値
はオリジナルデータとイノベーション時系列の場合で各々 t > 4.1950 と t > 1.98 となった．オリジナルデータとイノベーション時系列に対応する t-map を同じ水準で比較するため，閾値にはより厳しい t > 4.1950 を採用した）．呼吸関連の神経賦活は吻側延髄腹外側（RVLM）（Onimaruに対応する pFRG の尾側から始まり，吻側に広がり，そして preBötC の尾側et al., 1988, 2003）
へと伝搬していることが認められる．図 3（a）はイノベーション時系列の代わりにオリジナル
データを用いて行った t 検定に対応する t-map である．これらの t-map は同じレベルの閾値を設定しているが，イノベーション時系列のものより有意性を示した領域は小さくなっている．
図 4 はより強く神経賦活が認められたピクセルの中から場所を変えて選んだ 4 点
（A∼D（
）図 3）
における t 値の時間変動を示している．点 A においては，オリジナルデータとイノベーション
時系列を用いた解析の両方で似たような変動を示し，C4VR で検出される吸息活動のオンセッ
トから約 0.28s 後に有意水準を超えている．点 B においては，イノベーション時系列を用いた解析でのみ有意水準を超え，この近傍での神経賦活は吸息活動のオンセットの約 0.1s 前から生じている．点 C ではオリジナルデータとイノベーション時系列を用いた解析結果の差が特に顕著に見られ，吸息活動のオンセットの約 0.1s 後にクリアな神経賦活のピークが認められる．点 D においては有意水準を超えたのはイノベーション時系列を用いた解析結果のみであり，pFRG と preBötC に相当する部位よりも神経賦活は遅く吸息活動のオンセットの約 0.8s 後にピークが認められる．

単一試行生体イメージングデータ解析のための AR 型モデルに基づく時空間フィルタリング法

図 4.

4.

155

図 3 の点 A∼D に対応するピクセルにおける t 値の時間変動．実線はイノベーション
時系列，破線はオリジナルデータに対応する t 値の時間変動を示す．水平の点線は図 3の t-map の閾値に対応（t =4.1950）．

考察

本研究で解析した膜電位イメージングデータのように S/N 比が 0.01–0.3％と極めて低いデー
タ（Ruangkittisakul et al., 2009）から信号を検出するには，S/N 比を向上させるため多重試行と加算平均が必要である． Oku et al.（2007, 2008）
や Okada et al.（2007）
が解析に用いた time-lagged
相互相関解析も加算平均により事前処理されたデータに適用している．しかし，この方法は単一試行データには適用できなく，先行研究（Miwakeichi et al., 2011）において提案した時空間フィルタリング法も統計検定は多重試行を前提としているので直接適用することはできない．そこで，本研究においては，NNAR モデルにより推定されたイノベーションの統計検定に sliding
time window を用いた方法を提案し， 29 回の試行データ各々で神経賦活の検出を試みたところ，約 40％の試行において呼吸関連部位である pFRG と preBötC の検出に成功した．これは，NNAR モデルが背景の定常振動や近傍のピクセルからの散乱光の影響を同定し，フィルターア
ウトプットであるイノベーション時系列からこれらが効率良く除去されていることを示している．図 2 において，3.0–6.0 Hz と 12.0–14.0 Hz の 2 つのピークに対応するパワーの分布マップ
では，特に尾側の方に強度が高い部位が認められるが，これは生理学的に呼吸関連の部位とは認められないので，計測システムの機械的振動や試料を in vivo 状態にしておくための灌流液の乱れに起因するものと思われる．
今回用いた sliding time window は時間軸上でのスムージング操作に対応するのであるが，
図 3 と図 4 で示したように，イノベーション時系列の代わりにオリジナルデータを用いた解析
では，呼吸関連部位を十分に検出することができていない．これは，アーチファクトに起因す
る背景の振動に神経活動が埋もれてしまい，sliding time window のみでは効率的な神経賦活の検出は十分に行えず，NNAR モデルによる時空間フィルタリングが重要な役割を担っているということが言える．
本研究では呼吸関連の多重試行データをベンチマークデータとして用い，各々の試行を単一
試行データとして扱った．多重試行データであっても，呼吸関連の神経賦活は自励活動であり，

156

統計数理

第 60 巻

第 1 号 2012

外部刺激による誘発活動よりは S/N 比が低いので，このデータに基づいて検証された本研究で提案した方法は，同じ自励活動でかつ多重試行が困難な睡眠脳波やてんかん性の脳神経活動な
どの脳の自励活動の解明にも有効であることが示唆される．

参 考 文 献Friston, K. J., Holmes, A. P., Worsley, K. J., Poline, J. B., Frith, C. D. and Frackowiak, R. S. J.

（1995）. Statistical parametric maps in functional imaging: A general linear approach, Human（4）, 189–210.
Brain Mapping, 2
Miwakeichi, F., Oku, Y., Okada, Y., Kawai, S., Tamura, Y. and Ishiguro, M.（2011）. Detection andvisualization method of dynamic state transition for biological spatio-temporal imaging data,（3）, 859–866.
IEEE Transaction on Medical Imaging, 30
Okada, Y., Masumiya, H., Tamura Y. and Oku, Y.（2007）. Respiratory and metabolic acidosis differentially aﬀect the respiratory neuronal network in the ventral medulla of neonatal rats,（10）, 2834–2843.
European Journal of Neuroscience, 26
Oku, Y., Masumiya, H. and Okada, Y.（2007）. Postnatal developmental changes in activation proﬁlesof the respiratory neuronal network in the rat ventral medulla, The Journal of Physiology,（1）, 175–186.
585. Spatiotemporal organization of frog
Oku, Y., Kimura, N., Masumiya, H. and Okada, Y.（2008）
respiratory neurons visualized on the ventral medullary surface, Respiratory Physiology &Neurobiology, 161（3）, 281–290.
. A novel functional neuron group for respiratory rhythm generOnimaru, H. and Homma, I.（2003）
ation in the ventral medulla, The Journal of Neuroscience, 23（4）, 1478–1486.
Onimaru, H., Arata, A. and Homma, I.（1988）. Primary respiratory rhythm generator in the medulla（2）, 314–324.
of brainstem-spinal cord preparation from newborn rat, Brain Research, 445
Ruangkittisakul, A., Okada, Y., Oku, Y., Koshiya, N. and Ballanyi, K.（2009）. Fluorescence imaging（1–2）, 26–38.
of active respiratory networks, Respiratory Physiology and Neurobiology, 168

Proceedings of the Institute of Statistical Mathematics Vol. 60, No. 1, 149–157 (2012)

157

A Spatio-temporal Filtering Method Based on an AR Type Model for
Analysis of Single-trial Biological Imaging Data
Fumikazu Miwakeichi1,2, Yoshitaka Oku3 , Yasumasa Okada4 , Shigeharu Kawai5 ,Yoshiyasu Tamura1,2 and Makio Ishiguro6,7
1 The Institute of Statistical Mathematics
2 Department of Statistical Science, School of Multidisciplinary Sciences, The Graduate University for

Advanced Studies
3 Department of Physiology, Hyogo College of Medicine
4 Laboratory of Electrophysiology, Clinical Research Center, Murayama Medical Center
5 The Institute of Statistical Mathematics, Visiting Research Fellow
6 The Institute of Statistical Mathematics, Professor Emeritus
7 The Graduate University for Advanced Studies, Professor Emeritus

Regression and cross correlation analyses have been widely used to detect neuralactivation in dynamic brain imaging data. These analyses require a preliminarily assumed
reference function, which reflects temporal changes in neural activation. In other words,
only those neural activations whose temporal patterns resemble the reference function canbe detected. In cases that reference functions are hardly defined, these analyses are notapplicable. In our previous study, we proposed a method of spatio-temporal filtering toovercome these disadvantages. This method enables us to detect the time and region when
and where dynamical state transition according to neural activation arises in repeatedly
recorded data (multiple trial data). However this method cannot be directly applied tosingle-trial data, such as recording of spontaneous brain activity. In the present study,
we have modified the spatio-temporal filtering method using a sliding time window, andshown its capability to detect neural activation in single-trial data.

Key words: Spatio-temporal ﬁltering, innovation approach, brain functional imaging, optical imaging.