社団法人

電子情報通信学会

               

 脳内電流源の変分ベイズ推定琢Ý ÝÝ

吉岡信学技報     佐藤雅昭ÝÝ ÝÝÝ梶原 茂樹ÝÝÝÝ外山 敬介ÝÝÝÝ

〒  奈良県生駒市高山町  〒  京都府相楽郡精華町光台 
ÝÝÝ 科学技術振興事業団
ÝÝÝÝ 株 島津製作所 基盤技術研究所 〒  京都府相楽郡精華町光台 
 Ý ! "#$ ÝÝ !% "#Ý 奈良先端科学技術大学院大学 情報科学研究科
ÝÝ ＡＴＲ人間情報科学研究所

あらまし

我々は、

 観測データから脳内電流源の分布を推定するための変分ベイズ法に基づく手法を提案している。本手法の事前知識として、電流源はマクロ的には脳内の複数の場所に局在しているという局在条件とミクロ的には電流分布は連続分布をしているという連続条件を同時に考慮に入れ、さらに出来る階層事前分布を導入している。本研究では、

  による情報も考慮に入れる事が 画像データから作成した皮質構造モデルを用いたシミュレーションを行い、その有効性を確かめたので、これを報告する。

 脳内電流源推定 変分ベイズ法 階層事前分布キーワード  Ý ÝÝ                  ÝÝ ÝÝÝ   ÝÝÝÝ    ÝÝÝÝÝ &'( !)%% %* +*%% !(!($ , +( %*ÝÝÝÝ(!)%%-!(!( ' (!)%%-   $+%$ , $ .#ÝÝ / +*%% !(!( 01%($ 2+(% /'$ ($ %$ 3%% $ .#ÝÝÝ$ .# !(!( ' (!)%%- %#%%
((!) 01%%$ )'4 %#%%  /'$ ($ %$ 3%%$ .# Ý ! "#$ ÝÝ !% "#                                                                                                                             !         "      #        $    
何らかの事前知識無しにはこれを解決することは不可能である．

 はじめに

脳磁計 は，脳内の神経細胞活動により発生する微弱な磁場を超伝導磁場センサにより観測し，脳内神経活動を測定する装置であり



や

 と比較して時間解像度が高いため，脳活動の時間変化を詳細に調べるのに適している．
しかし，
布を推定するこの問題点を解決するために，従来は電流源が点状の電流双極子で近似できることを仮定したり，あるいは生理学的な知見から電流源が存在する場所を限定して線形逆フィルタによる推定などを行っていた．しかしこれらの方法は電流源が複数ある複雑な場合に推定の信頼性が悪くなることが指摘され の観測磁場データから脳内電流源の３次元分 ている．
 逆推定問題は原理的な不定性を抱えており， 我々はこのような問題点を解決するため，脳内電流源に対す         

る事前知識として，電流源はマクロ的には脳内の複数の場所に局在しているという局在条件と，ミクロ的には電流分布は連続分布を導入し，変分ベイズ法により電流分布推定を行う方法を提案している．

 と空間分解能に優れたを組み合わせることにより，推定の精度を高める方法が 近年と  で得られる研究されている ．しかし，

  

一方，時間分解能に優れた脳活動は必ずしもいつでも対応するものではないことが最近分

 とかってきている．線形逆フィルタでを組み合わの観測にかからない神経活動をうまく推定せる方法は，することが困難である．逆に活動が有る場所では，たとえ神経電流活動がないときでも，誤って活動電流を推定してしまう可能性がある。

    分布をしているという連続条件を同時に考慮に入れた階層事前                   ここで   はガンマ分布を表す．       は 連続条件を表すために導入された内部変数である．   は事前分布の信頼度を表し，電流分布のスパース（局在）性を制  は事前分布における電   はその期待値である．

御するメタパラメータである． 流の分散 共分散 に逆比例し， り時間分解能も異なるという点を考慮した階層事前分布を提案

  を小さくすることは，双極子  には電流が流れやすい，という事前情報を与えることを意味する．ただし     の値しているを直接指定するのではなく，階層事前分布のメタパラメータ

これに対し我々は，による観測結果を事前知識として推定に利用する際に，と が異なる量を観測してお．

を用いて得られた脳構造モデルに基づいてシ

本稿では，ミュレーション用データを作成し，提案手法の有効性を検証する．特に事前情報として，局在条件，連続条件，情報の 種類の情報を用いた場合の有効性をシミュレーションによっ  の値を指定することで，報として階層事前分布に入れることができる．各双極子ごとに  を与えることによって，電流強度に関する事前情報を与えることができる．例えば，  データを測定した異なる て調べる。また，提案手法の実行時間，および収束性についてときと同じ実験条件において，調べた結果を報告する．

定されている場合，      ，脳内曲面上に仮定した電流双極子に流れる電流分布を          と表す．提案手法では，電流分布  が与えられたときに磁場が観測される確率     が，次のような¼点で測定した磁場をここで 

 はノイズ分散の逆数である． はリードフィールド行列と呼ばれ，双極子の位置，各双極子の電流ベクトルの方向，および磁場センサの位置と方向から計算される

．磁場を の事後分布は次式で与えられる．
              観測したときの電流分布

 

ここで

   は周辺尤度である．   は電流分布に対する我々の事前知識を表す事前分布である．提案手法では，次のような脳内電流源に対する局在条件と連続条件を組み合わせた階層事前分布を導入する．

                      を大きくすることができる．

 は空間フィルタを表す行列であり，次のように与えられる．
      

 は双極子  と  の，脳内曲面（皮質）に沿った距離である．

 をこのように与えることで，事前分布  は，脳内曲面上で近い距離にある双極子同士の電流強度は近くなるという連続条件を表す．連続性はフィルタ半径Ï によって制御され，このパラメータは生理学的な知見に基づいて決められる．
を恒等的に とした場合には，局在条件のみを入れることになる．

    ，  ，階層事前分布を決めるメタパラメータ（  ，  ，  ，   ）と観測磁場が与えられると，電流分布

を含む未知パラメータの事後分布は変分ベイズ法により推定される．変分ベイズ法は自由エネルギーを最大化する逐次アルゴリズムとして実現される．その詳細については文献

 を参照されたい．

 皮質構造データを用いたシミュレーション                 

活動強度の強い格子点では によって，解に対する事前情報を曖昧な情報として与えること

ここで正規分布に従うと仮定する．

        

により脳の活動強度が測を小さく，逆に活動強度が弱ければ  確率モデルと階層事前分布
脳表面近くの情報を曖昧さを持つ情シミュレーションのために，画像データから得られた，

，，）周辺の皮質構造を用いた．，，
周辺の  の格子点（  ）から構成される皮質構造モデルを推定用のモデルとして用いた（図 ）．各格子点から最も近い格子点までの距離は， !! から数 !! 程度である（図視覚野（

Cortex Radius (m)0.1

Cortex theta−phi−0.50.09

Elevation angle

V2, V50.08

V5
Radius

0.070.06

V10.05V5−1V2−1.50.04

V10.03

V20.020100200300400500600700Index−200.20.40.6Azimuth angle

0.81図  シミュレーション皮質モデルの構造情報．左は各格子点のインデックスと深さのプロット．右は方位角と仰角のプロット．
V1

V1=1.0, V2=1.0, V5=0.5

V1=1.0, V2=0.0, V5=0.5

推定用の皮質モデル．黒い点は  ， および  近傍の電流源

図に対応する，推定用の双極子はメッシュの格子点に配置される．
9080
V1=0.0, V2=1.0, V5=0.570

Frequency

605040302010001234Distance to nearest neighbor

56−3x 10図  各格子点から最も近い格子点までの距離のヒストグラム．

図  観測ノイズがない場合に．局在条件を用いて電流源を推定した
結果．

    ，        とし，学習ステップ数を）．これらの格子点上に電流双極子が配置されていると仮定して，サルバスの式を用いて  リードフィールド行列を作成し としたときの，三種類の電流強度分布に対応する推定結果を図  に示す．図の左上，右上，左下は真の電流分布としてそた．皮質の法線方向に電流が流れると仮定し，各双極子の電流れぞれ方向をの位置について正しく推定できていることが分かる．ただし，画像データに基づいて得られた皮質の法線ベクトルで与えた．また チャンネルの磁場センサ（  ）が頭蓋骨の周囲に配置されているとして，リードフィールド行列を計算した．推定モデルの各格子点の深さ，および方位角と仰 ， ， を与えた場合に対応する．図より，電流源推定された電流源は一点に集中しており，その広がりまでは正しく推定できていない．
次に，局在条件に加えて連続条件を導入した場合の結果を示

 に示す． と  は方位角と仰角は近いが，（図  右）， す．    と    は同じ値を用いた．フィルタ半径 を の方がかなり深い位置にあることが分かる（図  左）．この !!，および !! としたときの結果をそれぞれ図 ， に示ため， 電流と  電流を分離して推定することは困難であす．いずれのフィルタ半径を用いても電流源の位置を正しく推る．特に  電流を正しく推定することは難しい．そこで，こ定できるが，電流源の広がりはフィルタ半径に依存して若干のの点を検討するために  種類の電流分布を用意した．すなわち， 違いが見られる．図  は，電流強度分布  に対する推定電流，， に対応する格子点に対して     ， 強度のヒストグラムであり，左上，右上，右下はそれぞれ局在    ，     という三種類の電流強 条件のみ，  !!，  !! とした場合に対応する．局度を与えて，シミュレーション用  データを作成した．電在条件のみを用いた場合，電流強度の高い数点を除いたすべて流源として用いた格子点は，推定用モデルよりも高い解像度をの点で電流が  になっていることが分かる．連続条件を加える
角を図用いた．このためノイズが全くない場合でもシミュレーション
用 データを推定用モデルで完全に再現することはできな¿º½

いことに注意する．
観測ノイズがない場合¿º¾と，推定電流強度の分布が広がることが分かる．
観測ノイズがある場合

次に，観測ノイズとして"# 比  ，および   のガウスノイズを加えた場合のシミュレーションを行なった．推定には局在最初に観測ノイズがない場合の推定結果を示す．まず，局条件と連続条件を組み合わせた階層事前分布を用いた．階層事在階 層事前分 布を用い て推定を 行った．メタパ ラメータ を前分布のメタパラメータ，および学習回数はノイズがない場合

V1=1.0, V2=1.0, V5=0.5

V1=1.0, V2=0.0, V5=0.5

V1=1.0, V2=1.0, V5=0.5

V1=0.0, V2=1.0, V5=0.5

V1=1.0, V2=0.0, V5=0.5

V1=0.0, V2=1.0, V5=0.5

図  局在条件と連続条件を用いて電流源を推定した結果（  ）．

比  のノイズを加えたシミュレーション用観測データに対図して，局在条件と連続条件を用いて電流源を推定した結果．
V1=1.0, V2=1.0, V5=0.5

V1=1.0, V2=0.0, V5=0.5
V1=1.0, V2=1.0, V5=0.5

V1=1.0, V2=0.0, V5=0.5

V1=0.0, V2=1.0, V5=0.5
V1=0.0, V2=1.0, V5=0.5

図  局在条件と連続条件を用いて電流源を推定した結果（  ）．
比  のノイズを加えたシミュレーション用観測データに対400して，スパース事前分布と連続条件を用いて電流源を推定した

600300結果．

Frequency

Frequency

図800

400
2000−10200

電流強度に関する事前知識を導入して推定を行った．電流源付

1000102030

Estimated current

0−2024

Estimated current

200

Frequency

150
10050

!! 以内に含まれる格子点（図 ）に対して階層事     とし，それ以外の格子点では       とした．その結果， にある電流源の位置と広がりを正しく推定することができた（図 ）．
このシミュレーション結果から次のようなことが分かる．

近から前分布のメタパラメータを 

真の活動領域よりも広い領域に活動があるという情報を入れて0−2024

Estimated current

図推定電流強度のヒストグラム．

も，

 データから電流の広がりが正しく推定できた．

電流活動がない場所に活動があるという誤った情報を与えても（  と  ），電流が無い場所に誤って電流を推定することはな

!! とした．
図 ， にその推定結果を示す．ノイズが小さい場合（"# ）のシミュレーションと同じにし，フィルタ半径はは電流源が正しく推定されている．しかし，ノイズが大きい場合（ある"# ）は，真の電流分布  ， に対して，深い位置に 上の電流源が正しく推定されていない．そこで， 

観測に対応するデータが得られたという状況を想定して，

¿º¿

かった．
アルゴリズムの実行時間と収束性

推定アルゴリズムの計算量は双極子数の三乗のオーダーであり，双極子数が大きい場合には多くの計算量を必要とする．推 ステップの計算を  $% の &'()! で測定した．双極子数が  と  のモデルに対する平均計算時間は，それぞれ   秒，  秒であった．したがって，
定アルゴリズムのsparse, SN=0−32x 10

R=5mm, SN=01

1.5
10.50.500A
100

200
300
400R=9mm, SN=0500

B000.110.05

0.5100200
300
400R=5mm, SN=0.1

C00100200
300
400R=5mm, SN=0.2

D500

001

10.50.5100200300
400R=5mm, fMRI, SN=0.2E00

図  事前電流強度を大きくした（ を小さくした）双極子の位置．

100200300

400

V1=1.0, V2=0.0, V5=0.5500

F50000100200300

400500

図  推定アルゴリズムの収束性を示すプロット．横軸は学習回数，縦軸は  回学習後の推定結果に対する正規化二乗誤差を表す． から  は，図  から

V1=1.0, V2=1.0, V5=0.5500に結果を示したシミュレーションにそれぞれ対応する．
sparse, SN=0R=5mm, SN=02600

2600

2400

2400

2200

2200

2000

2000

1800
0V1=0.0, V2=1.0, V5=0.5

1800

A
100

200
300
400
R=9mm, SN=0500

2600

2600

2400

2400

2200

2200

2000

100200300
400
R=5mm, SN=0.1500

2000

1800
0B0C
100

200
300
400
R=5mm, SN=0.2500

2600

1800

D

0100
200
300
400
R=5mm, fMRI, SN=0.2500

2600
E

2400

2400

2200

2200

2000

2000

1800

比  のノイズを加えたシミュレーション用観測データに図

対して，局在条件，連続条件，電流強度情報を用いて電流源を

01800100200300400500F

0100200300400500

図  シミュレーションにおける自由エネルギーの変化．

推定した結果．

 回の学習には約  時間かかる．しかも，
 時刻における  観測データに対する推定にかかる回までには学習がほぼ収束している．ノイズがある場合の変分点のモデルに対するベイズ推定では，ノイズ分散以上に誤差を小さくしようとしなこれはいので，ノイズが無い場合よりも収束が早くなると考えられる．

時間であり，時系列データに対して推定を行う場合には，解析する時間ステップ数に比例した時間が必要になる．そこで，妥当な推定結果が得られるまでに必要な学習回数について調べた．

 ま と め
本稿では，我々が提案する変分ベイズ法による

 電源推 回の学習で妥当な結果が得られているので，推定電流強度が  回学定の有効性を検証するために行ったシミュレーション結果につ習後の推定電流強度に収束するまでの挙動を調べる．各学習回

局在条件と連続条件を組み合わせた階層事前分布用いること上で示したシミュレーション結果では，いずれも数における推定電流強度と 回学習後の推定電流強度との正 に示す．また，学習回数に対する自由エ に示す．ノイズが無い場合に局在条件のみを用いて推定した場合，約  回で収束している．一方，フィルタ半径 !! でノイズが無い場合が最も収束が遅い．局いて説明した．
によって，から得られた皮質構造モデルを用いたシミュ

規化二乗誤差を，図レーションにおいて，電流源の位置とその広がりを正しく推定

ネルギーの収束性を図することができた．また，在条件のみを用いた場合，電流強度が強い数個の双極子のみを用いた推定になるので収束が早くなり，一方フィルタ半径を大きくすると解の不定性のために（データを完全に復元する解はないことに注意する）性能的に似たような解が多くあり，緩やかに収束していると解釈できる．ノイズがある場合には，

による情報を模した電流強度に対する事前情報を曖昧な形で階層事前分布に与えることにより，観測ノイズが大きい場合でも，皮質の深い位置にある電流

源の位置と広がりを正しく推定でき 特に以下の事が示された．

 真の活動領域よりも広い領域に活動があるという情報を データから電流の広がりが正しく推定できた．
 電流活動がない場所に活動があるという誤った情報を入れても，与えても，電流が無い場所に誤って電流を推定することはな



かった．
この事は


と

 で得られる脳活動が一致していな

いときでも 本手法を用いて安全に推定が出来ることを示唆し



ている。
実際の観測データを用いて本手法の有効性を検証することは今後の課題である．
謝辞

本研究に有益な議論をして下さった*
人間情報科学研究所、川人光男室長と銅谷賢治博士に感謝します。本研究は通信・放送機構の研究委託により実施したものである。
文          

献佐 藤雅昭，変 分ベイズ法 による  電源推定 ，信学技 報，   
佐藤雅昭，吉岡琢，梶原茂樹，外山敬介， 電源推定とベイズ推定事前分布，信学技報（本予稿），
  !
" #$%&''$( )'&&' $* '+,+'& -% .
&$ /'0.& -*&. 01 /'-$. 21&       !))  "3+',3.$* 4  $* &5$-67 8 9 !
" 4&')&,/&-$ ,-($/& $&+'-,&$,&      : ' ; !
" <$ 0'$= ($&, '&)-$& - &$-'1+.        :  >53'  <3$ ? -*  8$(6  8-1 >
!" @&$&' A.&'($&-&$,&)3.-(')31 -% /+. ,-',. ,/&    +0&*&. ;  ! "         )'$(&'&'.(-  ! " <$.$& -*&. &.&,-$ 2&* -$ 3& '-$. 21&     '/ 9 !
" 2, 3&,. $* &.&,'-($&,,-$,&) -% 3& 0-($&, $/&'& )'-0.& !
  :.$&$  ' ; #.-$& ; 9 >$++. 9
$* B-+$ <  !
" ($&-&$,&)3.-(')31
83&-'1 $'+&$-$ $* )).,-$ - $-$$//&+*& -% 3& 6-'7$( 3+$ 0'$ " # !   : 8-1 > ?-376 > ?-3* ? >-$*- 
8- ? 87$3 ? 5 ? $* ?-3C6 !
"  $&6 &3-* %-' ($&-&$,&)3.-(')31 !" D83'&& *&$-$. ($&-&&'). A.&' 1&     : @$( 9E @..-$  9 $* >+%$ B ! "($&, -+',& (& *&&'$&* 01  .&*A&.* $.1D
83& +$F+& $+$-' .&F+'& &-$ $$$%   ! $   : 2$*&$ 4 !" 83& &)-'. '&-.+-$ -% %+$,
-$. ;# $ &   " &* 01 --$&$ 8@$* 2$*&$ 4 )'$(&' : 
2&../&+ 9 @ >&$$&*1 G  ,>$'1 ;  2+,3
0$*&' 2 ; @&7-H ;  -3&$   &/& 9 2'*1 8 9 $* ;-&$ 2 ; !
" +$,-$. ))$(-% 3& 3+$ /+. ,-'&I 01 ($&, '&-$$,& ($(   : 
3+',3.$* 4  $* &5$-67 8 9 !
" 4&')&,/&-$ ,-($/& $&+'-,&$,&      : B-(-3& > 4+. 9 +(3  8'$3 8 $* <&.&'$$  ! " &+'-)31-.-(,. $/&(-$ -% 3&0 -% 3& %;# ($.    :  <(6  B&& 8 >1 ; 8$7 G@ ! " 2'$($&, '&-$$,& ($( 63 ,-$'*&)&$*&$ -I1(&$-$   !     : 3.%-'  4 )-$   G.&   2&../&+ 9@ B+  > >-'/&$-5  '$&$ 9 +-.$&$ 8--&.. ; 2  '-$&$  9 $* #.-$& ; 

9 !
" )-&)-'. ,/1 -%  ,-',. $&6-'7 %-')'-,&$( /+. --$ '&/&.&* 01  $* %;#     :
G.&   $* &'&$-  # ! " #)'-/&* .-,.C-$ -% ,-',. ,/1 01 ,-0$$(  $*  63
;# ,-',. +'%,& '&,-$'+,-$D  B$&' ))'-,3      : G.&   B+  > ,3. 2 ; 2+,3$&' ; B
2&../&+ 9 @ B&6$& 9 G $* .('&$  !"G1$, ,. )'&', ))$(D -0$$( %;#$*  %-' 3(3'&-.+-$ ($( -% ,-',. ,/1   : 
+57  176 8 &.&$  >J,3& 8; $*1+,3  !" $ %;#,-$'$&* -+',&
$.1 63 )'-,&*+'& %-' */*$( $* ('-+)$(    :
&$C&  9 $(+$  ; 2+',3&' @ $',3 ,3-.C  +$& 8  -  ,3&'(  9-3$$&
 +$*&3(&$  CC$(   $* ..1'*  ! " -0$&* ). $* &)-'. ($( -% 0'$ ,/1 *+'$( /+. &.&,/& &$-$ $ 3+$   :
B+  > 2&../&+ 9 @ $* G.&   !
")-&)-'. ($( -% 3+$ 0'$ ,/1 +$(%+$,-$. ;# ,-$'$&* ($&-&$,&)3.-(')31 *D
-$& '.- +.-$   !   
: 