以下、本発明に係る脳活動情報出力装置等の実施形態について図面を参照して説明する。なお、実施の形態において同じ符号を付した構成要素は同様の動作を行うので、再度の説明を省略する場合がある。
(実施の形態1)
本実施の形態において、脳活動とアーチファクト源の事前情報を用いて、アーチファクトの影響を取り除き、精度の高い脳活動の情報を出力できる脳活動情報出力システムについて説明する。特に、アーチファクト源の事前情報を適切な事前情報にすることにより、極めて精度高く、脳活動の情報を出力できる脳活動情報出力システムについて説明する。
また、本実施の形態において、電流レベルで脳活動の情報を推定する脳活動情報出力システムについて説明する。また、本実施の形態において、アーチファクト源が複数存在する場合について説明する。また、本実施の形態において、脳事前情報とアーチファクト事前情報との信号の大きさが異なることを想定している脳活動情報出力システムについて説明する。また、アーチファクト源は、生体である眼球、心臓、筋電などであったり、TVの電波などの外部のものであったりする。さらに、本実施の形態において、脳活動信号をチェックし、精度が一定以上であることを確認後、脳信号を出力について説明する。なお、アーチファクト事前情報とはノイズ事前情報と同意義である。
図1は、本実施の形態における脳活動情報出力システムのブロック図である。脳活動情報出力システムは、脳信号計測装置11、脳活動情報出力装置12を具備する。
脳信号計測装置11は、脳信号を計測する装置である。脳信号計測装置11は、例えば、脳磁計測装置である。脳磁計測装置は、脳活動に関する事象を観測する観測手段としてのブレインキャップ1101を含んでも良い。ブレインキャップ1101は、例えば、208チャンネルセンサを有する。ただし、ブレインキャップ1101のチャンネル数は問わない。なお、脳信号とは、磁場でも良いし、電流等でも良い。
脳活動情報出力装置12は、計測脳信号情報格納部121、脳事前情報格納部122、アーチファクト事前情報格納部123、脳活動情報取得部124、アーチファクト情報取得部125、出力脳活動情報取得部126、ノイズ除去精度情報取得部127、出力部128、アーチファクト事前情報決定部129を具備する。
アーチファクト事前情報決定部129は、アーチファクト信号強度情報推定手段1291、fMRI活動強度格納手段1292、脳活動信号強度推定手段1293とを具備する。
脳活動情報出力装置12は、脳信号計測装置11が取得した脳の信号に関する情報である計測脳信号情報から、アーチファクト除去を行い、精度の高い脳活動情報を出力する装置である。
計測脳信号情報格納部121は、計測された脳の信号に関する情報である計測脳信号情報を格納し得る。計測脳信号情報は、脳磁計測装置(脳信号計測装置11の一例)の測定結果でも良いし、当該計測結果から、低周波のドリフト除去を行った情報でも良い。また、計測脳信号情報は、MEGデータでも良いし、電流の情報でも良い。また、通常、計測脳信号情報は、複数地点(例えば、208チャンネル)の情報であり、計測脳信号情報格納部121には、複数の脳信号情報が格納されている。
計測脳信号情報は、例えば、被験者の脳神経活動の電気的変化に伴って、被験者の頭部周辺で観測された微弱な磁場である。例えば、神経電流を位置と方向が固定されたN個の電流ダイポールでモデル化し、その電流モーメントの強度をN次元ベクトルJで表すM個のMEGセンサ上での観測磁場をM次元ベクトルBで表すと、次式が成り立つ。
このM次元ベクトルBが、計測脳信号情報の例である。また、上記のN次元ベクトルJを計測脳信号情報の例としても良い。つまり、脳表面近くのM点での測定した磁場をB=(B1,...,BM)T、脳表面上とアーチファクト源に仮定したI点のダイポールに流れる電流分布をJ={Ji|i=1,...,N}と表している。ここで、Gはリードフィールド行列と呼ばれ、単位電流モーメントによって発生する磁場を表す行列である。脳を囲む頭蓋内部が均一な媒質の球であると仮定すると、リードフィールドの成分はSarvasの式によって解析的に計算できる。MEG信号を生成する主要な脳活動成分は錐体細胞に流れる電流であるという知見に基づき、皮質表面上に垂直な向きを持つ電流ダイポールを仮定する、とする。なお、Sarvasの式は、公知の式であるので、詳細な説明は省略する。
計測脳信号情報格納部121は、不揮発性の記録媒体が好適であるが、揮発性の記録媒体でも実現可能である。計測脳信号情報格納部121に計測脳信号情報が記憶される過程は問わない。例えば、記録媒体を介して計測脳信号情報が計測脳信号情報格納部121で記憶されるようになってもよく、脳信号計測装置11から通信回線等を介して送信された計測脳信号情報が計測脳信号情報格納部121で記憶されるようになってもよく、あるいは、入力デバイスを介して入力された計測脳信号情報が計測脳信号情報格納部121で記憶されるようになってもよい。
脳事前情報格納部122は、脳活動に関する事前の情報である脳事前情報を格納し得る。脳事前情報は、例えば、脳の電流に対する事前分布であり、事前知識に基づいて定められる情報である。脳事前情報は、例えば、後述するP0(J)である。脳事前情報は、例えば、階層ベイズ推定を用いて、図示しない導出手段により導出されたものである。階層ベイズ推定では、次のようなダイポールに対する階層事前分布を取得する。また、脳事前情報は、信号の強度に関する情報である強度情報を含む。さらに、脳事前情報は、fMRI強度情報を電流分散とし、当該電流分布に対する事前知識を表す事前分布であっても良い。
以下、階層事前分布の取得方法を説明する前に、階層ベイズ推定を用いて、ダイポール位置での電流を推定する方法について説明する。
つまり、脳表面近くのM点での測定した磁場をB=(B1,...,BM)T、脳表面上とアーチファクト源に仮定したI点のダイポールに流れる電流分布をJ={Ji|i=1,...,N}と表す、とする。また、階層ベイズ推定において、電流分布Jが与えられたときに磁場Bが観測される確率P(B|J)が、数式2の正規分布に従うと仮定する。
数式2において、βは観測ノイズ分散の逆数である。また、磁場Bを観測したときの電流Jの事後分布(P(J|B))は次の数式3で表される。
数式3において、P(B)は周辺尤度である。P0(J)は電流に対する事前分布であり、所定の事前知識に基づいて定められる。階層ベイズ推定では、以下の数式4から数式6のようなダイポールに対する階層事前分布が導入される。
数式4において、Aはその対角成分がα={αi|i=1,...,N}である対角行列であり、α?1 iがJiの事前電流分散に対応する。βは、MEGデータの観測ノイズ分散の逆数であり、Γ(γ)はガンマ関数である。
ここで、事前電流分散を表すパラメータを、以下の数式7に示すものとする。
階層事前分布のパラメータであるν0iバー(バーは、νの上の−を有する文字を意味する)とγ0iはハイパーパラメータと呼ばれている。これらの値はfMRIの活動強度情報(t値)、倍率のパラメータm0、信頼度のパラメータγ0・γ0iから定義されている。ν0iバーとγ0iはfMRIのt値に応じて設定され、背景電流分散をν0,baseとすると、ν0iバーは、数式8のように表される。なお、ハイパーパラメータは、信号の強度に関する情報である強度情報の例である。
ここで、tiハット(ハットは、tの上に^が存在することを示す)は最大値が1となるようにfMRIのt値を正規化した値である。したがって、ν0iバーの取り得る値の範囲は、ν0,baseからm0×ν0,baseである。
背景電流分散ν0,baseと観測ノイズ分散ρbaseは,タスク開始前のMEG信号から最小ノルム推定を用いて推定された。事前分散は、数式5から,ρ×νiと表される。
νiバーが大きいと、電流の分散が大きくなる確率が高くなる(i番目のダイポールに電流が流れやすい、という事前情報を与えることを意味する)。ハイパーパラメータγ0iは階層事前分布の信頼度を制御する。γ0iが小さいほど信頼度が低く、大きいほど信頼度が高いことを表す。
脳事前情報は、例えば、図示しない脳事前情報取得手段により、fMRI活動強度情報を数式8に代入し,m0=100、γ0i=100として事前分散の推定を行い、取得された情報である。なお、脳事前情報は、皮質上に配置したダイポールに対して、事前分散の推定を行い、取得された情報である。
脳事前情報格納部122は、不揮発性の記録媒体が好適であるが、揮発性の記録媒体でも実現可能である。脳事前情報格納部122に脳事前情報が記憶される過程は問わない。例えば、記録媒体を介して脳事前情報が脳事前情報格納部122で記憶されるようになってもよく、通信回線等を介して送信された脳事前情報が脳事前情報格納部122で記憶されるようになってもよく、あるいは、入力デバイスを介して入力された脳事前情報が脳事前情報格納部122で記憶されるようになってもよい。
アーチファクト事前情報格納部123は、アーチファクト源に関する事前の情報である1以上のアーチファクト事前情報を格納し得る。アーチファクト事前情報格納部123は、複数のアーチファクト事前情報を格納していることは好適である。アーチファクト源とは、例えば、生体の一部である眼球の位置2箇所、頸動脈付近の左右1箇所ずつ、心臓の位置などである。眼球は、定常的に電気的な極性があり、角膜側を正電位、網膜側が負電位となっている。このため、眼球は角膜を正極、網膜側を負極としたダイポールだと考えることができる。MRIの構造画像から眼球の中心位置を決め、ダイポールを配置することは好適である。心臓の位置は、(x,y,z)=(?5,0,?23)[cm]の位置程度が好適である。また、頸動脈の位置は、MRIの構造画像を参考に決定することは好適である。また、アーチファクト源とは、例えば、生体以外の外部のアーチファクト源であるテレビ、輸送物(自動車、トラック、電車等)の通過などでも良い。つまり、アーチファクト事前情報は、被験者の生体に含まれるアーチファクト源に関する事前の情報である生体アーチファクト事前情報を含む。また、アーチファクト事前情報は、被験者の生体に含まない外部のアーチファクト源に関する事前の情報である外部アーチファクト事前情報を含んでも良い。アーチファクト事前情報は、脳事前情報と同様に、例えば、階層ベイズ推定を用いて得られる情報であり、数式4から数式6のようなダイポールに対する階層事前分布である。ただし、アーチファクト事前情報と脳事前情報とは、電流分散の大きさ(信号の強度である強度情報)が異なる。例えば、各種のアーチファクト事前情報は、数式9、10のように表される。また、アーチファクト事前情報は、強度情報を含む。また、アーチファクト事前情報は、後述するアーチファクト事前情報決定部129が取得した情報である。また、脳事前情報の強度情報と、1以上の各アーチファクト事前情報の各強度情報は、異なることは好適である。さらに、各アーチファクト事前情報は、fMRI強度情報と比較して大きな電流分散を用い、かつ当該電流分布に対する事前知識を表す事前分布であることは好適である。
ここで,νeyeバー(バーとは、νの上の−を有することを示す)、νcardiacバー、νmasticatoryバーはそれぞれ眼球中心、心臓位置、咀嚼筋に配置したダイポールに関するハイパーパラメータν0iバーである。また、γeye、γcardiac、γmasticatoryは、ハイパーパラメータγ0iに対して同様の値を表す。
なお、νbrainバー:νeyeバー:νcardiacバー:νmasticatoryバーは、例えば、「1:50:5000:100」などの簡易な比でも良い。
また、これらの事前電流分散は、例えば、図示しないアーチファクト事前情報取得手段が、観測磁場と復元磁場のあてはまりの度合(Goodness−of−fit,GoF)である第一の指標、fMRIの活動強度情報と平均電流分散の相関係数である第二の指標、の2種類の指標をもとに決定したものでも良い。
具体的には、GoFは数式11により定義される量とした。
数式11において、チャンネルjでの計測磁場をBobs(j)、チャンネルjにおける復元磁場をBest(j)とする。Mはセンサ数とする。推定電流から計算される磁場(復元磁場)と観測磁場が一致しているならばGoFは大きくなる(100に近づく)。推定結果と計測結果が一致するようなハイパーパラメータを選択するための判断材料としてGoFを用いた。第二の指標については、fMRIの活動が期待される位置では、脳内での神経活動も大きくなることが期待される。そのため、fMRIの活動強度情報と平均電流分散には対応関係があると考えられる(ふたつの相関係数が最も大きくなることが期待される)。第一の指標と第二の指標の条件をできるだけ満たすようなハイパーパラメータの組を選択することでアーチファクト源の電流強度に関する粗い情報をあらかじめ設定し、アーチファクト事前情報取得手段は、ダイポールに関する事前分散(アーチファクト事前情報)を推定し、取得するものとする。
アーチファクト事前情報格納部123は、不揮発性の記録媒体が好適であるが、揮発性の記録媒体でも実現可能である。アーチファクト事前情報格納部123にアーチファクト事前情報が記憶される過程は問わない。例えば、記録媒体を介してアーチファクト事前情報がアーチファクト事前情報格納部123で記憶されるようになってもよく、通信回線等を介して送信されたアーチファクト事前情報がアーチファクト事前情報格納部123で記憶されるようになってもよく、あるいは、入力デバイスを介して入力されたアーチファクト事前情報がアーチファクト事前情報格納部123で記憶されるようになってもよい。
脳活動情報取得部124は、計測脳信号情報格納部121に格納されている計測脳信号情報と、脳事前情報格納部122に格納されている脳事前情報とを用いて、脳活動の情報である脳活動情報を推定し、取得する。脳活動情報取得部124は、複数のダイポール位置での電流についての情報である脳活動情報を、計測脳信号情報(例えば、B)と脳事前情報(例えば、P0(J))とを用いた階層ベイズ推定を用いて推定し、取得することは好適である。例えば、脳活動情報取得部124は、上述した数式3の情報を格納しており、当該数式3を用いて、電流Jの事後分布(P(J|B))を取得する。電流Jの事後分布(P(J|B))は、脳活動情報である。脳活動情報取得部124は、通常、MPUやメモリ等から実現され得る。脳活動情報取得部124の処理手順は、通常、ソフトウェアで実現され、当該ソフトウェアはROM等の記録媒体に記録されている。但し、ハードウェア(専用回路)で実現しても良い。
アーチファクト情報取得部125は、計測脳信号情報と、1以上のアーチファクト事前情報とを用いて、アーチファクト源の情報である1以上のアーチファクト情報を推定し、取得する。アーチファクト情報取得部125は、複数のアーチファクト源ごとに、アーチファクト情報を推定し、複数のアーチファクト情報を取得することは好適である。アーチファクト情報取得部125は、アーチファクト源のダイポール位置での電流についての情報である1以上のアーチファクト情報を、計測脳信号情報と、1以上のアーチファクト事前情報とを用いた階層ベイズ推定を用いて推定し、取得することは好適である。アーチファクト情報取得部125は、事前情報が異なるが、脳活動情報取得部124と同様の動作を行う。つまり、アーチファクト情報取得部125は、数式3を用いて、各アーチファクト源に対応する電流Jの事後分布(P(J|B))を取得する。アーチファクト情報取得部125は、通常、MPUやメモリ等から実現され得る。アーチファクト情報取得部125の処理手順は、通常、ソフトウェアで実現され、当該ソフトウェアはROM等の記録媒体に記録されている。但し、ハードウェア(専用回路)で実現しても良い。
出力脳活動情報取得部126は、計測脳信号情報と脳活動情報と1以上のアーチファクト情報とを用いて、計測脳信号情報からアーチファクトの影響を取り除く処理を行い、出力される脳活動の情報である出力脳活動情報を取得する。つまり、アーチファクトの影響を取り除くために、脳皮質のダイポールとアーチファクト源のダイポールの推定を同時に行うこととなる。なお、脳皮質のダイポールの推定とは、脳活動情報取得部124が脳活動情報を取得する処理である。アーチファクト源のダイポールの推定とは、アーチファクト情報取得部125が、アーチファクト源の情報である1以上のアーチファクト情報を取得する処理である。
また、出力脳活動情報取得部126は、計測脳信号情報と脳活動情報と複数のアーチファクト情報とを用いて、計測脳信号情報から複数のアーチファクトの影響を取り除く処理を行い、出力される脳活動の情報である出力脳活動情報を取得することは好適である。
さらに、出力脳活動情報取得部126は、計測脳信号情報から、脳活動情報のうちの脳の活動が無いと考えられる位置に配置されたダイポールに対応する脳活動情報とアーチファクト情報とを除く処理を行い、出力脳活動情報を取得することは好適である。なお、脳の活動が無いと考えられる位置に配置されたダイポールの情報は、予め設定されていても良いし、ユーザが手入力により与えるなどしても良い。
以下、出力脳活動情報取得部126等の具体的な動作例について説明する。今、脳皮質とアーチファクト源のダイポールをそれぞれJbrain,Jnoise,リードフィールド行列をGbrain,Gnoise,ダイポールの事前分散をA?1 brain,A?1 noiseとすると、以下の式12から14を用いて、脳活動情報取得部124、およびアーチファクト情報取得部125は、それぞれ脳皮質とアーチファクト源のダイポールJを推定する。
つまり、脳活動情報取得部124とアーチファクト情報取得部125により、階層ベイズ推定法でアーチファクト源と脳皮質のダイポールを同時に推定した後に、出力脳活動情報取得部126は、fMRIの活動がない位置に配置されたダイポールのみが作る磁場を、数式15にしたがって観測磁場Bから取り除き、Bprunedを算出し、メモリ上に配置する。なお、観測磁場Bは、計測脳信号情報の一例であり、Bprunedは、出力脳活動情報の一例である。
出力脳活動情報取得部126は、通常、MPUやメモリ等から実現され得る。出力脳活動情報取得部126の処理手順は、通常、ソフトウェアで実現され、当該ソフトウェアはROM等の記録媒体に記録されている。但し、ハードウェア(専用回路)で実現しても良い。
ノイズ除去精度情報取得部127は、出力脳活動情報取得部126が取得した出力脳活動情報と、計測脳信号情報とを用いて、当該出力脳活動情報におけるノイズ除去の精度に関する情報であるノイズ除去精度情報を取得する。ノイズ除去精度情報取得部127の具体的な処理の例は、後述する。ノイズ除去精度情報取得部127は、通常、MPUやメモリ等から実現され得る。ノイズ除去精度情報取得部127の処理手順は、通常、ソフトウェアで実現され、当該ソフトウェアはROM等の記録媒体に記録されている。但し、ハードウェア(専用回路)で実現しても良い。
出力部128は、出力脳活動情報を出力する。出力部128は、ノイズ除去精度情報が、予め格納されている閾値と比較して高い精度である場合のみ、出力脳活動情報を出力しても良い。ここで、出力とは、ディスプレイへの表示、プロジェクターを用いた投影、プリンタへの印字、音出力、外部の装置への送信、記録媒体への蓄積、他の処理装置や他のプログラム等への処理結果の引渡し等を含む概念である。出力部128は、ディスプレイやスピーカー等の出力デバイスを含むと考えても含まないと考えても良い。出力部128は、出力デバイスのドライバーソフトまたは、出力デバイスのドライバーソフトと出力デバイス等で実現され得る。
アーチファクト事前情報決定部129は、アーチファクト源の信号強度に関する事前分布を取得し、かつ、当該事前分布を用いて当該アーチファクト源の信号強度に関する事後分布を推定し、事前分布と事後分布との差を算出し、当該差が予め決められた範囲以内の小さな差である信号強度に関する情報を有する1以上のアーチファクト事前情報を取得し、アーチファクト事前情報格納部123に蓄積する。
なお、アーチファクト事前情報決定部129は、アーチファクト信号強度情報推定手段1291と脳活動信号強度推定手段1293の処理により、アーチファクト事前情報を取得し、アーチファクト事前情報格納部123に蓄積することはさらに好適である。つまり、アーチファクト事前情報決定部129は、以下の2つの基準を満たすアーチファクト事前情報を取得することは好適である。第一の基準は、統計モデルの妥当性であり、第二の基準は、fMRIと推定電流の信号強度情報との相関である。つまり、アーチファクト事前情報決定部129は、アーチファクト信号強度情報推定手段1291により、統計モデルが妥当なアーチファクト事前情報を取得する。また、アーチファクト事前情報決定部129は、脳活動信号強度推定手段1293により、fMRIと推定電流の信号強度情報との相関性が高いアーチファクト事前情報を取得する。
アーチファクト信号強度情報推定手段1291は、アーチファクト源の信号強度に関する事前分布を取得し、かつ、当該事前分布を用いて当該アーチファクト源の信号強度に関する事後分布を推定し、前記事前分布と前記事後分布との差を算出し、当該差が予め決められた範囲以内の小さな差である信号強度に関する情報を有する1以上のアーチファクト事前情報の候補を取得する。アーチファクト源の電流強度に関する情報(事前電流分散(上記のν0iバー)が妥当ならば、それに基づいて推定される電流強度(事後電流分散νiバー)も同じになる。そのため、アーチファクト信号強度情報推定手段1291は、アーチファクト源を表すそれぞれのダイポールの事前分散ν0iバー×ρbaseと信頼度のパラメータγ0iを変化させて電流推定を行い、以下の数式16に基づき、指標を計算する。そして、アーチファクト信号強度情報推定手段1291は、1以上のアーチファクト事前情報の候補を取得する。
数式16において、iはアーチファクト源のダイポール番号、Oi,priは事前分散のオーダー、Oi,postは事後分散のオーダー、Oi,diffは事前分散と事後分散のオーダーの差、ν0iバーは事前電流分散、νiバーは事後電流分散、ρbaseは観測ノイズ分散を表す。
アーチファクト信号強度情報推定手段1291は、例えば、事前分散とγ0iはそれぞれ1〜1011nAm2,1〜106の範囲で変化させ(事前分散12種類,γ0i13種類,合計156通り)、階層ベイズ法の計算を行う。そして、アーチファクト信号強度情報推定手段1291は、例えば、Oi,diffがゼロに近い値のハイパーパラメータの組合せを計算結果から選択し、アーチファクト事前情報の候補として取得する。
アーチファクト事前情報決定部129は、通常、MPUやメモリ等から実現され得る。アーチファクト事前情報決定部129の処理手順は、通常、ソフトウェアで実現され、当該ソフトウェアはROM等の記録媒体に記録されている。但し、ハードウェア(専用回路)で実現しても良い。
fMRI活動強度格納手段1292は、fMRIの活動強度を格納し得る。fMRI活動強度格納手段1292は、脳事前情報格納部122と同様であっても良い。かかる場合、脳事前情報は、fMRI強度情報を電流分散とし、当該電流分布に対する事前知識を表す事前分布である。fMRI活動強度格納手段1292は、不揮発性の記録媒体が好適であるが揮発性の記録媒体でも良い。
脳活動信号強度推定手段1293は、アーチファクト信号強度情報推定手段1291が取得した1以上の各アーチファクト事前情報の候補が有する1以上の信号強度を取得し、当該1以上の各信号強度と前記fMRIの活動強度との相関に関する情報を算出し、当該相関に関する情報が最大の信号強度を有するアーチファクト事前情報を、アーチファクト信号強度情報推定手段1291が取得した1以上のアーチファクト事前情報の候補から選択し、当該選択したアーチファクト事前情報をアーチファクト事前情報格納部123に蓄積する
さらに具体的には、脳活動信号強度推定手段1293は、以下の数式17から数式19を用いて、アーチファクト事前情報を選択する。
数式17から数式19は、fMRIの活動がある位置では、脳内での神経活動も大きくなると期待されることを利用した数式である。数式17から数式19は、fMRIの活動強度情報fi(fMRI活動強度格納手段1292に格納されている)と、1以上の各アーチファクト事前情報の候補が有する1以上の電流強度Jn,iの相関指標を算出する式である。数式17から数式19において、iは皮質のダイポール番号、nは試行番号、Tはサンプル数、Nは試行数を表す。脳活動信号強度推定手段1293は、単一試行毎に推定電流強度とfMRI信号強度の相関係数corrnを計算する。そして、脳活動信号強度推定手段1293は、その平均値corrを算出し、当該平均値corrが最大となるアーチファクト事前情報を取得する。また、アーチファクト事前情報決定部129は、それぞれのダイポールに対して選んだハイパーパラメータを使い、数式16による処理(妥当な統計モデルを選択する処理)と数式17から数式19による処理(fMRIと推定電流の信号強度情報との相関性が高いものを選択する処理)を、ハイパーパラメータの値が変化しなくなるまで繰り返し、最終的なハイパーパラメータの組を決定する。そして、このハイパーパラメータの組が、例えば、アーチファクト事前情報となる。また、数式18において、Iはダイポール番号の最大値(アーチファクト源の総数)である。
脳活動信号強度推定手段1293は、通常、MPUやメモリ等から実現され得る。脳活動信号強度推定手段1293の処理手順は、通常、ソフトウェアで実現され、当該ソフトウェアはROM等の記録媒体に記録されている。但し、ハードウェア(専用回路)で実現しても良い。
次に、脳活動情報出力システムを構成する脳活動情報出力装置12の動作について、図2のフローチャートを用いて説明する。
(ステップS201)脳活動情報取得部124は、計測脳信号情報格納部121から計測脳信号情報を読み出す。
(ステップS202)脳活動情報取得部124は、脳事前情報格納部122から脳事前情報を読み出す。
(ステップS203)脳活動情報取得部124は、予め格納している数式であり、脳活動情報等の推定のための数式を、メモリ上に読み出す。かかる数式は、例えば、階層ベイズ推定を実現するための式である。
(ステップS204)脳活動情報取得部124は、ステップS203で読み出した式に対して、ステップS201で読み出した計測脳信号情報と、ステップS202で読み出した脳事前情報を代入し、当該式を実行し、脳活動情報を取得し、メモリ上に配置する。
(ステップS205)アーチファクト情報取得部125は、カウンタiに1を代入する。
(ステップS206)アーチファクト情報取得部125は、アーチファクト事前情報格納部123に、i番目のアーチファクト事前情報が存在するか否かを判断する。i番目のアーチファクト事前情報が存在すればステップS207に行き、存在しなければステップS210に行く。
(ステップS207)アーチファクト情報取得部125は、アーチファクト事前情報格納部123から、i番目のアーチファクト事前情報を読み出す。
(ステップS208)アーチファクト情報取得部125は、ステップS203でメモリ上に読み出された数式に、ステップS201で読み出された計測脳信号情報と、ステップS207で読み出したアーチファクト事前情報を代入し、当該式を実行し、アーチファクト情報を取得し、メモリ上に配置する。
(ステップS209)アーチファクト情報取得部125は、カウンタiを1、インクリメントし、ステップS206に戻る。
(ステップS210)出力脳活動情報取得部126は、ステップS204で取得された脳活動情報から、fMRI活動のない位置のダイポールに対する情報を取り除く。
(ステップS211)出力脳活動情報取得部126は、ステップS210で残った情報から、ステップS208で取得された1以上のアーチファクト情報を取り除く。そして、出力脳活動情報取得部126は、出力脳活動情報を得る。
(ステップS212)出力部128は、ステップS211で取得された出力脳活動情報を出力し、処理を終了する。
なお、図2のフローチャートにおいて、ステップS210、ステップS211における処理の順序、タイミングは問わない。ステップS210、ステップS211における処理は、例えば、上述した数式15を用いて行われる。
また、図2のフローチャートにおいて、ノイズ除去精度情報取得部127の処理は行わなかった。しかし、図2のフローチャートにおいて、出力部128は、ノイズ除去精度情報が、予め格納されている閾値と比較して高い精度である場合のみ、出力脳活動情報を出力しても良い。
また、以下、脳活動情報出力装置12が、適切なアーチファクト事前情報をアーチファクト事前情報格納部123に蓄積する動作について、図3のフローチャートを用いて説明する。
(ステップS301)アーチファクト事前情報決定部129は、カウンタiに1を代入する。
(ステップS302)アーチファクト信号強度情報推定手段1291は、i番目のアーチファクト源が存在するか否かを判断する。i番目のアーチファクト源が存在すればステップS303に行き、存在しなければステップS312に行く。
(ステップS303)アーチファクト信号強度情報推定手段1291は、i番目のアーチファクト源の信号強度に関する事前分布を、アーチファクト事前情報格納部123から取得する。ここで、信号強度に関する事前分布は、例えば、電流強度に関する情報(事前電流分散)である。
(ステップS304)アーチファクト信号強度情報推定手段1291は、ステップS303で取得した事前分布を用いて、i番目のアーチファクト源の信号強度に関する事後分布を、1以上取得する。ここで、事後分布は、例えば、i番目のアーチファクト源の事前電流分散に基づいて推定される電流強度(事後電流分散)である。なお、この処理を、事後分布取得処理という。
(ステップS305)アーチファクト信号強度情報推定手段1291は、カウンタjに1を代入する。
(ステップS306)アーチファクト信号強度情報推定手段1291は、ステップS304で取得された事後分布のうち、j番目の事後分布が存在するか否かを判断する。j番目の事後分布が存在すればステップS307に行き、存在しなければステップS311に行く。
(ステップS307)アーチファクト信号強度情報推定手段1291は、i番目のアーチファクト源の信号強度に関する事前分布と、j番目の信号強度に関する事後分布との差を算出する。
(ステップS308)アーチファクト信号強度情報推定手段1291は、ステップS307で算出した差が、予め決められた閾値(例えば、「0」や「0.1」など)以内であるか否かを判断する。予め決められた閾値以内であればステップS309に行き、予め決められた閾値より大きければステップS310に行く。
(ステップS309)アーチファクト信号強度情報推定手段1291は、i番目のアーチファクト源の信号強度に関する事前分布を、バッファに一時格納する。なお、ここで、アーチファクト信号強度情報推定手段1291は、i番目のアーチファクト源の信号強度に関する事前分布と、信頼度のパラメータ(γ0i)とを、バッファに一時格納することは好適である。また、ここで、ループを抜けて、ステップS311に行っても良い。また、ここでバッファに格納されたアーチファクト源の信号強度に関する事前分布は、最終的にアーチファクト事前情報格納部123に格納されるアーチファクト事前情報(または、その一部)の候補となる情報である。
(ステップS310)アーチファクト信号強度情報推定手段1291は、カウンタjを1、インクリメントする。ステップS306に戻る。
(ステップS311)アーチファクト事前情報決定部129は、カウンタiを1、インクリメントする。ステップS302に戻る。
(ステップS312)脳活動信号強度推定手段1293は、カウンタiに1を代入する。
(ステップS313)脳活動信号強度推定手段1293は、ステップS309でバッファに格納されたアーチファクト事前情報の候補に、i番目の候補が存在するか否かを判断する。i番目の候補が存在すればステップS314に行き、存在しなければステップS322に行く。
(ステップS314)脳活動信号強度推定手段1293は、カウンタjに1を代入する。
(ステップS315)脳活動信号強度推定手段1293は、j番目の試行が存在するか否かを判断する。j番目の試行が存在すればステップS316に行き、存在しなければステップS320に行く。
(ステップS316)脳活動信号強度推定手段1293は、i番目のアーチファクト事前情報の、j番目の試行の信号強度(例えば、電流強度)を取得する。ここで、信号強度は、数式17により算出される。
(ステップS317)脳活動信号強度推定手段1293は、i番目のダイポールのfMRI信号強度を取得する。なお、本ステップはjのループの中にあるが、jのループの外に出しても良い。
(ステップS318)脳活動信号強度推定手段1293は、j番目の試行の相関係数(相関に関する情報の一例)を算出し、バッファに一時格納する。相関係数は、数式18により算出される。
(ステップS319)脳活動信号強度推定手段1293は、カウンタjを1、インクリメントする。ステップS315に戻る。
(ステップS320)脳活動信号強度推定手段1293は、ステップS318において、バッファに一時格納した相関係数の平均値を算出する(数式19参照)。そして、脳活動信号強度推定手段1293は、相関係数の平均値をバッファに一時格納する。
(ステップS321)脳活動信号強度推定手段1293は、カウンタiを1、インクリメントする。ステップS313に戻る。
(ステップS322)脳活動信号強度推定手段1293は、ステップS320でバッファに格納した相関係数の平均値のうちで、最大の相関係数の平均値に対応するアーチファクト事前情報(事前電流分布、または事前電流分布と信頼度のパラメータ等)を取得し、バッファに一時格納する。上位処理にリターンする。
また、以下、ステップS304の事後分布取得処理について、図4のフローチャートを用いて説明する。
(ステップS401)アーチファクト信号強度情報推定手段1291は、カウンタiに1を代入する。
(ステップS402)アーチファクト信号強度情報推定手段1291は、i番目の事前分布が存在するか否かを判断する。i番目の事前分布が存在すればステップS403に行き、存在しなければ上位処理にリターンする。なお、アーチファクト信号強度情報推定手段1291は、例えば、事前分散を、1〜1011nAm2の範囲で12種類(iは1から12まで)、変化させる。
(ステップS403)アーチファクト信号強度情報推定手段1291は、i番目の事前分布を取得する。
(ステップS404)アーチファクト信号強度情報推定手段1291は、カウンタjに1を代入する。
(ステップS405)アーチファクト信号強度情報推定手段1291は、j番目の信頼度が存在するか否かを判断する。j番目の信頼度が存在すればステップS406に行き、存在しなければステップS409に行く。なお、アーチファクト信号強度情報推定手段1291は、例えば、信頼度を、1〜106の範囲で13種類(jは1から13まで)、変化させる。
(ステップS406)アーチファクト信号強度情報推定手段1291は、j番目の信頼度を取得する。
(ステップS407)アーチファクト信号強度情報推定手段1291は、階層ベイズ法により、i番目の事前分布とj番目の信頼度とを用いて、事後分散(例えば、事後電流分散)を取得する。
(ステップS408)アーチファクト信号強度情報推定手段1291は、カウンタjを1、インクリメントする。ステップS405に戻る。
(ステップS409)アーチファクト信号強度情報推定手段1291は、カウンタjを1、インクリメントする。ステップS402に戻る。
以下、本実施の形態における脳活動情報出力システムを用いた実験結果について説明する。被験者は、正常な視覚を有する健常男性3名である。被験者3名(X,Y,Z)が、MEG計測実験とfMRI計測実験の両方に参加した。実験に先だって、被験者には実験内容についての説明を行い、同意を得た上で実験を行った。
MEG計測実験とは、以下のような実験である。脳信号計測装置11は、全頭型の208チャンネルセンサを有する市販の脳磁計測装置である。そして、脳信号計測装置11により、1000HzでMEGデータが記録した。なお、MEGデータは、計測脳信号情報の一例であるとも考えられる。
そして、実験課題として、被験者は運動する視標を、固視点を注視したまま、心的に追従する課題(covert pursuit task)を行った。視標は、正弦波を基礎として4分の1周期毎に速度が0になるとその周波数が変更されるようにした。ここで、実験で被験者に提示した視標の速度波形を、図5に示す。また、視標運動の1試行あたり平均周波数は0.5Hzであり、赤色の点とした。具体的には、以下の(a)から(d)の流れで課題を遂行した。
(a)被験者は、準備が出来たら、タスク開始のボタンを押す。(b)固視点と視標が画面中央に現れる(1000±500ms間待つ)。(c)その後、視標が水平方向に動き始める。被験者は、固視点を見つめながら視標を心の中で追従するように指示された。視標は4000ms間動き続けた。(d)視標が消え、被験者は休憩する。そして、(a)から(d)の流れを1試行とした。指標の運動パターンは3種類あり、それぞれ30試行ずつランダムに提示した(3種類×30試行=合計90試行)。これを1セッションとして、被験者は合計で3セッション(1セッション90試行×3セッション=合計270試行)課題を遂行した。
次に、fMRI計測実験について説明する。fMRI計測には市販の1.5TのMRIスキャナを用いた。視覚刺激は眼前から370mm離れた位置に設置したスクリーンに提示した。被験者は、正弦波の動きをする指標を、目を動かさずに心的に追従するテストブロックと、静止する固視点を注視しつづけるコントロールブロックを、それぞれ15秒間ずつ繰り返した。テストブロックでは、MEG計測実験と同じ視標を用いた。
まず、最初に、脳信号計測装置11が計測し、取得したMEGデータは、リファレンスセンサの値を用いて、図示しない低周波ドリフト除去手段(脳信号計測装置11が有する)により、低周波のドリフトの影響が除去された。次に、図示しないデータ切出手段(脳信号計測装置11が有する)により、指標の提示時刻を基準にして−400〜+4000msのデータを切り出し、これを1試行のデータとした。データ切出手段により、切り出されたMEGデータの最初の400ms間は指標が動き出す前の状態であり(pretrigger期間)、図示しない補正手段(脳信号計測装置11が有する)により、この期間の平均値を用いて基線補正を行った。上記の前処理をしたMEGデータが、計測脳信号情報格納部121に格納された計測脳信号情報の例であり、上述した階層ベイズ推定により、後述する電流推定に用いられた。
また、解析の際、ひとつ以上のセンサのデータが、計測可能な値の範囲の95%以上の値をとった場合に、その試行データは除外した。また、実験中に被験者が眠ったセッションのデータも除外した。そのため、被験者A,B,Cから計測されたデータのうち、それぞれ169試行(全試行の62.59%)、251試行(92.96%)、267試行(98.89%)のデータが解析に用いられた.
fMRIデータの解析はSPM5(The Wellcome Department of Cognitive Neurology)を用いて行われた。静止する固視点のみが提示され、それを注視している条件と、固視点を注視したままその上側で運動する指標を心の中で追従している条件での脳活動の有意差検定をt検定で行った。
あらかじめ計測しておいた被験者のMRI構造画像をもとに、Brain Voyager(Brain Innovation)を用いて大脳表面モデルを作成した。その全脳モデル上に、その全脳モデル上にダイポールが配置された(図6の皮質上(cortex)の点)。皮質表面モデル上に配置された格子点の数は、被験者A,B,Cでそれぞれ、13,755点、13,160点、12,370点であった。
また、本実験では、主なアーチファクト源として、以下に示す理由により、眼球付近の位置2箇所、頸動脈付近の左右1箇所ずつ、心臓の位置を仮定し、それぞれの場所にダイポールを仮定した。
まず、眼球は、定常的に電気的な極性があり、角膜側を正電位、網膜側が負電位となっている。このため、眼球は角膜を正極、網膜側を負極としたダイポールだと考えることができる。また、MRIの構造画像から眼球の中心位置を決め、ダイポールを配置した。
また、心拍成分のR波をトリガ時刻として、計測したMEGデータの時間軸を揃え直し,加算平均波形を求めた。
そのデータを用いてムービング・ダイポール法で電流源推定を行った。前交連(anterior commissure,AC)を原点として,前交連と後交連(posterior commissure,PC)を結ぶ線をy軸とする座標系を考え(x軸は正中矢状面に垂直な左右方向,z軸はこれらに垂直な上下方向)、探索範囲は被験者の上半身付近(x=?20〜20cm,y=?20〜20cm,z=?50〜20cm)とした。電流源推定の結果、心臓付近と頸部付近に推定された。その結果に基づいて、ダイポールを心臓位置に1箇所と頸動脈付近に左右1箇所ずつ配置した(図6参照)。心臓の位置は、(x,y,z)=(?5,0,?23)[cm]の位置である。頸動脈の位置は、MRIの構造画像も参考に決定した。
計測脳信号情報格納部121に格納されているMEGデータから、上述した階層ベイズ推定によって電流推定を行う際、fMRI活動強度情報(t値)を階層ベイズ推定の事前知識として用いた。つまり、脳事前情報格納部122には、fMRI活動強度情報(t値)を含む脳事前情報が格納されている。脳事前情報として、fMRI活動強度の大きいダイポールには大きな事前分散を与え、逆に活動強度が小さければ小さい事前分散を与えた。また、皮質上に配置したダイポールに対するハイパーパラメータは、皮質に対する事前分散の倍率のパラメータ「m0=500」、各ダイポールに対する信頼度のパラメータ「γ0i=500)を用いた。
アーチファクトの信号強度は未知である。アーチファクト源の電流分散に対するハイパーパラメータを決定するため、上述したように、二つの基準を設けた。一つは(a)統
計モデルの妥当性であり、もう一つは(b)fMRIと推定電流の信号強度情報との相関である。
まず、第一の基準を満たすために、以下のようにした。アーチファクト源の電流強度に関する情報(事前電流分散ν0iバー)が妥当ならば、それに基づいて推定される電流強度(事後電流分散νiバー)も同じになる。アーチファクト源を表すそれぞれのダイポールの事前分散ν0iバー×ρbaseと信頼度のパラメータγ0iを変化させて電流推定を行い、上述した数式17から数式19を用いて、指標を計算した。
アーチファクト信号強度情報推定手段1291は、事前分散とγ0iはそれぞれ1〜1011nAm2,1〜106の範囲で変化させ(事前分散12種類,γ0i13種類,合計156通り)、階層ベイズ法の計算を行った。そして、アーチファクト信号強度情報推定手段1291は、例えば、Oi,diffがゼロに近い値のハイパーパラメータの組合せを計算結果から選択し、アーチファクト事前情報の候補として取得した。その結果、アーチファクト事前情報の候補を、図7に示す線(1)上の点の集合のように得た。
次に、脳活動信号強度推定手段1293は、以下のように動作した。つまり、脳活動信号強度推定手段1293は、fMRIの活動強度情報fiと電流強度Jn,iの相関指標を、数式17から数式19を用いて、算出した。これは、fMRIの活動がある位置では、脳内での神経活動も大きくなると期待されることを用いたものである。
そして、アーチファクト事前情報決定部129は、(a)統計モデルの妥当性、(b)fMRIと推定電流の信号強度情報との相関、の2つの基準を同時に満たすハイパーパラメータの組を選択した。具体的には、アーチファクト事前情報決定部129は、図8の黒点(1)を選んだ。
そして、その後、アーチファクト事前情報決定部129は、それぞれのダイポールに対して選んだハイパーパラメータを使い、数式17から数式19の計算をやり直した。アーチファクト事前情報決定部129は、これらの手順をハイパーパラメータの値が変化しなくなるまで繰り返した。そして、アーチファクト事前情報決定部129は、最終的なハイパーパラメータの組(アーチファクト事前情報)を取得した。繰り返した回数は、被験者A,B,Cでそれぞれ5回,3回,5回であった。つまり、アーチファクト事前情報決定部129は、1以上のアーチファクト事前情報の推定を、所定の条件に合致するまで繰り返すことは好適である。ここで、所定の条件とは、アーチファクト事前情報(例えば、ハイパーパラメータの値)が所定以内の変化(変化しない場合も含む)になることである。
また、アーチファクト源のハイパーパラメータを変化させることが電流強度に与える影響を調べるため、数式20を用いて計算した。
ここで、図9はハイパーパラメータの組ごとに、数式20を計算した値をマップとして表している。ハイパーパラメータを変化させることで、アーチファクト源の電流強度は変化する。そのことによって差し引く磁場の大きさも変化するため、ハイパーパラメータの選択がMEGアーチファクト除去の性能に大きな影響を及ぼす。
また、3種類の視標運動パターン遂行中のMEGデータのうち、1種類を用いて、アーチファクト源のハイパーパラメータ(アーチファクト事前情報)を決定した。
また、ハイパーパラメータは、推定電流とfMRI活動強度情報との相関が高くなるように決めているため、パラメータ推定に使ったデータセットを同じデータセットに適用して、アーチファクト除去の精度を推定するのは不適切である。そこで、パラメータの決定に用いなかった残り2種類のMEGデータを使って、アーチファクト除去の精度を評価した。
以上により、計測脳信号情報格納部121、脳事前情報格納部122、およびアーチファクト事前情報格納部123に、それぞれの情報が格納された。
次に、出力脳活動情報取得部126等は、数式12から数式15を用いて、上述したようなアーチファクト除去を行った。つまり、階層ベイズ推定法でアーチファクト源と大脳皮質の電流モーメントを同時に推定した後、出力脳活動情報取得部126は、数式15に従って、アーチファクトのみが作る磁場を観測磁場Bから取り除き、Bprunedを求めた。
以下、従来のアーチファクト除去方法について述べる。また、従来のアーチファクト除去方法を用いた場合の実験と、脳活動情報出力装置12によるアーチファクト除去方法を用いた場合の実験について述べる。
従来のアーチファクト除去方法とは、主成分分析(PCA)と独立成分分析(ICA)によるアーチファクト除去である。従来方法により、計測したMEGデータは、セッションと視標速度の条件毎にデータセットに分け、それぞれに対してPCAとICAを適用した。
まず、計測したMEGデータにPCAを適用し、その第一主成分が心拍波形に類似していることを目視で確認した上で、その成分を除去した。図10は、観測したMEGデータにPCAを適用したときの第一主成分の時間波形を表す図である。図11は、第一主成分を頭部マップ上に示した図である。
次に、計測した208chのMEGデータにInfomax ICAを適用し、208個の独立成分を抽出した。それらの中で、心拍波形に似た特徴や眼球付近のセンサでのみ変化がある独立成分を目視で確認した場合、それらをアーチファクトとみなした。図12、図13は、観測したMEGデータにICAを適用したときの独立成分の時間波形を表す。図14、図15は、各々図12、図13に対応し、独立成分を頭部マップ上に示した図である。図10から図15は、被験者Bの1セッション目の一つの視標運動遂行中のMEGデータにPCAを適用したときの結果である。図10、図12、図13では83番目のセンサの値を示した。
計測したMEGデータからアーチファクトを表す独立成分を差し引き、その値をアーチファクト除去したMEGデータとした。被験者A,B,Cでそれぞれ2個,2.00±0.50個,1.11±0.33個の独立成分が除去された。
ICAでアーチファクトとみなした独立成分のうち、心拍アーチファクトの影響が特に大きかった83番目のセンサ値からR波のピーク時刻が検出された。その時刻をトリガとして、観測磁場から−150〜+450msの時間区間が切り出された。被験者A,B,Cでそれぞれ795,876,1100サイクルの心拍が検出され、その加算平均波形を解析に用いた。なお、ICAで心拍アーチファクトとみなした成分は、R波のピーク時刻検出にのみ用い、加算平均波形の計算には観測磁場を用いた。なお、83番目のセンサの位置を、図16に示す(図16の(1)の点)。
次に、脳活動情報出力装置12によるアーチファクト事前情報の取得について述べる。図17は、心拍成分のR波をトリガとして観測磁場を切り出し、加算平均した値を示す図である。図18、図19、図20は、図17と同じR波のピーク時刻を基準に、心臓、頸動脈(左)、頸動脈(右)に配置したダイポールの推定電流を加算平均した波形である。図18、図19、図20において、各ダイポールのx,y,z軸成分のそれぞれの波形を示している。観測磁場の心拍成分は、R波とS波の区間でインパルス状の波形を示している。推定電流波形も、これらと同じ時間区間で同様の特徴を示した。また、Q波とT波の時間区間では、観測磁場は緩やかに変化し、推定電流も同様の変化を示した。このように観測磁場に見られる心拍成分の特徴は、推定電流の時間波形でも確認された。眼球には、角膜を正極、網膜側を負極とした定常的な電位がある(WR. Miles, "The steady polarity potential of the human eye," Proc. Natl. Acad. Sci. U. S. A., 25,1, pp.25-36, 1939. 参照)。被験者が点滅する光刺激を見ているときの頭部付近の磁場を計測し、その磁場から眼球位置に配置したダイポールの電流モーメントが推定されている。その値は約10nAmであったことを示している(T. Katila, R. Maniewski, T. Poutanen, and T.Varpula, "Magnetic fields produced by the human eye," J. Appl. Phys., 52,2, pp. 2565-2571, 1981. 参照)。本実験の観測磁場から推定された眼球位置での電流強度を、数式17を用いて試行毎に計算した結果,試行平均値は1.88−8.54nAmであり、先行研究と同じオーダーの値が推定された。
追従性眼球運動遂行中のMEGとEOGを計測した従来の実験がある(Y. Fujiwara, O. Yamashita, D. Kawawaki, K. Doya,M. Kawato, K. Toyama, and MA. Sato, "A ierarchical Bayesian method to resolve an inverse problem of MEG contaminated with eye movement artifacts," NeuroImage, 45, 2, pp. 393-409, 2009.参照)。従来の実験において、動かす眼球の振幅角は4度,振動周波数は0.7Hzであり,その実験中のEOGデータから眼球アーチファクトの順モデルを求め、観測磁場を実現する電流モーメントを推定した。その推定電流の試行平均値は21.7-36.3nAmであった。
脳活動情報出力装置12による本実験(本実験タスク)の場合は、眼球を動かさない実験である。そのため、眼球位置で推定される電流の主な要因は固視微動(tremor)に起因するものである。tremorの振幅角は約0.0042度で,振動周波数は最大90Hzである。従来の実験と本実験タスクにおいて、単位時間当りの電荷の移動量(電流) の比は、4×0.7:0.0042×90=1:0.135となる.従来の実験の結果から換算すると、本実験の場合の眼球位置での推定電流の大きさは2.93−4.90nAmだと推測できる。この値の大きさも、本実験で示した値のそれと一致していた。
(実験結果1)
実験結果1における実験は、センサーレベルでのアーチファクト除去の実験である。本実験において、推定電流からアーチファクトの磁場を復元し、観測磁場から差し引くことによってアーチファクト除去を試みた。図21Aは、観測磁場とアーチファクトを除去した磁場を比較し、頭部マップ上に示した図である。アーチファクトを除去する前(黒色の線)は頭部外周部で波形が乱れている。しかし、アーチファクト除去後(灰色の線)は波形の乱れは小さい(図21C)。また、頭部中心付近は、外側部分に比べてアーチファクトの影響が小さいと予想されるが、そのような箇所では、アーチファクト除去の前後で波形は似通っている(図21B)。本実験において、脳活動の影響を主に反映しているセンサ位置の磁場は歪めずに、アーチファクトの影響の大きい位置でその影響を小さくすることができた。
観測磁場を試行全体に渡って加算平均された観測磁場はアーチファクトの影響が抑えられている。この値(Bハット[Bの上に^])は、真の脳活動の試行平均値と仮定できる。単一試行の観測磁場Bは、アーチファクトの影響もしくは試行毎の脳活動の揺らぎによってBハットからばらついている、と考えられる。各センサ位置での磁場のばらつきSDを、数式21を用いて計算した。
数式21において、Tはサンプル数、Nは試行数を表す。
図22は典型的な被験者のデータから計算した指標を頭部マップ上に表示した図である。観測磁場(図22A)では、アーチファクトの影響を受けやすい頭部外周部での値が大きいが、ノイズ除去した場合(PCA,ICA,extra−dipole法)は、その値が小さい(図22B,C,D)。また、頭部中心付近では、脳活動を反映した成分が観測磁場に含まれているため、ノイズ除去の前後に関わらず、ある程度大きな値が期待される。しかし、PCA(図22B)では、頭部中心付近の値も外周部と同様に小さい。PCA(図22B)では、アーチファクトだけではなく、脳活動を反映した成分も除去されている可能性がある。なお、図22Dのextra−dipole法は、脳活動情報出力装置12の方法である。
図23は頭部外周部(71個)のセンサの指標を観測値とPCA,ICA,extra−dipole法を用いた場合で比較した図である。4種類のデータから計算した指標に対して乱塊法による1要因分散分析を行った結果、有意差がみられた(F(3,1914)=496.64,p<0.0001)。多重比較(TukeyのHSD検定)を行ったところ、extra−dipole法を使ってアーチファクト除去をした場合が、他の解析法に比べて指標が有意に小さかった(図23)。つまり、脳活動情報出力装置12のノイズ除去方法を用いることで,アーチファクトの影響のみを最も効果的に除去できた、ことを示せた。
タスク遂行中の皮質電流の信号強度を、数式17を用いて計算し、その試行平均値を皮質マップ上に示した(図24)。被験者A,B,Cに共通して電流強度が大きかった領域として、中心前溝(precentral cortex)、内側上前頭回(medial superior frontal cortex)、頭頂間溝(intraparietal cortex)、外側後頭側頭野(Lateral occipitto−temporal cortex)が挙げられる。これらの領域はそれぞれ、前頭眼野(FEF)、補足眼野(SEF)、頭頂眼野(PEF)、MT/V5野といった眼球運動に関わる皮質領野である。また、これらの領域は目を動かさずに視野内の視覚対象に注意を向けるときにも同様に活動することが知られている。大きな電流活動が推定された部位は、心的な視標運動時に関係のある領域と一致していた。
(実験結果2)
実験結果2における実験は、推定電流レベルでのアーチファクト除去の実験である。本実験において、観測磁場は、複数のアーチファクトと皮質の神経活動を主に反映している。その磁場を説明できる個々の電流を推定することは、アーチファクトの影響を分離した脳活動を推定することにつながる。ここで、複数のアーチファクトと皮質電流を同時に推定することにより、電流レベルで推定精度が改善しているか否かを検討する。もし、脳活動とアーチファクトの影響を良く分離できていれば、推定された皮質電流は、実際の神経活動により近いはずである。ここでは、真の皮質電流の信号強度はfMRIの信号強度と相関があると仮定し、fMRIの活動強度(t値)と推定電流の信号強度(数式17で計算される値)の相関係数を、数式18を用いて計算した。また、他のノイズ除去方法の結果と比較するため、ここでは、4種類の方法で皮質電流を推定した。一つ目はダイポールを皮質のみに配置し観測磁場を用いて皮質電流を推定する方法、二つ目はPCAを用いてアーチファクトを除去した観測磁場から推定する方法、三つめはICAを用いてアーチファクトを除去した観測磁場から推定する方法、四つ目はダイポールを皮質表面だけではなくアーチファクト源にも配置し観測磁場を用いて電流推定を行う方法である。アーチファクト除去の精度評価には、アーチファクト源の電流分散分布ハイパーパラメータの決定には用いなかった残り2種類のデータセットを使った。推定された電流の善し悪しを調べるため、4種類の方法で求めた推定電流とfMRIの信号強度との相関係数を求め、乱塊法による1要因分散分析を行った。その結果、これらの相関係数の間に有意差がみられた(F(3,1368)=5928.48,p<0.0001)。多重比較(TukeyのHSD検定)を行った結果、拡張ダイポールを用いて推定した場合(脳活動情報出力装置12による方法の場合)にもっとも相関係数が大きかった(図25)。これは、皮質電流だけでなく、アーチファクトも同時に推定することで精度が向上したことを示している。
以上、本実施の形態における脳活動情報出力装置12によれば、脳活動とアーチファクトの振舞いを同時に予測し、ノイズ除去を行うことができる。そのため、精度の高い脳活動の情報を出力できる。
また、本実施の形態によれば、脳活動とその信号を歪める複数のアーチファクト源の電流波形を同時に推定し、実データからアーチファクトの影響を分離することで脳活動由来の信号を抽出することができた。また、本実施の形態によれば、適切な事前情報を設定することで、アーチファクト成分と皮質電流をより良く推定できた。
また、脳活動情報出力装置12のアーチファクト除去の手法を用いることで、センサーレベルでのアーチファクト除去に有効であること、さらに電流推定の精度が向上することが示せた。このとき、配置したダイポールの位置は左右の眼球、心臓、左右の頸動脈付近の5箇所であることは好適である。また、そのうち、心拍成分を表現するために配置したダイポールの位置はわずか3箇所である。それにもかかわらず、効果的にアーチファクトの影響を抑えることができた。
また、従来の方法(Y. Fujiwara, O. Yamashita, D. Kawawaki, K. Doya,M. Kawato, K. Toyama, and MA. Sato, "A hierarchical Bayesian method to resolve an inverse problem of MEG contaminated with eye movement artifacts," NeuroImage, 45, 2, pp. 393-409, 2009.)において、眼球位置での電流からEOGへの順モデルを事前に求め、その後、観測磁場を説明できる皮質電流を推定している。しかし、この方法を用いるためには、順モデルを求めるための情報(例えば、EOGなど)を別に用意する必要がある。生体由来のアーチファクトには様々なものが想定され、それらのMEGへの混入の仕方も複雑である。想定される様々なアーチファクト源の電位変化をすべて記録し、それらの順モデルを必要とする上記の手法は、複数のアーチファクトの影響を除去する手法としては現実的でない。それに対して、脳活動情報出力装置12の手法はアーチファクト源の特性が未知であっても適用することができ、また、その数が複数でも可能である。さらに、例えば、電源ノイズのような電位変化を計測できないアーチファクトに対しても有効である。
また、本実施の形態によれば、アーチファクト源や脳活動の信号強度は、それぞれ異なることを考慮し、信号強度に関する情報を事前知識として個別に設定し、電流の推定を行うことにより、より精度の高い脳活動の情報を出力できる。
また、本実施の形態によれば、推定した電流からアーチファクトの磁場を復元し、計測脳信号情報から差し引くことによってアーチファクト除去を行うことができ、精度の高い脳活動の情報を出力できる。
また、本実施の形態によれば、計測脳信号情報は、主として、脳磁図(MEG)であった。ただし、計測脳信号情報は、脳波(EEG)でも良いことは言うまでもない。
さらに、本実施の形態によれば、拡張ダイポールを用いてアーチファクト除去を行い、従来の方法と比べて合理的かつ効果的にセンサーレベルでノイズ除去ができる。また、アーチファクトも含めて脳電流を推定することによって、電流レベルでのクリーニングもできる。そのため、精度の高い脳活動の情報を出力できる。
なお、本実施の形態において、上述した数式の情報は、通常、当該数式を用いて演算する部、手段等が、予め保持しており、当該数式を読み出して、実行する。
また、出力脳活動情報は、数式15により算出される例を示した。しかし、出力脳活動情報は、以下に示すように算出されても良い。つまり、脳神経活動の電気的変化に伴って、頭部周辺では微弱な磁場が観測される。神経電流を位置と方向が固定されたIbrain個の電流ダイポールでモデル化し、その電流モーメントをIbrain次元ベクトルJbrainで表す。また、計測される磁場には、脳活動を反映した成分だけでなく、様々なアーチファクトに起因する成分も含まれている。それらのアーチファクト源をIartifact個のダイポールでモデル化し、その電流モーメントをIartifact次元ベクトルJartifactで表す。時刻τでの電流ダイポールと観測磁場の関係は数式22のように表すことができる。
数式22において、B(τ)はM個のMEGセンサ上での観測磁場を表す。Gbrainは皮質に関するリードフィールド行列であり、単位電流モーメントによって発生する磁場を表す行列である。脳を囲む頭蓋内部が均一な媒質の球であると仮定すると、リードフィールドの成分はSarvasの式によって解析的に計算できる。MEG信号を生成する主要な脳活動成分は錐体細胞に流れる電流であるという知見に基づき、皮質表面上に垂直な向きを持つ電流ダイポールを仮定した。アーチファクト源のダイポールに関するリードフィールド行列Gartifactは、アーチファクト源での磁場がx,y,z方向のベクトル量として表されると考え、Bio−Savartの式によって計算した。また、観測ノイズε(τ)は平均ゼロのガウス分布に従うと仮定した。
さらに、本実施の形態における処理は、ソフトウェアで実現しても良い。そして、このソフトウェアをソフトウェアダウンロード等により配布しても良い。また、このソフトウェアをCD−ROMなどの記録媒体に記録して流布しても良い。なお、このことは、本明細書における他の実施の形態においても該当する。なお、本実施の形態における脳活動情報出力装置を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、コンピュータを、アーチファクト源の信号強度に関する事前分布を取得し、かつ、当該事前分布を用いて当該アーチファクト源の信号強度に関する事後分布を推定し、前記事前分布と前記事後分布との差を算出し、当該差が予め決められた範囲以内の小さな差である信号強度に関する情報を有する1以上のアーチファクト事前情報を取得し、記憶媒体に蓄積するアーチファクト事前情報決定部と、記憶媒体に格納されている計測脳信号情報と脳事前情報とを用いて、脳活動の情報である脳活動情報を推定し、取得する脳活動情報取得部と、前記計測脳信号情報と、記憶媒体に格納されている1以上のアーチファクト事前情報とを用いて、アーチファクト源の情報である1以上のアーチファクト情報を推定し、取得するアーチファクト情報取得部と、前記計測脳信号情報と前記脳活動情報と前記1以上のアーチファクト情報とを用いて、前記計測脳信号情報からアーチファクトの影響を取り除く処理を行い、出力される脳活動の情報である出力脳活動情報を取得する出力脳活動情報取得部と、前記出力脳活動情報を出力する出力部として機能させるためのプログラム、である。
また、上記プログラムにおいて、前記アーチファクト事前情報決定部は、アーチファクト源の信号強度に関する事前分布を取得し、かつ、当該事前分布を用いて当該アーチファクト源の信号強度に関する事後分布を推定し、前記事前分布と前記事後分布との差を算出し、当該差が予め決められた範囲以内の小さな差である信号強度に関する情報を有する1以上のアーチファクト事前情報の候補を取得するアーチファクト信号強度情報推定手段と、fMRIの活動強度を格納し得るfMRI活動強度格納手段と、前記アーチファクト信号強度情報推定手段が取得した1以上の各アーチファクト事前情報が有する1以上の信号強度を取得し、当該1以上の各信号強度と前記fMRIの活動強度との相関に関する情報を算出し、当該相関に関する情報が最大の信号強度を有するアーチファクト事前情報を、前記アーチファクト信号強度情報推定手段が取得した1以上のアーチファクト事前情報の候補から選択し、当該選択したアーチファクト事前情報を前記アーチファクト事前情報格納部に蓄積する脳活動信号強度推定手段とを具備するように、コンピュータを、機能させるためのプログラム、であることは好適である。
また、上記プログラムにおいて、前記アーチファクト情報取得部は、複数のアーチファクト源ごとに、アーチファクト情報を推定し、複数のアーチファクト情報を取得し、前記出力脳活動情報取得部は、前記計測脳信号情報と前記脳活動情報と前記複数のアーチファクト情報とを用いて、前記計測脳信号情報から複数のアーチファクトの影響を取り除く処理を行い、出力される脳活動の情報である出力脳活動情報を取得するように、コンピュータを、機能させるためのプログラム、であることは好適である。
また、上記プログラムにおいて、前記脳事前情報の強度情報と、前記1以上の各アーチファクト事前情報の各強度情報は、異なることは好適である。
また、上記プログラムにおいて、前記1以上のアーチファクト事前情報は、被験者の生体に含まれるアーチファクト源に関する事前の情報である生体アーチファクト事前情報を含むことは好適である。
また、上記プログラムにおいて、前記1以上のアーチファクト事前情報は、被験者の生体に含まない外部のアーチファクト源に関する事前の情報である外部アーチファクト事前情報を含むことは好適である。
また、上記プログラムにおいて、前記脳事前情報は、fMRI強度情報を電流分散とし、当該電流分布に対する事前知識を表す事前分布であり、前記1以上の各アーチファクト事前情報は、fMRI強度情報と比較して大きな電流分散を用い、かつ当該電流分布に対する事前知識を表す事前分布であり、前記脳活動情報取得部は、複数のダイポール位置での電流についての情報である脳活動情報を、前記計測脳信号情報と前記脳事前情報とを用いた階層ベイズ推定を用いて推定し、取得し、前記アーチファクト情報取得部は、アーチファクト源のダイポール位置での電流についての情報である1以上のアーチファクト情報を、前記計測脳信号情報と、前記1以上のアーチファクト事前情報とを用いた階層ベイズ推定を用いて推定し、取得し、前記出力脳活動情報取得部は、前記計測脳信号情報から、前記脳活動情報のうちの脳の活動が無いと考えられる位置に配置されたダイポールに対応する脳活動情報と前記アーチファクト情報とを除く処理を行い、出力脳活動情報を取得するように、コンピュータを、機能させるためのプログラム、であることは好適である。
また、上記プログラムにおいて、前記計測脳信号情報は、頭部周辺で観測された磁場に関する情報であり、前記出力脳活動情報取得部は、前記脳活動情報のうちの脳の活動が無いと考えられる位置に配置されたダイポールに対応する脳活動情報である電流の情報から磁場に関する情報である背景脳磁場情報を取得し、かつ、前記アーチファクト情報である電流の情報から磁場に関する情報であるアーチファクト磁場情報を取得し、前記計測脳信号情報から、前記背景脳磁場情報と前記アーチファクト磁場情報とを除く処理を行い、出力脳活動情報を取得するように、コンピュータを、機能させるためのプログラム、であることは好適である。
さらに、上記プログラムにおいて、前記出力脳活動情報取得部が取得した出力脳活動情報と、前記計測脳信号情報とを用いて、当該出力脳活動情報におけるノイズ除去の精度に関する情報であるノイズ除去精度情報を取得するノイズ除去精度情報取得部をさらに具備し、前記出力部は、前記ノイズ除去精度情報が、予め格納されている閾値と比較して高い精度である場合のみ、前記出力脳活動情報を出力するように、コンピュータを、機能させるためのプログラム、であることは好適である。
また、図26は、本明細書で述べたプログラムを実行して、上述した実施の形態の脳活動情報出力装置等を実現するコンピュータの外観を示す。上述の実施の形態は、コンピュータハードウェア及びその上で実行されるコンピュータプログラムで実現され得る。図26は、このコンピュータシステム340の概観図であり、図27は、コンピュータシステム340のブロック図である。
図26において、コンピュータシステム340は、FDドライブ3411、CD−ROMドライブ3412を含むコンピュータ341と、キーボード342と、マウス343と、モニタ344とを含む。
図27において、コンピュータ341は、FDドライブ3411、CD−ROMドライブ3412に加えて、MPU3413と、CD−ROMドライブ3412及びFDドライブ3411に接続されたバス3414と、ブートアッププログラム等のプログラムを記憶するためのROM3415と、MPU3413に接続され、アプリケーションプログラムの命令を一時的に記憶するとともに一時記憶空間を提供するためのRAM3416と、アプリケーションプログラム、システムプログラム、及びデータを記憶するためのハードディスク3417とを含む。ここでは、図示しないが、コンピュータ341は、さらに、LANへの接続を提供するネットワークカードを含んでも良い。
コンピュータシステム340に、上述した実施の形態の脳活動情報出力装置等の機能を実行させるプログラムは、CD−ROM3501、またはFD3502に記憶されて、CD−ROMドライブ3412またはFDドライブ3411に挿入され、さらにハードディスク3417に転送されても良い。これに代えて、プログラムは、図示しないネットワークを介してコンピュータ341に送信され、ハードディスク3417に記憶されても良い。プログラムは実行の際にRAM3416にロードされる。プログラムは、CD−ROM3501、FD3502またはネットワークから直接、ロードされても良い。
プログラムは、コンピュータ341に、上述した実施の形態の脳活動情報出力装置等の機能を実行させるオペレーティングシステム(OS)、またはサードパーティープログラム等は、必ずしも含まなくても良い。プログラムは、制御された態様で適切な機能(モジュール)を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいれば良い。コンピュータシステム340がどのように動作するかは周知であり、詳細な説明は省略する。
また、上記プログラムを実行するコンピュータは、単数であってもよく、複数であってもよい。すなわち、集中処理を行ってもよく、あるいは分散処理を行ってもよい。
また、上記各実施の形態において、各処理(各機能)は、単一の装置(システム)によって集中処理されることによって実現されてもよく、あるいは、複数の装置によって分散処理されることによって実現されてもよい。
本発明は、以上の実施の形態に限定されることなく、種々の変更が可能であり、それらも本発明の範囲内に包含されるものであることは言うまでもない。