TOP > 国内特許検索 > 脳内電流源推定方法、脳内電流源推定プログラム、脳内電流源推定プログラムを記録した記録媒体および脳内電流源推定装置 > 明細書

明細書 :脳内電流源推定方法、脳内電流源推定プログラム、脳内電流源推定プログラムを記録した記録媒体および脳内電流源推定装置

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第3730646号 (P3730646)
登録日 平成17年10月14日(2005.10.14)
発行日 平成18年1月5日(2006.1.5)
発明の名称または考案の名称 脳内電流源推定方法、脳内電流源推定プログラム、脳内電流源推定プログラムを記録した記録媒体および脳内電流源推定装置
国際特許分類 A61B   5/05        (2006.01)
A61B   5/0476      (2006.01)
FI A61B 5/05 A
A61B 5/04 322
請求項の数または発明の数 20
全頁数 38
出願番号 特願2003-557403 (P2003-557403)
出願日 平成14年12月27日(2002.12.27)
国際出願番号 PCT/JP2002/013739
国際公開番号 WO2003/057035
国際公開日 平成15年7月17日(2003.7.17)
優先権出願番号 2001400519
優先日 平成13年12月28日(2001.12.28)
優先権主張国 日本国(JP)
審査請求日 平成15年9月29日(2003.9.29)
特許権者または実用新案権者 【識別番号】503360115
【氏名又は名称】独立行政法人科学技術振興機構
【識別番号】393031586
【氏名又は名称】株式会社国際電気通信基礎技術研究所
発明者または考案者 【氏名】佐藤 雅昭
個別代理人の代理人 【識別番号】100064746、【弁理士】、【氏名又は名称】深見 久郎
【識別番号】100085132、【弁理士】、【氏名又は名称】森田 俊雄
【識別番号】100096781、【弁理士】、【氏名又は名称】堀井 豊
【識別番号】100064746、【弁理士】、【氏名又は名称】深見 久郎
【識別番号】100085132、【弁理士】、【氏名又は名称】森田 俊雄
【識別番号】100083703、【弁理士】、【氏名又は名称】仲村 義平
【識別番号】100096781、【弁理士】、【氏名又は名称】堀井 豊
【識別番号】100098316、【弁理士】、【氏名又は名称】野田 久登
【識別番号】100109162、【弁理士】、【氏名又は名称】酒井 將行
審査官 【審査官】門田 宏
参考文献・文献 特表2001-514031(JP,A)
特開平7-327946(JP,A)
特開平5-197767(JP,A)
特開平6-217954(JP,A)
特開平10-211181(JP,A)
特開平10-124851(JP,A)
調査した分野 A61B 5/05
A61B 5/04
JICSTファイル(JOIS)
特許請求の範囲 【請求項1】
頭蓋外部において観測された電磁場に基づいて、脳内に存在する前記電磁場の源となる電流源の位置を推定するための方法であって、
前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定するステップと、
各前記仮想曲面上において、それぞれ前記観測された電磁場を復元するための電流分布を推定するステップと、
前記複数の仮想曲面上において推定された電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定するステップとを備える、脳内電流源推定方法。
【請求項2】
前記電流分布を推定するステップは、
事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定するステップを含む、請求項1記載の脳内電流源推定方法。
【請求項3】
前記電流分布を推定するステップは、
前記複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、前記脳表面に最も近く、かつ事後確率が極大になる第1の仮想曲面を特定するステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記第1の仮想曲面において、前記電流分布の極大点に対応する局在電流分布を特定するステップと、
前記局在電流分布を各々が含む複数の局所面を切り出すステップと、
前記複数の局所面のうち、特定対象となる局所面以外を固定し、前記特定対象となる局所面を深さ方向に動かして、前記事後確率が極大となる位置を、前記脳表面に近い側から順次特定するステップとを含む、請求項2記載の脳内電流推定方法。
【請求項4】
前記電流分布を推定するステップにおいては、第1の空間分解能で、前記電流分布を推定し、
前記複数の仮想曲面の深さ方向の分解能を高め、かつ、前記第1の空間分解能よりも高い第2の空間分解能で、前記電流分布の再推定を行うステップをさらに備える、請求項3記載の脳内電流源推定方法。
【請求項5】
前記電流分布を推定するステップは、前記電流分布をベイズ推定する際に、前記電流源の推定のための前記電磁場の観測とは独立した他の観測方法による観測データを用いて、前記ベイズ推定における階層事前分布を設定するステップを含む、請求項1記載の脳内電流源推定方法。
【請求項6】
頭蓋外部において観測された電磁場に基づいて、脳内に存在する前記電磁場の源となる電流源の位置を推定するコンピュータのためのプログラムであって、
前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定するステップと、
各前記仮想曲面上において、それぞれ前記観測された電磁場を復元するための電流分布を推定するステップと、
前記複数の仮想曲面上において推定された電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定するステップとを実行させるためのプログラム。
【請求項7】
前記電流分布を推定するステップは、
事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定するステップを含む、請求項6記載のプログラム。
【請求項8】
前記電流分布を推定するステップは、
前記複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、前記脳表面に最も近く、かつ事後確率が極大になる第1の仮想曲面を特定するステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記第1の仮想曲面において、前記電流分布の極大点に対応する局在電流分布を特定するステップと、
前記局在電流分布を各々が含む複数の局所面を切り出すステップと、
前記複数の局所面のうち、特定対象となる局所面以外を固定し、前記特定対象となる局所面を深さ方向に動かして、前記事後確率が極大となる位置を、前記脳表面に近い側から順次特定するステップとを含む、請求項7記載のプログラム。
【請求項9】
前記電流分布を推定するステップにおいては、第1の空間分解能で、前記電流分布を推定し、
前記複数の仮想曲面の深さ方向の分解能を高め、かつ、前記第1の空間分解能よりも高い第2の空間分解能で、前記電流分布の再推定を行うステップをさらに備える、請求項8記載のプログラム。
【請求項10】
前記電流分布を推定するステップは、前記電流分布をベイズ推定する際に、前記電流源の推定のための前記電磁場の観測とは独立した他の観測方法による観測データを用いて、前記ベイズ推定における階層事前分布を設定するステップを含む、請求項6記載のプログラム。
【請求項11】
頭蓋外部において観測された電磁場に基づいて、脳内に存在する前記電磁場の源となる電流源の位置を推定するための脳内電流源推定装置であって、
前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定する仮想曲面設定手段と、
各前記仮想曲面上において、それぞれ前記観測された電磁場を復元するための電流分布を推定する電流分布推定手段と、
前記複数の仮想曲面上において推定された電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定する電流源特定手段とを備える、脳内電流源推定装置。
【請求項12】
前記電流分布推定手段は、
事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求める事後確率導出手段を含み、
前記電流源特定手段は、
前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定する仮想曲面特定手段を含む、請求項11記載の脳内電流源推定装置。
【請求項13】
前記電流分布推定手段は、
前記複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、前記脳表面に最も近く、かつ事後確率が極大になる第1の仮想曲面を特定する最浅仮想曲面特定手段を含み、
前記電流源特定手段は、
前記第1の仮想曲面において、前記電流分布の極大点に対応する局在電流分布を特定する局在電流分布特定手段と、
前記局在電流分布を各々が含む複数の局所面を切り出す局所面抽出手段と、
前記複数の局所面のうち、特定対象となる局所面以外を固定し、前記特定対象となる局所面を深さ方向に動かして、前記事後確率が極大となる位置を、前記脳表面に近い側から順次特定する局所面位置特定手段とを含む、請求項12記載の脳内電流源推定装置。
【請求項14】
前記電流分布推定手段は、初期的に第1の空間分解能で、前記電流分布を推定した後、前記複数の仮想曲面の深さ方向の分解能を高め、かつ、前記第1の空間分解能よりも高い第2の空間分解能で、前記電流分布の再推定を行う、請求項13記載の脳内電流源推定装置。
【請求項15】
前記電流分布推定手段は、前記電流分布をベイズ推定する際に、前記電流源の推定のための前記電磁場の観測とは独立した他の観測方法による観測データを用いて、前記ベイズ推定における階層事前分布を設定する手段を含む、請求項11記載の脳内電流源推定装置。
【請求項16】
頭蓋外部において観測された電磁場に基づいて、脳内に存在する前記電磁場の源となる電流源の位置を推定するコンピュータのためのプログラムを記録した記録媒体であって、前記プログラムは、
前記脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各前記仮想曲面上に格子点を設定するステップと、
各前記仮想曲面上において、それぞれ前記観測された電磁場を復元するための電流分布を推定するステップと、
前記複数の仮想曲面上において推定された電流分布の広がりと、前記電流分布に基づいて復元される電磁場と前記観測された電磁場との誤差とに基づいて、前記複数の仮想曲面のうちから前記広がりと誤差とが極小となる仮想曲面を、真の前記電流源の存在する曲面として特定するステップとを備える記録媒体。
【請求項17】
前記電流分布を推定するステップは、
事前分布と前記電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記仮想曲面のうち、対応する前記事後確率が最大となる仮想曲面を特定するステップを含む、請求項16記載の記録媒体。
【請求項18】
前記電流分布を推定するステップは、
前記複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、前記脳表面に最も近く、かつ事後確率が極大になる第1の仮想曲面を特定するステップを含み、
前記真の電流源の存在する曲面として特定するステップは、
前記第1の仮想曲面において、前記電流分布の極大点に対応する局在電流分布を特定するステップと、
前記局在電流分布を各々が含む複数の局所面を切り出すステップと、
前記複数の局所面のうち、特定対象となる局所面以外を固定し、前記特定対象となる局所面を深さ方向に動かして、前記事後確率が極大となる位置を、前記脳表面に近い側から順次特定するステップとを含む、請求項17記載の記録媒体。
【請求項19】
前記電流分布を推定するステップにおいては、第1の空間分解能で、前記電流分布を推定し、
前記複数の仮想曲面の深さ方向の分解能を高め、かつ、前記第1の空間分解能よりも高い第2の空間分解能で、前記電流分布の再推定を行うステップをさらに備える、請求項18記載の記録媒体。
【請求項20】
前記電流分布を推定するステップは、前記電流分布をベイズ推定する際に、前記電流源の推定のための前記電磁場の観測とは独立した他の観測方法による観測データを用いて、前記ベイズ推定における階層事前分布を設定するステップを含む、請求項16記載の記録媒体。
発明の詳細な説明 技術分野
本発明は、脳内活動に対応して、頭蓋外部に磁場または電場を生成する脳内電流源の位置や分布を推定する方法、脳内電流源推定プログラム、当該プログラムを記録した記録媒体および脳内電流源推定装置に関する。
背景技術
近年の生体計測技術の進歩はめざましく、従来測定が困難で誤差も大きかった脳から発生する微弱電場(脳波)や微弱磁場(脳磁波)の計測精度は、年々向上している。
すなわち、外部からの刺激を受けて、脳内の神経細胞は電流を発生する。この電流により、上述したような微弱な電場や微弱な磁場が発生することになる。このうち「脳波」とは、この神経細胞の電流により、脳から発生する電場のことである。「脳波計(Electroencephalogram:EEG)」は、このような脳波を測定する手法である。
一方、「脳磁波」とは、この神経細胞の電流により、脳から発生する磁場のことである。「脳磁計(Magnetoencephalography:MEG)」は、このような脳磁波を測定する手法である。脳磁波計測の最大の利点は、磁場が容積導体による影響をほとんど受けないので、頭蓋外からの磁気計測により、脳内電流位置を3次元的に比較的正確に推定できることが期待されることである。
このような脳磁波の解析においては、発生した磁場を脳の外部から測定することで非侵襲的に脳の活動部位を推定する。
しかしながら、この磁場は非常に微弱であるために、地磁気等の外部の磁場の影響を大きく受けていまう。そこで、外部の磁場を遮蔽するシールドを作り、その中で超伝導を用いた測定機器である超伝導量子干渉素子(SQUID:Superconducting quantum interference device)により微弱な磁場を測定することが行われる。
ただし、このような「脳内電流源位置推定」のアルゴリズムの研究においても、初期のモデルから様々な変法が使用されているにも関わらず、決定的な手法は今だ存在しないのが現状である。
たとえば、このような「脳内電流源位置推定」のアルゴリズムの1つである「双極子推定法」が、文献1:Mosher.J.C.,Lewis P.S.and Leahy R.M.IEEE Trans.Biomed.Engng.〈1992〉vol.39,pp.541-557に開示されている。しかしながら、この「双極子推定法」は、脳内の電流源が1つか数個の電流双極子で表せると仮定して観測磁場から双極子の位置を推定する方法であり、双極子の数を幾つにすれば良いかを決めるのが難しいという欠点を持っている。
一方、他のアルゴリズムとしては、「空間フィルター法」が、文献2:Toyama.K,Yoshizawa K,Yoshida Y,Kondo Y,Tomita S,Takanashi Y,Ejima Y,and Yoshikawa S.Neuroscience.1999;91(2):pp.405-415に開示されている。「空間フィルター法」は、生理学的な知見から、脳内電流源の場所を限定し、双極子の分布を推定する方法である。しかしながら、この方法の欠点は電流源の深さを正確に推定できない事である。
発明の開示
本発明の目的は、観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能な脳内電流源推定方法を提供することである。
この発明の他の目的は、複数の脳内電流源が存在する場合にも、観測データから複数の脳内電流源の位置を推定することが可能な脳内電流源推定方法を提供することである。
この発明のさらに他の目的は、脳内の電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データが同時に利用できる場合に、複数の観測方法によるデータを組み合わせて、より精度の高い推定を行う脳内電流源推定方法を提供することである。
この発明のさらに他の目的は、観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能な脳内電流源推定プログラム、当該プログラムを記録した記録媒体を提供することである。
この発明のさらに他の目的は、複数の脳内電流源が存在する場合にも、観測データから複数の脳内電流源の位置を推定することが可能な脳内電流源推定プログラム、当該プログラムを記録した記録媒体を提供することである。
この発明のさらに他の目的は、脳内の電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データが同時に利用できる場合に、複数の観測方法によるデータを組み合わせて、より精度の高い推定を行う脳内電流源推定プログラム、当該プログラムを記録した記録媒体を提供することである。
この発明のさらに他の目的は、観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能な脳内電流源推定装置を提供することである。
この発明のさらに他の目的は、複数の脳内電流源が存在する場合にも、観測データから複数の脳内電流源の位置を推定することが可能な脳内電流源推定装置を提供することである。
この発明のさらに他の目的は、脳内の電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データが同時に利用できる場合に、複数の観測方法によるデータを組み合わせて、より精度の高い推定を行う脳内電流源推定装置を提供することである。
係る目的を達成するために本願発明に係る脳内電流源推定方法は、頭蓋外部において観測された電磁場に基づいて、脳内に存在する電磁場の源となる電流源の位置を推定するための方法であって、脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各仮想曲面上に格子点を設定するステップと、各仮想曲面上において、それぞれ観測された電磁場を復元するための電流分布を推定するステップと、複数の仮想曲面上において推定された電流分布の広がりと、電流分布に基づいて復元される電磁場と観測された電磁場との誤差とに基づいて、複数の仮想曲面のうちから広がりと誤差とが極小となる仮想曲面を、真の電流源の存在する曲面として特定するステップとを備える。
好ましくは、電流分布を推定するステップは、事前分布と電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、真の電流源の存在する曲面として特定するステップは、仮想曲面のうち、対応する事後確率が最大となる仮想曲面を特定するステップを含む。
好ましくは、電流分布を推定するステップは、複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、脳表面に最も近く、かつ事後確率が極大になる第1の仮想曲面を特定するステップを含み、真の電流源の存在する曲面として特定するステップは、第1の仮想曲面において、電流分布の極大点に対応する局在電流分布を特定するステップと、局在電流分布を各々が含む複数の局所面を切り出すステップと、複数の局所面のうち、特定対象となる局所面以外を固定し、特定対象となる局所面を深さ方向に動かして、事後確率が極大となる位置を、脳表面に近い側から順次特定するステップとを含む。
好ましくは、電流分布を推定するステップにおいては、第1の空間分解能で、電流分布を推定し、複数の仮想曲面の深さ方向の分解能を高め、かつ、第1の空間分解能よりも高い第2の空間分解能で、電流分布の再推定を行うステップをさらに備える。
好ましくは、電流分布を推定するステップは、電流分布をベイズ推定する際に、電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データを用いて、ベイズ推定における階層事前分布を設定するステップを含む。たとえば、他の観測方法による観測データにより脳内電流源が存在する場所が限られた領域内にあることが分かっている場合には、深さ方向の探索を省略して、限定された領域内における電流分布をベイズ推定により求めることができる。
この発明の他の局面に従うと、頭蓋外部において観測された電磁場に基づいて、脳内に存在する電磁場の源となる電流源の位置を推定するコンピュータのためのプログラムであって、脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想曲面と、各仮想曲面上に格子点を設定するステップと、各仮想曲面上において、それぞれ観測された電磁場を復元するための電流分布を推定するステップと、複数の仮想曲面上において推定された電流分布の広がりと、電流分布に基づいて復元される電磁場と観測された電磁場との誤差とに基づいて、複数の仮想曲面のうちから広がりと誤差とが極小となる仮想曲面を、真の電流源の存在する曲面として特定するステップとを実行させる。
好ましくは、電流分布を推定するステップは、事前分布と電磁場の観測データとから、ベイズ推定法により事後確率を求めるステップを含み、
真の電流源の存在する曲面として特定するステップは、仮想曲面のうち、対応する事後確率が最大となる仮想曲面を特定するステップを含む。
好ましくは、電流分布を推定するステップは、複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、脳表面に最も近く、かつ事後確率が極大になる第1の仮想曲面を特定するステップを含み、真の電流源の存在する曲面として特定するステップは、第1の仮想曲面において、電流分布の極大点に対応する局在電流分布を特定するステップと、局在電流分布を各々が含む複数の局所面を切り出すステップと、複数の局所面のうち、特定対象となる局所面以外を固定し、特定対象となる局所面を深さ方向に動かして、事後確率が極大となる位置を、脳表面に近い側から順次特定するステップとを含む。
好ましくは、電流分布を推定するステップにおいては、第1の空間分解能で、電流分布を推定し、複数の仮想曲面の深さ方向の分解能を高め、かつ、第1の空間分解能よりも高い第2の空間分解能で、電流分布の再推定を行うステップをさらに備える。
好ましくは、電流分布を推定するステップは、電流分布をベイズ推定する際に、電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データを用いて、ベイズ推定における階層事前分布を設定するステップを含む。たとえば、他の観測方法による観測データにより脳内電流源が存在する場所が限られた領域内にあることが分かっている場合には、深さ方向の探索を省略して、限定された領域内における電流分布をベイズ推定により求めることができる。
この発明のさらに他の局面に従うと、頭蓋外部において観測された電磁場に基づいて、脳内に存在する電磁場の源となる電流源の位置を推定するための脳内電流源推定装置であって、脳内に脳表面からの深さが互いに異なり、かつ互いに交わらない形状を有する複数の仮想局面と、各仮想曲面上に格子点を設定する仮想曲面設定手段と、各仮想曲面上において、それぞれ観測された電磁場を復元するための電流分布を推定する電流分布推定手段と、複数の仮想曲面上において推定された電流分布の広がりと、電流分布に基づいて復元される電磁場と観測された電磁場との誤差とに基づいて、複数の仮想曲面のうちから広がりと誤差とが極小となる仮想曲面を、真の電流源の存在する曲面として特定する電流源特定手段とを備える。
好ましくは、電流分布推定手段は、事前分布と電磁場の観測データとから、ベイズ推定法により事後確率を求める事後確率導出手段を含み、電流源特定手段は、仮想局面のうち、対応する事後確率が最大となる仮想曲面を特定する仮想曲面特定手段を含む。
好ましくは、電流分布推定手段は、複数の仮想曲面のうち、脳表面側の仮想曲面から順次深い側の仮想曲面に移行するのに応じて、脳表面に最も近く、かつ事後確率が極大になる第1の仮想曲面を特定する最浅仮想曲面特定手段を含み、電流源特定手段は、第1の仮想曲面において、電流分布の極大点に対応する局在電流分布を特定する局在電流分布特定手段と、局在電流分布を各々が含む複数の局所面を切り出す局所面抽出手段と、複数の局所面のうち、特定対象となる局所面以外を固定し、特定対象となる局所面を深さ方向に動かして、事後確率が極大となる位置を、脳表面に近い側から順次特定する局所面位置特定手段とを含む。
好ましくは、電流分布推定手段は、初期的に第1の空間分解能で、電流分布を推定した後、複数の仮想曲面の深さ方向の分解能を高め、かつ、第1の空間分解能よりも高い第2の空間分解能で、電流分布の再推定を行う。
好ましくは、電流分布推定手段は、電流分布をベイズ推定する際に、電流源の推定のための電磁場の観測とは独立した他の観測方法による観測データを用いて、ベイズ推定における階層事前分布を設定する手段を含む。たとえば、他の観測方法による観測データにより脳内電流源が存在する場所が限られた領域内にあることが分かっている場合には、深さ方向の探索を省略して、限定された領域内における電流分布をベイズ推定により求めることができる。
本発明の脳内電流源推定方法に従えば、MEGやEEGの観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能である。しかも、このような深さ方向の推定は、複数の電流源が存在する場合にも適用可能である。さらに、電流源が電流双極子のように局在している場合にも、広がりを持つ電流源の場合にも適用できるだけでなく、電流源の広がり具合を推定することも可能である。
また、初期分解能での推定の後に、位置の推定分解能を逐次的に上げる処理を追加して行うことが可能である。
【図面の簡単な説明】
図1は、脳磁計システムの構成の一例を示す概念図である。
図2は、電流源が作る電磁場を適当な曲面上で観測している場合の概念図である。
図3は、脳内の電流源と複数の仮想曲面との関係を示す概念図である。
図4は、本発明の脳内電流源推定方法の手順の全体的な流れを示すフローチャートである。
図5は、初期分解能を用いた電流源の初期推定の処理を説明するためのフローチャートである。
図6は、最も脳表面に近い電流源の特定を行なう処理を説明するためのフローチャートである。
図7は、各局所面に対応する電流源の深さの特定処理を説明するための第1のフローチャートである。
図8は、各局所面に対応する電流源の深さの特定処理を説明するための第2のフローチャートである。
図9は、空間分解能を上げて電流源の位置を再推定する処理を説明するための第1のフローチャートである。
図10は、空間分解能を上げて電流源の位置を再推定する処理を説明するための第2のフローチャートである。
図11は、人間の脳を半球状であると仮定し、この半球の表面上で観測された磁場分布を上面から見た上面図である。
図12(a)、図12(b)および図12(c)は、それぞれ、半径R5.0、半径R=6.0および半径R=7.0の場合において、初期分解能を用いた電流源の初期推定結果を示す図である。
図13(d)、図13(e)および図13(f)は、それぞれ、半径R=7.5、半径R=8.0および半径R=9.0の場合であって、脳表面により近い場合の初期分解脳を用いた電流源の初期推定結果を示す図である。
図14は、自由エネルギーを最大とするような電流分布を各々求める際に得られる、各仮想半球面上での自由エネルギーの大きさを示す図である。
図15(a)は半径R=5.0、図15(b)は半径R=6.0、図15(c)は半径R=7.0の場合において、極大点を含む局所面についてさらに電流源の特定処理を行なった場合の処理をそれぞれ示す図である。
図16(d)は半径R=7.5、図16(e)は半径R=8.0、図16(f)は半径R=9.0の場合において、極大点を含む局所面についてさらに電流源の特定処理を行なった場合の処理をそれぞれ示す図である。
図17は、局所面の深さを変えつつ自由エネルギーを求めた結果を示す図である。
図18は、脳表面を半球と仮定した場合に、電流源が2つ脳内に存在する場合の脳表面における磁場分布を示す上面図である。
図19(a)は半径R=5.0、図19(b)は半径R=6.0、図19(c)は半径R=7.0の場合において、初期分解能を用いた電流源の初期推定結果を示す図である。
図20(d)は半径R=7.5、図20(e)は半径R=8.0、図20(f)は半径R=9.0の場合において、初期分解能を用いた電流源の初期推定結果を示す図である。
図21は、各仮想半球面ごとに自由エネルギーが最大となるように電流分布を求めた場合の、自由エネルギーの半径依存性を示す図である。
図22(a)は半径R=5.0、図22(b)は半径R=6.0、図22(c)は半径R=7.0の場合の局所面上の電流分布をそれぞれ示し、最も脳表面に近い電流源を特定するために、事後確率が最大になる1番目の局所面の深さを求める処理を説明する図である。
図23(d)は、半径R=7.5、図23(e)は半径R=8.0、図23(f)は半径R=9.0の場合の局所面上の電流分布をそれぞれ示し、最も脳表面に近い電流源を特定するために、事後確率が最大になる1番目の局所面の深さを求める処理を説明する図である。
図24は、仮想局所面を移動させた場合の自由エネルギーの局所面位置(半径)依存性を示す図である。
図25(a)は半径R=5.0、図25(b)は半径R=5.5、図25(c)は半径R=6.0の場合の局所面上の電流分布を示し、特定された面上に一方の局所面を固定したままで、他方の電流源に対応する局所面を移動させた場合の電流分布を示す図である。
図26(d)は、半径R=6.5、図26(e)は半径R=7.0、図26(f)は半径R=7.5の場合の局所面上の電流分布をそれぞれ示し、特定された面上に一方の局所面を固定したままで、他方の電流源に対応する局所面を移動させた場合の電流分布を示す図である。
図27は、仮想局所面を移動させた場合の自由エネルギーの局所面位置(半径)依存性を示す図である。
図28は、局在条件のみを考慮したシミュレーション結果を示す第1の図である。
図29は、局在条件のみを考慮したシミュレーション結果を示す第2の図である。
図30は、連続条件および局在条件の双方を考慮したシミュレーションを行った結果を示す第1の図である。
図31は、連続条件および局在条件の双方を考慮したシミュレーションを行った結果を示す第2の図である。
図32は、連続条件および局在条件の双方を考慮し、かつ電流源が存在する可能性の高い領域内で電流が推定されやすくなる階層事前分布を用いた場合のシミュレーションを行った結果を示す第1の図である。
図33は、連続条件および局在条件の双方を考慮し、かつ電流源が存在する可能性の高い領域内で電流が推定されやすくなる階層事前分布を用いた場合のシミュレーションを行った結果を示す第2の図である。
発明を実施するための最良の形態
以下、図面を参照して本発明の実施の形態について説明する。
なお、上述のとおり、脳内活動を外部から観測する手法としては、脳磁計(Magnetoencephalography:MEG)や脳波計(Electroencephalogram:EEG)がある。以下では、これらのうち、主にMEGについて説明することにする、ただし、MEGにおける磁場を電場に置き換えればEEGについても、本発明を適用することが可能である。
したがって、本明細書においては、本来「電場」と「磁場」とが複合して存在する「電磁場」に対して、「観測対象」となる物理量が「電場」である場合、または「磁場」である場合とを総称して、「電磁場を観測する」と呼ぶことにする。
[実施の形態1]
図1は、脳磁計システムの構成の一例を示す概念図である。
脳磁計12には、マルチチャンネルのSQUID磁束計(超高感度磁気測定器)アレイを備えており、被験者10の頭蓋表面での磁場を測定する。この脳磁計12での測定結果を受け取って、コンピューターシステム20が、測定結果を解析し、脳内の電流源の位置の推定処理を行う。
本発明は、このコンピュータシステム20でのソフトウェア処理に関するものである。ただし、処理の一部ないし全部は、処理時間の高速化のためにハードウェアにより行われる構成であってもよい。また、このソフトウェアは、特に限定されないが、CD-ROM(Compact Disc Read Only Memory)やDVD-ROM(Digital Versatile Disc Read Only Memory)のような記録媒体22上に記録されて、コンピュータシステム20にインストールされてもよいし、あるいは、通信回路を介して配信されてもよい。
[脳内電流源推定の原理]
以下、本発明の「脳内電流源推定方法」の詳細について、説明する前に、本発明の推定方法の原理、および概略についてまとめておく。
(曲面上の電流分布による電磁場の復元)
電流源を金属と磁性体とで作られた理想的な電磁シールド面で囲み込めば、シールド面の外側に電磁場が存在しなくなることは電磁シールドの原理として良く知られている.すなわち、シールド面上を流れる電流と磁性体を構成する微小磁石の集まりが作る電磁場により、電流源がシールド面の外側に作る電磁場が完全に打ち消されることになる。また、微小磁石は仮想的な渦電流で置き換えることができる。これは見方を変えれば、このシールド面に適当な電流を流してやれば電流源が作る電磁場と全く同じ電磁場をシールド面の外側に作れることを示している.逆にシールド面の外側にある電流源が作る電磁場はこのシールド面にどのような電流を流しても完全には復元できない.この事実を以後「電磁場復元の原理」と呼ぶ。
図2は、電流源が作る電磁場を適当な曲面上で観測している場合の概念図である。
図2に示すように、観測曲面と電流源の間に仮想曲面を用意すると、上に述べた「電磁場復元の原理」から仮想曲面上に適当な電流を流してやれば電流源が観測面上に作る電磁場を復元できることが分かる。
図3は、脳内の電流源と複数の仮想曲面との関係を示す概念図である。
図3に示すとおり、電流が作る電磁場は距離の二乗に逆比例して減衰することを考慮すれば、電流源と等価な仮想曲面(仮想曲面2)上の電流の広がりは仮想曲面が電流源から遠ざかるほど広くなってゆく。したがって、仮想曲面1上の電流の広がりは、仮想曲面2上の電流の広がりよりも大きい。
一方、電流源に対して観測曲面と反対側にある仮想曲面(仮想曲面3)では電流源が観測面に作る電磁場を完全に復元する事が出来ない。
本発明では、以上の原理を元に観測面上での電磁場の観測データから電流源を推定する。
(電流源推定の原理)
図2に示すように、脳内で発生した電流源から作られる磁場(または電場)を脳の表面近くの観測面で観測していると仮定する。
脳内に仮想曲面を考え、観測された磁場を復元する仮想曲面上の電流分布を求める。この仮想曲面を脳表面から半径を縮めながら脳の内部へ移動させてゆくと、観測磁場を復元する電流分布は、その広がりが小さくなってゆき真の電流源を仮想曲面が含むときに最小になる。また仮想曲面を電流源より深く移動させていくと仮想曲面上での電流分布の広がりは再び広がってゆき、この電流が生ずる磁場と観測磁場との誤差も大きくなってゆく。
すなわち、観測磁場を復元する電流分布の広がりと磁場の復元誤差を見ることにより電流源の深さが特定できる。またこのようにして特定された深さの仮想曲面上の電流分布を求めれば電流源の広がりもわかることになる。以上が本発明の原理である。
(ベイズ推定法による電流源の深さの特定)
上述した電流源推定の原理を具体的に実施するにあたり、本発明では、以下のような確率論に基づく手続きをとる。その手続きの概略を以下にまとめることにする。
MEGやEEGで観測できるのは脳表面近くの数十カ所から数百カ所の磁場(電場)である。また仮想曲面上の電流分布を近似するために、仮想曲面上に格子点を設定し各格子点に電流双極子(または適当な電流源モデル)を割り当てる。電流分布を高い分解能で推定するためには格子点の数を増やして格子点の密度を大きくする必要がある。
一方、仮想曲面上の格子点の数を観測点の数より増やすと解が一意に決まらない。また観測点より推定点の数の方が多いと、観測データ数よりも推定すべきパラメータの数の方が多くなるので、電流源より深い位置にある仮想曲面でも観測磁場の復元性が良くなる。
このような困難を解決するために、本発明ではベイズ推定理論を利用して仮想曲面上の電流分布を推定する。ベイズ推定を行う際に、電流源に関する事前知識を表す事前分布を用いる。具体的には、脳内電流源は複数の場所に局在していると考えられるので、電流分布の広がりができるだけ小さく推定されるような事前分布、すなわち、観測データだけからではその強さが良く特定できないような格子点上の電流双極子は、その強さがゼロに近くなるような事前分布を導入する.これは有効信号路自動決定(Automatic Relevance Determination)事前分布(文献:R.M.Neal,Bayesian Learning for Neural Networks,Springer Verlag.,1996)と呼ばれる階層的事前分布を導入することで実現できる。この事前分布を以後、「局在条件事前分布」と呼ぶことにする。また、局在条件以外の事前知識を導入した事前分布の選び方については実施の形態2において説明をすることにする。
しかし、局在条件事前分布と観測データから事後確率分布を解析的に計算することは不可能である.そこで、本発明では、後に説明するように、事後確率分布を近似的に計算する手法として変分法的ベイズ推定法(Variational Bayes method)(文献:Attias,H.Proc.15th Conference on Uncertainty in Artificial Intelligence pp.21-30,1999およびMasa-aki Sato,Naural Computation,13,1649-1681,2001)を用いる。ただし、モンテカルロ法等の他の近似手法を用いることも可能である。
この局在条件事前分布を用いたベイズ推定法を行うことにより、観測データを復元し、かつ出来るだけ広がりの少ない仮想局面上の電流分布が得られる。また深さの異なる仮想曲面を使って推定を行ったときのモデル事後確率を比べることにより電流源の深さが特定できる。
モデル事後確率の対数は対数尤度項とモデル複雑度項の和で表される。対数尤度項は復元誤差が小さいほど大きくなる。
一方、モデル複雑度項は有効な(すなわち適当なしきい値以上の強さを持つ)格子点上の電流双極子の数が少ないほど大きくなる.上述したように、復元誤差と仮想曲面上の電流分布の広がり(有効な電流双極子の数)は仮想曲面の深さが電流源と一致したときに最小になるので、この時モデル事後確率が最大になることが分かる。すなわちモデル事後確率が最大になる深さの仮想曲面上に電流源が存在することが分かる。
すなわち、以上をまとめると、本発明は、脳磁波(MEG)や脳波(EEG)の観測データから脳内電流源の位置を深さ方向まで含めて推定する方法を提供する。
本発明の基本的な考え方は、上述のとおり、電流源と観測面の間にある曲面上に適当な電流を流すことによって、電流源が発生する電磁場を復元することができ、またこの曲面が電流源に近づくほど曲面上の電流分布の広がりが小さくなるという点に基づく。この考え方をベースに、本発明は、観測データを復元する曲面上の電流分布に対するベイズ推定を行う際、モデル事後確率は曲面が電流源を含むときに最大になること、すなわちモデル事後確率を調べることにより電流源の位置を深さ方向まで含めて推定することを1つの特徴とする。
(複数個の電流源の存在する場合)
以上では、電流源が一つの場合について説明したが、この方法は電流源が複数ある場合にも適用できる。
電流により生ずる電磁場は距離の二乗に逆比例して減衰するので、脳表面に最も近い電流源が脳表面の観測磁場に最も影響を与える.そこで脳表面に近い電流源から順番に特定することが出来る。
仮想曲面を脳表面からだんだん深くしてゆくと脳表面に最も近い電流源(これを第1電流源と呼ぶ)の近くでモデル事後確率が極大になる.電流源が2つ以上ある場合にはこの仮想曲面上の電流分布には電流源の数に対応した局所的な電流双極子の集合が複数個出来る。
この局所的な電流双極子の集合を局在電流分布と呼ぶ。個々の局在電流分布を含む局所面を切り出し、これを深さ方向に動かしてモデル事後確率を求める。第1電流源に対応した局所面を動かしたときにはモデル事後確率が第1電流源の深さで極大になる。しかし、これ以外の局所面を動かしてもモデル事後確率は第1電流源の深さで極大になることはない。このことから第1電流源の位置、すなわち最も表面に近い電流源の位置が特定できる。
2番目に深い電流源を求めるために、第1電流源に対応する局所面を固定し、残りの局所面の深さを同じにして深くしてゆく。すると今度は2番目に深い電流源のところでモデル事後確率が極大になる。ここでまた個々の局所面を深さ方向に動かすと第2電流源に対応する局所面を動かしたときにのみ第2電流源の深さでモデル事後確率が極大になる.こうして第2電流源の位置を特定できる。3番目以降の電流源も同様にして順次特定することが出来る。
この方法は個々の局所面の深さを独立に動かして調べる方法に比べて、計算時間がずっとすくなくて済むという利点を持つ。
(逐次的に解像度を上げる方法)
以上のようにして、脳表面に近い側から各電流源の位置を特定することが可能であることに加えて、本発明の脳内電流源推定方法では、電流源推定の解像度を逐次的に上げることが出来る。
すなわち、まず低い解像度で(仮想曲面上の格子点の数を少なくし、深さ方向の標本点の数も少なくして)電流源の位置を大雑把に推定する。この段階で各電流源に対応した局所面の位置が大体定まっている。
次により高い分解能で推定を行う。この時、局所面の面積は元の仮想曲面の面積に比べて小さくなっているので、同じ数の格子点を使っても分解能は高くなる。深さ方向の分解能もあげて電流源の位置の推定精度を上げることが出来る。分解能を上げたときに電流分布がより局在する場合には局所面をより小さくして再び分解能を上げて推定できる。
一方、分解能を上げても電流分布の広がりが余り変わらなければ、この電流源が広がっていることを示している。このようにして電流源の広がり具合に応じて分解能を調節することが可能である。
[電流源推定の具体的手続き]
以下では、以上概説した手続きにしたがって、具体的に「脳内電流源」の位置を特定する手続きについて、さらに詳しく説明する。
[(I)電流源推定のための準備]
以下に説明する電流源推定のための処理手順を行なうに当たって、まず、その準備となる手続きについて説明する。
すなわち、電流分布の推定のために、以下の手続きを予め行っておく必要がある。
(I-1) 仮想曲面形状と電流モデルの決定
(I-2) さらに、仮想曲面を移動させてその仮想的な電流分布を推定するに当たり、仮想曲面上のサンプル点と、深さ方向のサンプル地点の決定を行なう必要がある。
(I-3) または、以下で詳しく説明するように、電流分布推定のための階層事前分布の分布形状を指定するためのメタパラメータを初期値として予め決定しておく必要がある。
以下では、上述した(I-1)仮想曲面形状と電流モデルの決定の手続きについて、さらに詳しく説明する。
(I-1-1) 仮想曲面形状の決定
最も単純な仮想曲面形状は、脳を半球とみなし、半径の異なるいろいろな半球面を仮想曲面とすることにより得られるものである。
核磁気共鳴画像(Magnetic Resonance Imaging:MRI)等の他の測定方法により、大脳皮質等の脳内電流が存在する可能性のある場所がわかっている場合には、それらの情報をもとに仮想曲面形状を決定することもできる。
特に高精度で大脳皮質の形状がわかっていて、かつ、その他の領域には脳内電流源が存在していないと考えられる場合は、以下で説明する深さ方向の探索を省略して大脳皮質点の電流分布推定を行うこともできる。また、一部の領域だけ深さ方向の探索をすることも可能である。
この場合は、深さの異なる仮想曲面の形状が一般に異なっていてもよいが、互いに交わらないように形状を決めておく必要がある。
(I-1-2) 電流モデルの決定
仮想曲面上の電流モデルとして、ここでは電流双極子モデルを考える。ただし、他の電流モデルを考えることも可能である。
仮想曲面上に適当な格子点(サンプル点){Xn|n=1,…,N}を考える。各格子点上に電流の強さがjnの電流双極子を考える。このとき、格子点Xnの電流双極子jnが脳表面上での観測点Yi(i=1,…,I)に作る磁場は、以下のビオサバールの式で表わされる。
JP0003730646B2_000002t.gifここで、μは透磁率であり、たとえば、ベクトルXnに対し、記号|Xn|は、ベクトルXnの絶対値を表す。より正確な式として、脳を導電性溶液がつまった球体とみなした時のサルバスの式(文献:Sarvas J.Phys.Med.Biol.32,11-22,1987)、あるいは脳の細かな構造まで考慮に入れて有限要素法や境界要素法等の数値解を用いることもできる。以下では、説明や表記を簡単にするために、上記のビオサバールの式を用いて説明する。
このとき、仮想曲面上の電流双極子{jn|n=1,…,N}が観測点Yi(i=1,…,I)に作る磁場は、次式で表わされる。
JP0003730646B2_000003t.gif観測点Yiで観測できる磁場の方向をベクトルSi,c(c=1,…,C)とすると、この点での磁場のベクトルSi,c方向成分Bi,cは、次式のようになる。
JP0003730646B2_000004t.gifまた、観測磁場として磁場の局所的勾配を測定する場合には、上式のYiによる微分を観測することになる。
格子点Xn(位置ベクトル)での電流双極子の方向に制限をつける場合には、電流双極子の可能な独立方向をbn,d(d=1,…,D)として、電流双極子は以下の式で表わされる。この式において、D=3の場合は方向に制限がない場合に対応している。
JP0003730646B2_000005t.gif以上まとめると、仮想曲面上の格子点{Xn|n=1,…,N}での電流双極子が観測点Yに作る磁場は次式のように書ける。
JP0003730646B2_000006t.gifここで、jn,dは、格子点Xnでの電流双極子の独立な成分である。また、関数G(i,c;n,d)は格子点Xnでの電流双極子jn,dが観測点Yiに作る磁場のSi,c方向成分である。
電流分布推定の問題とは、観測された磁場{Bi,c|i=1,…,I;c=1,…,C}から仮想曲面上の電流分布{jn,d|n=1,…,N;d=1,…,D}を推定する問題であるといえる。
上述した観測点における磁場の測定値を、行列表現を用いて表わすこととすると、以下のように書くことができる。なお、Gはリードフィールド行列と呼ばれる。ビオサバールの式ではなく、サルバスの式、あるいは境界要素法等の数値解を用いる場合には、このリードフィールド行列の値が異なるだけであり、以下の説明のその他の部分は同様に適用することができる。
JP0003730646B2_000007t.gif(電流源確率モデル)
このように電流モデルを決定した上で、さらにこのような電流分布の推定のために、以下に説明する「電流源の確率モデル」を導入する。
(I-1-3) 電流源確率モデル
上述した仮想曲面上の電流モデルに対する確率モデルを以下のように考える。観測される磁場は仮想曲面上の電流分布Jから作られる磁場と観測ノイズの和として表わされるとする。また、観測ノイズは各測定点で独立な分数σ2を持つガウスノイズとして表わされるとする。
より一般的には、各測定点のノイズの相関が共分散行列の形で表される多次元正規分布に従うガウスノイズを考えることもできるが、以下では説明の簡単のために等方等分散のノイズモデルを用いて説明する。
すなわち、観測磁場が以下の式のように表わされるものとする。
JP0003730646B2_000008t.gifこの電流モデルに対する確率分布は以下のように与えられる。
まず、ある特定の深さに対する仮想曲面、あるいは各局所面の深さが指定された局所面の集合をMで表わす。この仮想曲面M上の電流分布Jが与えられたときに、観測される磁場がBである確率P(B|J,β,M)は、以下の式で表わされる。ただし、以下の式において、β=1/σ2である。
JP0003730646B2_000009t.gif(I-1-4) 階層事前分布
仮想曲面M上の電流分布Jに対する事前分布として、上述したように、局在条件階層事前分布を仮定することにする。
観測を行なう前の電流分布Jに対する事前分布(Jが実現する確率)は局在条件階層事前分布の仮定の下で以下の式のように表わされる。
JP0003730646B2_000010t.gifただし、Γ(…)は、ガンマ分布を表わし、次式のように定義される。次式において、Γ(γ0)は、ガンマ関数であり、その定義式も併せて示す。
JP0003730646B2_000011t.gifまた、αとτは、電流分布Jとノイズ逆分散βに対する事前分布をモデル化するために導入されたハイパーパラメータであり、このαとτに対する階層事前分布P0(α|M)とP0(τ|M)は以下のように与えられる。
JP0003730646B2_000012t.gifなお、上式において、αとτに対する階層事前分布P0(α|M)とP0(τ|M)の分布形状を決定するメタパラメータ バーα0、γα0、τ0、バーγτ0については、後に詳しく説明する。なお、変数名の前に「バー」との語句をつけたものは、その変数の頭部に記号「-」がついていることをあらわす。
(I-1-5) ベイズ推定
磁場のデータBを観測したときに、電流分布がJである事後確率分布P(J|B,M)は、ベイズの定理を使って以下のように求まる。
JP0003730646B2_000013t.gifここで、事後確率分布P(J,β,α,τ|B,M)は、以下の式で表される。
JP0003730646B2_000014t.gifこの事後確率分布を用いて、電流分布の期待値E[J|B,M]は以下のように与えられる。
JP0003730646B2_000015t.gifまた、P(B|M)は仮想曲面モデルMの周辺尤度である。電流源の深さを推定するとき、深さの異なるいくつかのモデルを比べる。このとき、深さに関する事前の知識がないものとする。すなわち、以下では、P(M)=一定であるものとする。
観測データBが与えられたときに、モデルMが真のモデルである確率、すなわち、モデル事後確率P(M|B)は、上の仮定のもとでモデル周辺尤度P(B|M)に比例し、以下の関係が成り立つ。
JP0003730646B2_000016t.gif(I-1-6) 変分法的ベイズ推定
確率モデルと階層事前分布が与えられたときに、モデル周辺尤度を解析的に求めることは一般的にはできない。
そこで、モデル周辺尤度を近似的に計算する手法として、変分法的ベイズ法を用いる。モンテカルロ法や、ラプラス近似法等、他の近似方法を使うことも可能である。
「変分法的ベイズ法」について、具体的な手続きは後述するものとして、以下に簡単に説明しておく。
真の事後分布P(J,β,α,τ|B,M)を近似的に求めるために、試験事後分布Q(J,β,α,τ)を考える。
2つの確率分布P(J,β,α,τ|B,M)とQ(J,β,α,τ)の近さは、次式で表されるようなK-L距離を用いて計算することができる。
JP0003730646B2_000017t.gifK-L距離は、2つの分布が等しいときに限りゼロになり、それ以外のときは常に正の値をとる。
上式において試験事後分布Qに対する自由エネルギーF(Q)を定義すると、以下の不等式が得られることになる。
JP0003730646B2_000018t.gifすなわち、自由エネルギーF(Q)を最大にする分布Q(J,β,α,τ)は、真の事後分布P(J,β,α,τ|B,M)に等しくなり、このとき自由エネルギーは対数周辺尤度に等しくなる。
変分法的ベイズ法では、Q(J,β,α,τ)の関数形を以下の形に制限して自由エネルギーの最大化を行なう。
JP0003730646B2_000019t.gif上式の右辺第2項Qαを固定して、右辺第1項QJに関してF(Q)を最大化するステップと、第1項QJを固定して、第2項Qαに関してF(Q)を最大化するステップを交互に繰返すことにより、自由エネルギーF(Q)を極大にする分布Q*が得られることになる。
[(II)電流源推定の手順]
以上のような準備の下に、以下電流源を推定する具体的な手続きについて、図にしたがって説明する。
図4は、本発明の脳内電流源の位置の推定方法の手順の全体的な流れを示すフローチャートである。
図4を参照して、まず、脳内電流源の位置の推定処理が開始されると(ステップS100)、上述したように、電流分布推定のための階層事前分布の分布形状を指定するためのメタパラメータの値と推定すべき変数の初期値を決定する(ステップS102)。
続いて、初期分解能を用いた電流源の初期推定が行われ、電流源の候補の抽出が行われる(ステップS104)。
さらに、このようにして抽出された電流源の候補のうち、まず、最も脳表面に近い電流源の位置の特定が行われる(ステップS106)。その後、他の電流源の深さの特定が順次行われる(ステップS108)。
以上のようにして、初期分解能による電流源の位置の特定が行われた後、空間分解能を上げて、電流源の位置の再推定が行われ(ステップS110)、処理が終了する(ステップS112)。
以下、図4で説明した各ステップの処理について、さらに詳しく説明する。
(II-1) 初期分解能を用いた電流源の初期推定(電流源の候補の抽出)
図5は、図4で説明した処理のうち、ステップS104の処理、すなわち、初期分解能を用いた電流源の初期推定の処理を説明するためのフローチャートである。
図5を参照して、まず、初期分解能での深さ方向のサンプル地点を{Rk|k=1,…,K}とする。各深さRkにおける仮想曲面に対して電流分布の推定を行なう。さらに、このようにして推定された電流分布に基づいて、各深さRkに対する事後確率を求める。この事後確率が後に説明する自由エネルギー値に相当する(ステップS202)。
以上のようにして求められた事後確率が最大になる深さRMでの電流分布に対して、電流の強さの極大点を求める(ステップS204)。
さらに、各極大点を含む局所面を決定する。局所面の数をL個とすると、各局所面が局在化された電流源の候補となる(ステップS206)。
以下、図4中のステップS102の処理および、それに続く図5中のステップS202の処理について、さらに詳しく説明する。
(II-1-1) 電流分布推定のための階層事前分布の分布形状を指定するためのメタパラメータの値と推定すべき変数の初期値の決定
上述したとおり、電流源の推定は、変分法的ベイズ法を用いて行われる。したがって、図4のステップS102の処理は、このような変分法的ベイズ法の適用にあたって、以下のように行われる。
階層事前分布の分布形状を指定するメタパラメータを、たとえば、以下のように決定する。
JP0003730646B2_000020t.gifさらに上述した式に基づいて、各変数の初期設定を以下にように行なう。
JP0003730646B2_000021t.gif(II-1-2) 電流分布推定のための変分法的ベイズ法の具体的処理
(1) パラメータJ,βの期待値の計算(J-ステップの処理)
ここでは、パラメータJ,βの期待値の計算を行なう。この処理により、QJに関して自由エネルギーF(Q)を最大化する処理が行なわれる。
対角行列Aを以下のように定義することで、QJが以下の式に基づいて導出される。なお、以下の式において、変数の期待値には、その変数名の上に、”-”をつけて表現している。
JP0003730646B2_000022t.gif(2) ハイパーパラメータ期待値、すなわち、パラメータα,τの期待値の計算(α-ステップの処理)
Jステップに続いて、ハイパーパラメータα,τの期待値の計算を行なう。すなわち、この処理では、Qαに関して自由エネルギーF(Q)を最大化する処理を行なうことになる。
この手続は以下の式のように表わされる。
JP0003730646B2_000023t.gif(3) 自由エネルギーの計算
以上説明したようなJステップおよびαステップに基づいて計算されたQJおよびQαを用いて、自由エネルギーを以下のようにして計算する。
JP0003730646B2_000024t.gifこのようにして、自由エネルギーFの値が収束するまで、上述したJステップの処理から自由エネルギーの計算までの処理を繰返す。収束後の自由エネルギーF(Q)の値が、対数周辺尤度logP(B|M)の近似値を与える。
また、対数モデル事後確率は、logP(B|M)と定数分だけの違いしかないので、モデル事後確率が最大になるモデルは、上述したような近似の下では自由エネルギーが最大になるモデルに等しいことになる。
以上のような手続により、深さRkに対する事後確率を求めることができる。これに応じて、上述した図5中のステップS204およびS206の処理を行なうことにより、各局所面に局在化された電流源の候補を求めることができる。
(II-2) 最も脳表面に近い電流源の特定
次に、図4のステップS106の処理、すなわち、以上のようにして求められた電流源の候補から、まず、最も脳表面に近い電流源の特定を行なう処理について説明する。
図6は、最も脳表面に近い電流源の特定を行なう処理を説明するためのフローチャートである。
まず、初期値として変数lの値を1とする(ステップS302)。つづいて、変数lと、変数lのとり得る最大値Lとの比較が行われ(ステップS304)、変数lが最大値L以下であれば、処理はステップS306に移行し、変数lが最大値Lを超えていれば、処理はステップS324に移行する(ステップS304)。すなわち、ステップS306からS322までの処理を、l=1からLまで繰返す。
ステップS306では、l番目の局所面以外の局所面の深さを図5のステップS204で求められた深さRMに固定する。
変数kの値を1とする(ステップS308)。つづいて、変数kと、変数kのとり得る最大値Kとの比較が行われ(ステップS310)、変数kが最大値K以下であれば、処理はステップ312に移行し、変数kが最大値Kを超えていれば、処理はステップS320に移行する。
ステップS312では、l番目の局所面の深さをまずRkとする。
次に、L個の局所面の集合に対して電流分布の推定を行なう(ステップS314)。
さらに、この局所面配置に対する事後確率(自由エネルギー)を求める(ステップS316)。変数kの値を1だけインクリメントして(ステップS318)、処理はステップS310に復帰する。
このようにして、k=1からk=Kまでの深さの各局所面についての処理を行なった後、事後確率が最大になるl番目の局所面の深さRM(l)を求める(ステップS320)。変数lの値を1だけインクリメントして(ステップS322)、処理はステップS304に復帰する。
このようにして、各変数lについて求めた局所面の深さRM(l)の中で、最も脳表面に近い(浅い)lを求め、これをlとする(ステップS324)。
以上の処理により、l番目の局所面が最も脳表面に近い電流源に対応し、その深さの初期推定値がRM(ll)として求められ、処理が終了する(ステップS326)。
(II-3) 各局所面に対応する電流源の深さの特定
図7および図8は、図4のステップS108の処理、すなわち、各局所面に対応する電流源の深さの特定処理を説明するためのフローチャートである。
図7を参照して、初期値として変数sの値を1とする(ステップS402)。つづいて、変数sと、変数sのとり得る最大値Lとの比較が行われ(ステップS404)、変数sが最大値L未満であれば、処理はステップS406に移行し、変数sが最大値L以上であれば、処理はステップS434に移行する(ステップS404)。すなわち、ステップS406からS432までの処理を、s=1からLまで繰返す。
ステップS406では、これまでに特定された局所面の深さを{RM(ll),…,RM(ls)}とするとき、これらの局所面の深さを、それぞれ{RM(ll),…,RM(ls)}に固定する(ステップS406)。
まず、lが{ll,…,ls}のいずれとも等しくなく、かつlは{1,…,L}の集合に属しているものとする(ステップS408)。次に、変数lの値として、可能な値を全て尽くしたか否かを判断し(ステップS410)、全て尽くしている場合は処理はステップS430に移行する。一方、全てを尽くしていない場合は、処理はステップS412に移行する。
ステップS412では、l番目の局所面以外で、かつ{ll,…,ls}でもない局所面の深さをRMに固定する。
変数kの値を1とする(ステップS414)。つづいて、変数kと、変数kのとり得る最大値Kとの比較が行われ(ステップS416)、変数kが最大値K以下であれば、処理はステップS418に移行し、変数kが最大値Kを超えていれば、処理はステップS426に移行する。すなわち、ステップS418からS424までの手続を、k=1からk=Kまで繰返す。
ステップS418では、1番目の局所面の深さをRkにする。
次に、L個の局所面の集合に対して電流分布推定を行なう(ステップS420)。
さらに、この局所面配置に対する事後確率(自由エネルギー)を求める(ステップS422)。
図8を参照して、ステップS426では、以上のようにしてk=Kまでの処理が終わった後に、事後確率が最大になるl番目の局所面の深さRM(l)を求める(ステップS426)。
{ll,…,ls}のいずれとも等しくなく、かつ{l,…,L}の集合に属している、まだ処理を行っていない他の変数lを選択し(ステップS428)、ステップS410に復帰する。
以上の手続により、処理を行った全ての変数lに対応するRM(l)の中で、脳表面に最も近いlを求め、これをls+1とする(ステップS430)。
さらに、sにs+1を代入し(ステップS432)、処理はステップS404に復帰する。
全ての可能なsについて処理が終了することで、各局所面に対応する電流源の深さの特定が終了する(ステップS434)。
(II-4) 分解能を上げて電流源の再推定
図9および図10は、図4のステップS110の処理、すなわち、空間分解能を上げて電流源の位置を再推定する処理を説明するためのフローチャートである。図9および図10を参照して、まず、初期分解能を用いて推定した電流源と対応する局面の番号を脳表面に近いものから1,2,3,…,Lとなるように番号の並べ替えを行なう(ステップS502)。次に、各局所面での電流分布の拡がりに対応して局所面を小さくする(ステップS504)。
図4のステップS108までの処理で求めた局所面の深さRM(l)をRMold(l)とする(ステップS506)。
新しい深さ方向の分解能と探索幅を、それぞれΔRおよび(KL・ΔR)とする。また局所面上での格子点の分解能も上げる(ステップS508)。
初期値として変数lの値を1とする(ステップS510)。つづいて、変数lと、変数lのとり得る最大値Lとの比較が行われ(ステップS512)、変数lが最大値L以下であれば、処理はステップS514に移行し、変数lが最大値Lを超えていれば、処理はステップS534に移行する(ステップS512)。すなわち、ステップS514からS532までの処理を、l=1からl=Lまで繰返す。
ステップS514では、l番目以外の局所面l′の深さをRMold(l′)に固定する。
その上で、さらに、k=0,±1,…,±KLについて、以下のステップS518からステップS526までの処理を行なう。
まず、変数kの値を0とする(ステップS516)。つづいて、変数kの絶対値と、変数kの絶対値のとり得る最大値KLとの比較が行われ(ステップS518)、変数kの絶対値が最大値KL以下であれば、処理はステップS520に移行し、変数kの絶対値が最大値KLを超えていれば、処理はステップS528に移行する(ステップS518)。
ステップS520では、l番目の局所面の深さを(RMold(l)+k・ΔR)とする。
次に、L個の局所面の集合に対して、現在の分解能を用いて、電流分布推定を行なう(ステップS522)。
さらに、この局所面配置に対する事後確率(自由エネルギー)を求める(ステップS524)。その上で、kの値を{±1,…,±KL}のうちの次の値とし、処理はステップS518に復帰する。
以上のような処理がk=±KLまで行なわれた後に、ステップS528では、事後確率が最大となるkを求め、これをkMとする。
その上で、RMold(l)にRMold(l)+kM・ΔRを代入する(ステップS530)。変数lの値を1だけインクリメントして(ステップS532)、処理はステップS512に復帰する。
以上のような処理が変数lの値が最大値Lとなるまで行われた時点で、分解能が最終分解能に達していれば(ステップS534)、処理を終了する(ステップS538)。
一方、分解能が最終分解能に達していなければ(ステップS534)、局所面上での格子点の分解能と深さ方向の分解能を上げて(ステップS536)、処理は、ステップS510に復帰する。
以上説明したような「脳内電流源推定方法」を用いれば、MEG(またはEEG)の観測データから脳内電流源の位置を深さ方向まで含めて推定することが可能となる。さらに、このような深さ方向の推定は、複数の電流源が存在する場合にも適用可能である。しかも、電流源が電流双極子のように局在している場合にも、広がりを持つ電流源の場合にも適用できるだけでなく、電流源の広がり具合を推定することも可能である。
また、図9および図10で説明したような処理を、処理分解能での推定の後に追加して行えば、位置の推定分解能を逐次的に上げることが可能となる。このことは、また、生理学的な知見や他の観測方法による観測データなどにより探索範囲を限定して調べることが可能なことも意味する。
[シミュレーション結果]
以下では、以上説明した脳内電流源推定方法に従って、モデル的に仮定した電流源位置を推定するシミュレーションを行なった結果について説明する。
図11は、人間の脳を、半径10.0(任意単位)の半球状であると仮定し、この半球の表面上で観測された動径方向の磁場分布を上面から見た上面図である。
図11においては、後に説明するように、中心から半径7.0の位置に、単一の電流源が置かれている場合の半球表面上での磁場分布を示している。以下のシミュレーションでは磁場の観測データにS/N比が0.1のノイズを加えている。
図12は、初期分解能を用いた電流源の初期推定結果を示す図である。
図12では、図4のステップS102~ステップS104の処理において説明したとおり、初期分解能を用いて電流源の初期推定を行なうために、深さ(半径R)をステップ状に変化させつつ、各深さの半球上において自由エネルギーが最大となるような電流分布を変分法的ベイズ法に基づいて算出した結果を示す。
図12(a)は半径R5.0の場合を示し、図12(b)は半径R=6.0の場合を示し、図12(c)は半径R=7.0の場合を示す。
図13は、さらに、半径がより大きい(脳表面により近く)場合の初期分解能を用いた電流源の初期推定結果を示す図である。
図13(d)は半径R=7.5の場合であり、図13(e)は半径R=8.0の場合であり、図13(f)は半径R=9.0の場合である。
図14は、このように、脳内に想定した様々な深さの仮想半球面に対して、自由エネルギーを最大とするような電流分布を各々求める際に得られる、各仮想半球面上での自由エネルギーの大きさを示す図である。
図14に示すとおり、半球面全体を仮想曲面として計算した場合、半径R=7から半径R=8の間に自由エネルギーの最大点が存在していることがわかる。
図13(d)に示すとおり、半径R=7.5の仮想半球面上には、電流分布に1つの極大点が存在する。
このようにして、初期分解能を用いて電流源が初期推定されたので、自由エネルギーが最大となっている仮想曲面において、電流分布の極大点について、その極大点を含む局所面を求め、さらにこの局所面について電流源の特定処理が行なわれる。
図15および図16は、このように、電流分布の極大点を含む局所面について、局所面上の格子点の分解能と深さ方向の分解能を上げて、さらに電流源の特定処理を行なった場合の処理を示す図である。
図15(a)は半径R=5.0の、図15(b)は半径R=6.0の、図15(c)は半径R=7.0の場合の局所面における電流分布を示している。同様に、図16(d)は、半径R=7.5の、図16(e)は半径R=8.0の、図16(f)は半径R=9.0の電流分布をそれぞれ示している。
図17は、このようにして局所面の深さを変えつつ自由エネルギーを求めた結果を示す図である。
図17に従えば、局所面を用いて計算した結果の自由エネルギーは、半径R=7付近に極大値を持ち、結果から、半径R=7の位置に1つの電流源が存在していることを特定することが可能である。
このときの電流源における電流分布は、図15(c)のとおりであり、上述したとおり、電流分布も最も拡がりの狭いものとなっていることがわかる。
図18は、さらに、脳表面を半球と仮定した場合に、電流源が2つ脳内に存在する場合の脳表面における磁場分布を示す上面図である。
図18に示した例においては半径R=6.0およびR=8.0の位置にそれぞれ電流源が置かれているものとする。
図19および図20は、初期分解能を用いた電流源の初期推定結果を示す図である。
すなわち、図19および図20では、2個の電流源が脳内に存在している場合に、図4に示した初期分解能を用いた電流源の初期推定を行なう処理(ステップS104)を説明している。
図19(a)は半径R=5.0の場合の、図19(b)は半径R=6.0の場合の、図19(c)は半径R=7.0の場合の仮想半径面上の電流分布を示している。
図20(d)は、半径R=7.5の場合の、図20(e)は半径R=8.0の場合の、図20(f)は半径R=9.0の場合の仮想半球面上の電流分布をそれぞれ示している。
図21は、このようにして仮想半球面全体について、各仮想半球面ごとに自由エネルギーが最大となるように電流分布を求めた場合の、自由エネルギーの半径依存性を示す図である。
図21に示すとおり、半球面全体にわたって自由エネルギーを計算した場合、半径R=7から半径R=8の間の半径7.5付近に自由エネルギーの最大値が存在していることがわかる。
図20(d)に示すとおり、半径R=7.5の仮想半球面上には、電流分布に2つの極大点が存在し、2つの局所面が決定される。
次に、図4におけるステップS106の処理に対応する計算結果について説明する。
図22および図23は、このようにして、最も脳表面に近い電流源を特定するために、事後確率が最大になるl番目の局所面の深さRM(l)を求める処理を説明するための図である。
したがって、図22および図23は、半球面全体にわたって計算した自由エネルギーが最大となる面上に他方の局所面を固定したままで、一方の電流源に対応する局所面を移動させた場合の電流分布を示す図である。
図22(a)は半径R=5.0の場合の、図22(b)は半径R=6.0の場合の、図22(c)は半径R=7.0の場合の局所面上の電流分布を示している。
図23(d)は、半径R=7.5の場合の、図23(e)は半径R=8.0の場合の、図23(f)は半径R=9.0の場合の局所面上の電流分布をそれぞれ示している。
図24は、図22および図23で示したように仮想局所面を移動させた場合の自由エネルギーの局所面位置(半径)依存性を示す図である。
図24に示すとおり、半球面全体にわたって計算した自由エネルギーが最大となる面上に他方の局所面を固定したままで、一方の電流源に対応する局所面を半径R=8の位置にしたときに、自由エネルギーが最大となる。
したがって、図24に示した結果から、最も脳表面に近い電流源の位置としては、半径R=8.0の位置であることがわかる。
続いて、他方の側の電流源の深さの特定を行なう。
この場合、一方の電流源に対応する局所面の深さは、半径R=8.0に固定しておき、他方の電流源に対応する局所面の深さを変化させていく。
図25および図26は、特定された面上に一方の局所面を固定したままで、他方の電流源に対応する局所面を移動させた場合の電流分布を示す図である。
図25(a)は半径R=5.0の場合の、図25(b)は半径R=5.5の場合の、図25(c)は半径R=6.0の場合の局所面上の電流分布を示している。
図26(d)は、半径R=6.5の場合の図26(e)は半径R=7.0の場合の、図26(f)は半径R=7.5の場合の局所面上の電流分布をそれぞれ示している。
図27は、図25および図26で示したように仮想局所面を移動させた場合の自由エネルギーの局所面位置(半径)依存性を示す図である。
図27に示すとおり、この場合半径R=6の位置で自由エネルギーが極大となることがわかる。
このような処理により、2つの電流源が脳内に存在している場合も、各々の電流源の深さを特定することが可能であることがわかる。
[実施の形態2]
以下では、上述した実施の形態1における「脳内電流源推定法」、「脳内電流源推定プログラム」および「脳内電流源推定装置」の構成において、ソフトウェア(または、一部ハードウェアにより処理されてもよい)により実現される処理のうち、「(I)電流源推定のための準備」から「(3)自由エネルギーの計算」までの処理を変更した本発明の実施の形態2について説明する。
まず、実施の形態2における変更点を説明する前提として、実施の形態1の処理の概念を改めて以下にまとめるとともに、その基本概念のうち、実施の形態2において変更される概念を説明しておく。
(実施の形態1における仮想曲面推定の概念)
実施の形態1においては、仮想曲面上の電流モデルとして、格子点における電流双極子を想定し、観測面上の観測磁場Bから電流分布Jを推定するという逆問題を解くという処理を行った。このとき、このような逆問題は、一般に、(電流双極子の存在する格子点数)が(観測点数)よりも多いために、いわゆる「不良設定問題」となり、厳密に解を求めることが困難である。
実施の形態1では、このような問題を解くために、「ベイズ推定」の手法を用いた。
ベイズ推定においては、仮想曲面Mと観測磁場Bが与えられた時の事後確率分布が以下のように表現されることを利用していた。
JP0003730646B2_000025t.gifここで、P(J|B,M)は、磁場Bが観測されたと言う条件での、仮想曲面M上での電流の確率分布であって、「事後確率分布」を表す。また、P(B|J,M)は「データ尤度」であり、P(J|M)は事前分布であり、P(B|M)は周辺尤度である。
周辺尤度は、以下の式で表される。
JP0003730646B2_000026t.gif事後確率P(M|B)は、以下の式で表される。
JP0003730646B2_000027t.gif深さに関する事前の知識がないものとすると、P(M)=一定となるので、実施の形態1でも説明したとおり、観測データBが与えられたときに、モデルMが真のモデルである確率、すなわち、モデル事後確率P(M|B)は、上の仮定のもとでモデル周辺尤度P(B|M)に比例する。
すなわち、実施の形態1において、電流モデルの選択にあたっては、事後確率、すなわち、周辺尤度が最大となるモデルを選ぶという処理を行った。
このような処理を行うために、以下の式で表される対数周辺尤度logP(B|M)を計算の対象とする。
JP0003730646B2_000028t.gif上式に示すとおり、対数周辺尤度を最大にするということは、概念的には、誤差と実効的なパラメータ自由度とを最小にすること、言いかえると、誤差と電流分布の広がりを最小とすることに相当する。
実施の形態1においては、上記のような対数周辺尤度の最大化の問題を、試験事後分布Qを事後分布の近似として導入し、自由エネルギーF(Q)を最大化するように事後分布を決定するという、いわゆる「変分法的ベイズ法」を用いて解いていた。
(実施の形態1における事前分布の概念)
実施の形態1においては、脳内電流源に対する事前知識として、電流源は脳内の複数の場所に局在しているという局在条件を用い、これを階層的事前分布の形で表わした。すなわち、ハイパーパラメータαn(n=1,…,N)を導入して、次式の階層事前分布の形を仮定した。
JP0003730646B2_000029t.gif上式からハイパーパラメータαnは電流Jnの逆分数(分散の逆数)に比例していることがわかる。もし、電流Jnの分散とノイズの逆分散βの値が予めわかっている場合には、αnの値を電流の逆分散とノイズの逆分散の値とから計算される値に固定した事前分布を用いればよい。しかし、電流の分散の値は通常わからないので、実施の形態1においては、ハイパーパラメータαnに対する階層事前分布P0(αn|バーα0,γα0)を仮定してαnに対するベイズ推定を行なっていた。
ここで、バーα0とγα0は、階層事前分布の分布形状を決めるメタパラメータであり、すべてのαnに対して共通の値を仮定していた。この局所条件階層事前分布を用いてベイズ推定を行なった場合、前述したようにできるだけ少ない数の電流双極子で観測データを復元するようなベイズ推定が行なわれる。
一方、MEGの電流源推定においてしばしば用いられる最小ノルム推定法(文献:Hamalainen M.S.and Ilmoniem R.J.,Med.& Biol.Eng.& Comput.,32,35-42,1994)は、ベイズ推定の立場からは、事前分布としてαnの値をすべての点において共通の値にし、しかもこの値を直接適当な値に定めることに相当する。
最小ノルム推定では、すべての点における電流強度の総和を最小にしようとするので、推定解には小さな電流強度を持つ点が多く現れる。このため、これらの小さな電流分布が真の電流分布を与えているのか、それとも推定誤差のためのノイズとみなしてよいのかを判定することが難しいという欠点を持っている。
局在条件階層事前分布は、最小ノルム推定法のこのような欠点を解消してくれる。
(実施の形態2における変更点)
実施の形態2においては局所条件のみならず連続条件も考慮に入れた事前分布を仮定する。空間分解能がミリメートル単位の高精度の電流源分布推定を行なう場合には、電流分布の連続性を考慮する必要がある。実際、大脳皮質の神経細胞は半径10mm程度のコラム構造を形成しており、神経活動が起こる場合には半径10mm程度の領域内における神経細胞が同時に発火する確率が高い。このような点を考慮して、以下に説明するように実施の形態2においては、局所条件と連続条件を組合せた階層事前分布を仮定する。
連続条件を導入するために、新たに内部変数Zを導入し、電流Jは内部変数Zを平滑化フィルタWで平滑化したものとして表現する。ここで、このような平滑化フィルタとしては、特に限定されないが、たとえばガウスフィルタ等を用いることができる。このとき、事前分布は以下のように表わせる。
JP0003730646B2_000030t.gif上記の事前分布においてZに対する積分は容易に実行できて、以下のように書ける。
JP0003730646B2_000031t.gif但し、AとΛは、それぞれ{αn|n=1,…,N}と{λn|n=1,…,N}を対角成分に持つ対角行列である。
すなわち、内部変数Zを導入して事前分布P0(J|Z,α,β)P0(Z|λ)を仮定することは、電流Jに対する共分散行列β-1(A-1+W・Λ-1・W′)を持つ事前分布を仮定することに等しい。但し、内部変数Zを用いて事前分布を表わした方が変分法的ベイズ法を適用しやすいので、以下では内部変数Zを用いて事前分布を表わす。また、ハイパーパラメータαnと同様に、ハイパーパラメータλnの値も事前にはわからないので、以下の形の階層事前分布を仮定して、αnとλnに対するベイズ推定を行なう。
JP0003730646B2_000032t.gifこのように拡張された階層的事前分布は、その特別の場合として実施の形態1の局在条件事前分布を含む。すなわち、平滑化フィルタ行列Wが恒等的に0の場合、上記の連続条件と局在条件を組合せた階層的事前分布は、実施の形態1の局在条件事前分布に一致する。
また、階層的事前分布の分布形状を決めるメタパラメータ バーα0n,γα0n,バーλ0n,γλ0nは、より一般的に格子点の場所に依存してもよいと仮定している。MEG(あるいはEEG)の観測データしかない場合は、これらの値に対する事前の知識は特にないので、後述するように、実施の形態1と同様に、実質的にはこれらのメタパラメータを格子点の場所に依存せずに共通に決める。しかし、実施の形態3で説明するように、他の観測手段による観測データが利用できる場合には、それらの条件を用いてメタパラメータの値を格子点毎に決めることができる。
(実施の形態2の具体的な計算手続き)
以下で説明する計算手続きは、その大筋において、実施の形態1と同様である。したがって、以下では、実施の形態1からの変更点を主として説明することにする。
まず、実施の形態1と同様にして、仮想曲面上の電流分布Jが与えられたときに、観測される磁場がBであるモデル確率分布P(B|J,β)は、対数表現を用いると、以下の式で表わされる。
JP0003730646B2_000033t.gif[階層事前分布]
実施の形態2の階層事前分布は、以下の式で表される。
JP0003730646B2_000034t.gifここで、βは実施の形態1と同様に、ノイズの逆分散を示す。
[試験事後分布の計算]
このとき、試験事後分布は、Q=QJ(J,β)QZ(Z)Qα(α,λ,τ)のように、QJとQZとQαとの積の形に表されると仮定する。
そして、自由エネルギーF(Q)を順番にQJとQZとQαの各々に関して最大化する。まず、第1ステップで、QZとQαとを固定して、QJに関してF(Q)を最大化すると、以下の式が得られる。
JP0003730646B2_000035t.gif次に、第2ステップで、QJとQαとを固定して、QZに関してF(Q)を最大化すると、以下の式が得られる。
JP0003730646B2_000036t.gifさらに、第3ステップで、QJとQZとを固定して、Qαに関してF(Q)を最大化すると、以下の式が得られる。
JP0003730646B2_000037t.gifなお、上式における未知変数の具体的計算方法については、以下に説明する。
実施の形態1と同様にして、実施の形態2においても、繰り返し計算を用いて自由エネルギーF(Q)を極大化する試験事後分布を求める。
具体的には、以下に説明する手続きを用いて、上述した試験事後分布Qを以下の第1~第3のステップで計算し、自由エネルギーF(Q)の値が収束するまでこの手続きを繰り返す。なお、その際の定数項の決め方については、後に説明する。
(1)第1のステップ(J-ステップと呼ぶ)の手続き
電流Jと逆分散βの期待値および変数γβについて以下の式に従い計算する。
JP0003730646B2_000038t.gif(2)第2のステップ(Z-ステップ)の手続き
続いて、内部変数Zの期待値は以下のとおり計算される。
JP0003730646B2_000039t.gif(3)第3のステップ(α-ステップ)の手続き
続いて、ハイパーパラメータαn、γαn、λn,γλn、γτについて、以下の計算を行う。
JP0003730646B2_000040t.gif以上の(1)~(3)の手続きにより、試験事後分布Qの計算を行うことができる。
[自由エネルギーの計算]
実施の形態2の自由エネルギーは以下の式で計算される。
JP0003730646B2_000041t.gif期待対数尤度LPおよびモデル複雑度項HJ,HZ,Hα,Hλ,Hβ,Hτは、以下のように計算される。
JP0003730646B2_000042t.gif実施の形態2でも、事後分布P(J|B)は、F(Q)を最大化するQJに対応し、対数周辺尤度log(P(B))が最大化されたF(Q)に対応する。
以上の手続きにより、ある深さの仮想曲面についての対数周辺尤度log(P(B))を求めることができる。
他の手続きは、実施の形態1と同様であるので、その説明は繰り返さない。
[階層事前分布のメタパラメータ]
階層事前分布の計算において、メタパラメータの選択について説明する。
他の観測方法による観測データが利用できないときには、実施の形態1と同様にメタパラメータの値は格子点によらずに共通に選ぶ。
すなわち、バーτ0,γβ0,γτ0は実施の形態1と同様に選ぶことができる。また、メタパラメータ バーα0n,γα0n,バーλ0n,γλ0nは、以下のように表わせる。
JP0003730646B2_000043t.gifここで、バーα0とγα0は実施の形態1と同様に選ぶことができる。
(実施の形態3)
実施の形態1および2においては他の観測方法による観測データが利用できないという前提に立って、脳内電流源の推定方法を説明してきた。
しかし、近年、脳の神経活動を観測するためのさまざまな手法が開発されており、同じ実験条件の下で複数の観測手段を用いた観測データが得られることも可能になってきた。そこで、実施の形態3においては、MRIやfMRIやPET等の他の観測方法による観測データが利用できる場合において、実施の形態1および2に対する変更点を説明する。
まず、同一人物に対するMRI観測データが得られる場合には、MRI画像から大脳皮質の位置やその他の脳部位の位置に関する情報が得られる。このような場合、実施の形態1および2において仮定された仮想曲面の形を大脳皮質の形や他の脳部位の形に合せることができる。
また、たとえば、ある実験条件下での神経活動が大脳皮質で起こっていることが確かなときには、深さ方向の探索を省略し、大脳皮質上にのみ格子点を置いて電流分布推定を行なうことができる。
fMRIやPETでは、神経活動の場所に関する情報を得ることができる。しかし、これらの情報の取扱には注意が必要である。なぜなら、fMRIでは血流量を、PETでは注入した放射性同位元素の量を測っており、これらは神経活動が高いところほど多くなると考えられるものの、神経活動による電流強度とは測定しているものが異なるからである。
fMRIやPETでの活動度と神経活動による電流強度はある程度の相関はあると考えられるが、その相関が線形的なものか、あるいは非線形な関係になっているかは明らかでない。また、MEGやEEGはミリ秒単位の時間分解能を持っているが、fMRIは秒単位の時間分解能であり、PETは数十秒単位の時間分解能しかない。すなわち、fMRIやPETで得られるのはMEGやEEGでの測定結果を長時間平均したものである。
Dale等は、文献:Dale A.M.et.Al.Neuron.26,pp.55-67(2000)において、MRIやfMRIやPETの観測結果をMEGやEEGの電流源推定に利用する方法を提案している。しかし、その方法は、fMRI等の活動情報を電流の分散情報として直接与えるものである。
これは、本発明におけるベイズ推定の立場からは、逆分散ハイパーパラメータαnの値を直接指定することに対応する。しかし、上述したようにfMRI等で得られる情報は電流そのものに対するものではなく、また時間的にも、MEGやEEGで測っている活動を長時間平均したものである。
これに対して、本発明において、fMRI等から得られる情報を電流の分散情報として直接指定するのではなく、階層事前分布のメタパラメータのレベルで指定する。メタパラメータγα0nとγλ0nは、階層事前分布により与える情報の信頼度を制御するメタパラメータである。すなわち、γα0nやγλ0nの値が小さいほど信頼性が低いことを表わし、γα0nやγλ0nの値が大きくなるほど信頼性が高いことを表わしている。
また、実施の形態2で説明したように、局在条件と連続条件を組合せた階層事前分布は、事前分布における電流の共分散行列としてβ-1(Λ-1+W・A-1・W′)を仮定していることに等しい。但し、AとΛはそれぞれ{αn|n=1,…,N}と{λn|n=1,…,N}を対角成分に持つ対角行列である。
すなわち、αnの値が小さいほど、電流の分散は大きくなり、λnの値が小さいほど電流の連続条件による共分散成分が大きくなる。メタパラメータ バーα0nとバーλ0nは、MEG(あるいはEEG)による観測データがないときのαnとλnの期待値を表わしている。
fMRIや他の観測方法により電流格子点nにおける神経活動が高いことが予想される場合は、メタパラメータ バーαUnとバーλUnの値を小さくして、この電流格子点における電流が大きい値をとりやすいようにバイアス情報を入れることができる。
このように、階層事前分布のメタパラメータのレベルで測っている量や時間スケールの異なる他の観測方法による情報を入れることで、曖昧な情報として他の観測手段の情報を取り入れることができる。
(シミュレーション結果)
以下、実施の形態2および3の手続によるシミュレーション結果について示す。以下では、ある人物のMRI画像から得られた大脳皮質の後頭部にある視覚野を含む領域を切り出したものを仮想曲面とし、深さ方向の探索は省略している。
シミュレーションデータとして、視覚野のV1,V2,V5の領域にそれぞれ電流強度が(1.0,1.0,0.5)(単位は任意)の電流を流した場合と、(1.0,0.0,0.5)の電流を流した場合について、これらの電流がMEGセンサ上に作る磁場を計算し、これにS/N比が0.1のガウスノイズを加えて観測データとした。
図28および図29は、局在条件のみを用いてベイズ推定した結果を示す。図中の各三角形の頂点は、この仮想曲面上に用意した格子点を表わし、この各々の格子点上に電流双極子が割当てられている。各格子点の間隔は平均約3mmであり空間分解能は高い。また、この仮想曲面は大脳皮質の一部を切り出したもので複雑な三次元形状をしている。図は、この三次元形状を二次元投影したものである。図中で、V1,V2,V5と示されているのは、それぞれ視覚野のV1,V2,V5領域を示している。図中で太い黒丸で示されている点は、推定電流強度が高い点を表わす。
図30および図31は、実施の形態2で説明した局在条件と連続条件を組合せた階層事前分布を用いてベイズ推定を行なった結果を示す。この推定結果は、ほぼ真の値に近い電流分布の位置と広がりを再現している。一方、図28と図29に示されているように、局在条件を用いてベイズ推定を行なった場合、電流の位置は正しく推定されているが、電流の広がりが過小評価されていることがわかる。このシミュレーションから空間分解能が高い場合に、連続条件を入れる有効性がわかる。
次に、電流分布はそのままにして、ノイズのS/N比を0.2にしてシミュレーションを行なった。この場合、局在条件と連続条件を組合せた階層事前分布を用いても、強いノイズの影響で電流分布の広がりが正確には再現されなかった。
続けて、このような劣悪の条件下で他の観測手段による観測データの情報が利用できる場合についてシミュレーションを行なった。シミュレーションでは、fMRIによる観測データからV1とV2とV5の領域の活動が高いという情報が得られたと仮定した。fMRIの活動領域と電流活動が必ずしも一致せず、またfMRIはMEG観測データの長時間平均に対応していることを考慮して、fMRIで高い活動を示している領域は、V1とV2とV5の各々の領域を半径20mm程度広げた近傍領域であるとした。また、V1,V2,V5の電流強度がそれぞれ(1.0,1.0,0.5)と(1.0,0.0,0.5)の両方の場合に同じfMRI情報を用いた。
以上の設定により、fMRIの活動領域と電流強度は必ずしも一致しないが、ある程度の相関はあるという状況をシミュレーションすることになる。このような設定の下で、fMRIの活性化領域に対応する格子点のメタパラメータをバーα0n=バーλ0n=2×Tr(G′・G)/ティルダNとし、そうでない格子点のメタパラメータをバーα0n=バーλ0n=50000×Tr(G′・G)/ティルダNとした。なお、変数の前に「ティルダ」との語句をつけた場合、その変数の頭部に「~」がつけえられていることを示す。
また、γα0nおよびγλ0nは全点共通にγα0n=γλ0n=0.1とした。
このときの推定結果を示したのが、図32と図33である。この図から、ノイズが非常に大きい場合でも他の観測手段の情報を加えることにより、正確な推定を行なえることがわかる。
この発明を詳細に説明し示してきたが、これは例示のためのみであって、限定となってはならず、発明の精神と範囲は添付の請求の範囲によってのみ限定されることが明らかに理解されるであろう。
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12(a)】
11
【図12(b)】
12
【図12(c)】
13
【図13(d)】
14
【図13(e)】
15
【図13(f)】
16
【図14】
17
【図15(a)】
18
【図15(b)】
19
【図15(c)】
20
【図16(d)】
21
【図16(e)】
22
【図16(f)】
23
【図17】
24
【図18】
25
【図19(a)】
26
【図19(b)】
27
【図19(c)】
28
【図20(d)】
29
【図20(e)】
30
【図20(f)】
31
【図21】
32
【図22(a)】
33
【図22(b)】
34
【図22(c)】
35
【図23(d)】
36
【図23(e)】
37
【図23(f)】
38
【図24】
39
【図25(a)】
40
【図25(b)】
41
【図25(c)】
42
【図26(d)】
43
【図26(e)】
44
【図26(f)】
45
【図27】
46
【図28】
47
【図29】
48
【図30】
49
【図31】
50
【図32】
51
【図33】
52