脳活動情報出力装置、脳活動情報出力方法、およびプログラム
【課題】従来、精度高く脳活動を推定することができなかった。
【解決手段】脳活動時光強度情報を複数の波長ごとに取得する脳活動時光強度情報取得部と、安静時光強度情報と脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、複数の波長ごと、およびプローブセットごとの光強度変化情報を、大脳皮質活動モデル情報と頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンとデオキシヘモグロビンの濃度変化を取得する推定部と、大脳皮質のオキシヘモグロビンとデオキシヘモグロビンの濃度変化を出力する出力部とを具備する脳活動情報出力装置により、精度高く脳活動を推定できる。
【解決手段】脳活動時光強度情報を複数の波長ごとに取得する脳活動時光強度情報取得部と、安静時光強度情報と脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、複数の波長ごと、およびプローブセットごとの光強度変化情報を、大脳皮質活動モデル情報と頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンとデオキシヘモグロビンの濃度変化を取得する推定部と、大脳皮質のオキシヘモグロビンとデオキシヘモグロビンの濃度変化を出力する出力部とを具備する脳活動情報出力装置により、精度高く脳活動を推定できる。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、脳の活動に関する情報を取得し、出力する脳活動情報出力装置等に関するものである。
【背景技術】
【０００２】
従来、ＮＩＲＳ（Ｎｅａｒ−Ｉｎｆｒａｒｅｄ Ｓｐｅｃｔｒｏｓｃｏｐｙ）データから推定される脳活動の表示技術があった。かかる技術は、計測した光強度変化をヘモグロビン濃度に変換し、２次元のトポグラフィとして表示する方法であった。
【０００３】
また、近年は、拡散光トモグラフィ（ＤＯＴ）を用いて推定したものを３次元のトモグラフィとして表示するやり方が広まってきている（非特許文献１、非特許文献２参照）。
【０００４】
また、これまでのＤＯＴでは、波長毎に単独で吸収係数の変化を推定し、それを変換することで、ヘモグロビンの濃度を推定していた（非特許文献２参照）。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】D. A. Boas and A. M. Dale, "Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function", Applied Optics 44 (2005) pp. 1957-1968.
【非特許文献２】B. W. Zeff, B. R. White, H. Dehghani, B. L. Schlaggar,and J. P. Culver, "Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography",PNAS 104 (2007) pp. 12169-12174.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、従来の技術においては、大脳皮質における血行動態反応は局所的であり、頭皮における血行動態反応は広域的であるという知識をうまく活用することができず、頭皮血流に起因するアーチファクトを適切に除去できないために、精度高く脳活動を推定することができなかった。また、従来の技術においては、波長毎の観測値からそれぞれ別々に吸収係数の変化を推定し、それを変換することで、ヘモグロビンの濃度を推定していたため、精度高く脳活動を推定することができなかった。なお、ここでの脳活動の推定とは、神経活動に伴う脳血行動態反応を推定することである。
【課題を解決するための手段】
【０００７】
本第一の発明の脳活動情報出力装置は、局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報を格納している大脳皮質活動モデル情報格納部と、広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報を格納している頭皮血流変化モデル情報格納部と、送光プローブと受光プローブの組であるプローブセットを１組以上有する近赤外分光計測装置を用いて、各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納している安静時光強度情報格納部と、近赤外分光計測装置を用いて、各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、安静時光強度情報格納部の安静時光強度情報と、脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、複数の波長ごと、およびプローブセットごとの光強度変化情報を、大脳皮質活動モデル情報および頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部とを具備する脳活動情報出力装置である。
【０００８】
かかる構成により、精度高く脳活動を推定できる。
【０００９】
また、本第二の発明の脳活動情報出力装置は、第一の発明に対して、推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式を格納し得る演算式格納手段と、演算式に、感度情報と光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する推定手段とを具備する脳活動情報出力装置である。
【００１０】
かかる構成により、極めて精度高く脳活動を推定できる。
【００１１】
また、本第三の発明の脳活動情報出力装置は、第一の発明に対して、推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の情報の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式であり、コスト関数である演算式を格納し得る演算式格納手段と、演算式に、感度情報と光強度変化情報を代入し、演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定手段とを具備する脳活動情報出力装置である。
【００１２】
かかる構成により、精度高く脳活動を推定できる。
【発明の効果】
【００１３】
本発明による脳活動情報出力装置によれば、精度高く脳活動を推定できる。
【図面の簡単な説明】
【００１４】
【図１】実施の形態１における脳活動情報出力システムの概念図
【図２】同脳活動情報出力システムのブロック図
【図３】同脳活動情報出力装置の動作について説明するフローチャート
【図４】同変分ベイズ推定法について説明するフローチャート
【図５】同プローブの配置を説明する図
【図６】同実験におけるプローブの配置を説明する図
【図７】同実験におけるチャネルを説明する図
【図８】同波長７８０ｎｍにおける光量の計測値等を示す図
【図９】同共分散行列を示す図
【図１０】同実験における各波長の光強度変化情報を示す図
【図１１】同観測値をビジュアルに示した図
【図１２】同推定手段が算出した頭皮（深さ１）のオキシヘモグロビンの濃度変化の推定値を示す図
【図１３】同推定手段が算出した大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の推定値を示す図
【図１４】同推定手段が算出した頭皮（深さ１）のデオキシヘモグロビンの濃度変化の推定値を示す図
【図１５】同推定手段が算出した大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の推定値を示す図
【図１６】同出力部の出力例を示す図
【図１７】同出力部の出力例を示す図
【図１８】同シミュレーションで仮定した真の頭皮（深さ１）のオキシヘモグロビンの濃度変化の値を示す図
【図１９】同シミュレーションで仮定した真の大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の値を示す図
【図２０】同シミュレーションで仮定した真の頭皮（深さ１）のデオキシヘモグロビンの濃度変化の値を示す図
【図２１】同シミュレーションで仮定した真の大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の値を示す図
【図２２】同脳活動情報出力装置の推定方法の概念について説明図
【図２３】同コンピュータシステムの概観図
【図２４】同コンピュータシステムのブロック図
【図２５】同波長７８０ｎｍにおける光量の計測値等を示す図
【図２６】同共分散行列を示す図
【図２７】同実験における各波長の光強度変化情報を示す図
【図２８】同観測値をビジュアルに示した図
【図２９】同推定手段が算出した頭皮（深さ１）のオキシヘモグロビンの濃度変化の推定値を示す図
【図３０】同推定手段が算出した大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の推定値を示す図
【図３１】同推定手段が算出した頭皮（深さ１）のデオキシヘモグロビンの濃度変化の推定値を示す図
【図３２】同推定手段が算出した大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の推定値を示す図
【図３３】同出力部の出力例を示す図
【図３４】同出力部の出力例を示す図
【発明を実施するための形態】
【００１５】
以下、脳活動情報出力装置等の実施形態について図面を参照して説明する。なお、実施の形態において同じ符号を付した構成要素は同様の動作を行うので、再度の説明を省略する場合がある。
【００１６】
（実施の形態１）
本実施の形態において、脳活動を推定する脳活動情報出力システムについて説明する。また、脳活動情報出力システムは、近赤外分光計測装置（ＮＩＲＳ）を具備する。なお、脳活動の推定とは、神経活動に伴う脳血行動態反応を推定することである。
【００１７】
また、本実施の形態における脳活動情報出力システムにおいて、拡散光トモグラフィ（ＤＯＴ）を用いる。ＤＯＴとは、頭部内での光子移動を拡散過程としてモデル化することにより、脳内活動と観測光強度変化の関係を物理的に正しく求め、それを利用して脳活動の推定の精度を高めようというものである。
【００１８】
また、本実施の形態における脳活動情報出力システムにおいて、特に、局所的な活動である大脳皮質の活動、および広域的な活動である頭皮の血流変化をモデル化しており、かかるモデルを用いて、脳活動を推定する方法を採用する。また、特に、本実施の形態において、頭皮血流をＤＯＴのモデルに組み込んで、大脳皮質における神経活動と頭皮血流に起因するアーチファクトとを同時推定することにより、頭皮血流に起因するアーチファクトを除去し、脳活動を推定する階層変分ベイズ法について説明する。なお、後述する数式に出現する上付きのcは大脳皮質(cortex),sは頭皮(scalp)を表す。
【００１９】
さらに、具体的には、本実施の形態における脳活動情報出力システムにおいて、近赤外分光計測装置における観測結果から、神経活動に伴う大脳皮質内のオキシヘモグロビン（Ｏｘｙ Ｈｂ）、デオキシヘモグロビン（Ｄｅｏｘｙ Ｈｂ）の濃度変化の場所と大きさを推定する。
【００２０】
図１は、本実施の形態における脳活動情報出力システムの概念図である。脳活動情報出力システムは、近赤外分光計測装置１、および脳活動情報出力装置２を具備する。
【００２１】
また、図２は、本実施の形態における脳活動情報出力システムのブロック図である。
【００２２】
近赤外分光計測装置１は、送光プローブ１１、受光プローブ１２、受付部１３、および制御部１４を具備する。近赤外分光計測装置１は、送光プローブ１１と受光プローブ１２の組であるプローブセット（チャネルとも言う）を１組以上有する。
【００２３】
脳活動情報出力装置２は、大脳皮質活動モデル情報格納部２１、頭皮血流変化モデル情報格納部２２、安静時光強度情報格納部２３、安静時光強度取得部２４、安静時光強度蓄積部２５、脳活動時光強度情報取得部２６、光強度変化情報取得部２７、推定部２８、および出力部２９を具備する。
【００２４】
また、推定部２８は、感度情報格納手段２８１、演算式格納手段２８２、および推定手段２８３を具備する。
【００２５】
近赤外分光計測装置１を構成する送光プローブ１１は、近赤外光を生体（ここでは、特に頭部）に照射する。
【００２６】
受光プローブ１２は、生体（ここでは、特に頭部）内部で吸収あるいは散乱され、生体内から出てくる光を受光する。なお、プローブの配置密度が３ｃｍ未満であることは好適である。また、プローブの配置密度が２ｃｍ以下であることはさらに好適である。また、プローブの配置密度が１．７ｃｍ以下であることは好適である。さらに、プローブの配置密度が１．５ｃｍ以下であることは極めて好適である。プローブの配置密度とは、送光プローブ１１と受光プローブ１２間の間隔である。なお、プローブの配置密度は密であればある程良いのは、頭皮血流をとらえるための短い距離のチャネルがたくさんできること、および脳活動をとらえるための３〜４ｃｍ程度の中距離のチャネルがたくさんできること等からである。
【００２７】
受付部１３は、外部からの指示を受け付ける。かかる指示は、近赤外光を生体に照射する指示である照射指示などである。
【００２８】
制御部１４は、受付部１３が照射指示を受け付けた場合に、送光プローブ１１に近赤外光の照射を指示する。なお、送光プローブ１１は、通常、制御部１４から照射指示を受け付けた場合に、近赤外光を生体に照射する。
【００２９】
脳活動情報出力装置２は、近赤外分光計測装置１が計測した光強度に関する情報を用いて、脳活動に関する情報である脳活動情報を取得し、出力する装置である。ここで、出力とは、ディスプレイへの表示、プロジェクターを用いた投影、プリンタへの印字、音出力、外部の装置への送信、記録媒体への蓄積、他の処理装置や他のプログラムなどへの処理結果の引渡しなどを含む概念である。
【００３０】
脳活動情報出力装置２を構成する大脳皮質活動モデル情報格納部２１は、大脳皮質活動モデル情報を格納している。大脳皮質活動モデル情報は、局所的な活動である大脳皮質の活動のモデルを示す情報である。大脳皮質活動モデル情報は、例えば、後述する局所的フィルタＷＣ、後述する正則化の処理を示すプログラムやフラグ、大脳皮質の各ボクセルにおける正則化（大きな値にペナルティ）に関するガンマ分布のシェイプパラメータ（後述するγα０ｉ）、大脳皮質の各ボクセルにおける正則化に関するガンマ分布の平均値（後述するα０ｉ−）等である。なお、大脳皮質活動モデル情報に対応する正則化は、ボクセル毎に異なった正則化である。また、ここでいう正則化とは、大きな正則化対象の値についてペナルティをかけることである。
【００３１】
頭皮血流変化モデル情報格納部２２は、頭皮血流変化モデル情報を格納している。頭皮血流変化モデル情報は、広域的な活動である頭皮の血流変化のモデルを示す情報である。頭皮血流変化モデル情報は、例えば、後述する広域的フィルタＷＳ、後述する正則化の処理を示すプログラムやフラグ、頭皮の正則化に関するガンマ分布のシェイプパラメータ（後述するγβ０）、頭皮の正則化に関するガンマ分布の平均値（後述するβ０−）等である。なお、頭皮血流変化モデル情報に対応する正則化は、全部のボクセルで共通な正則化である。
【００３２】
安静時光強度情報格納部２３は、複数の各波長に対応付けて安静時光強度情報を格納している。安静時光強度情報は、送光プローブ１１と受光プローブ１２の組であるプローブセットを１組以上有する近赤外分光計測装置１を用いて、各プローブセットを構成する各送光プローブ１１から、一の光強度の光を発した場合に、安静状態において測定された各受光プローブ１２の光強度に関する情報である
【００３３】
安静時光強度取得部２４は、安静状態の被験者が装着している近赤外分光計測装置１の受光プローブ１２から取得した光強度に関する情報である安静時光強度情報を、近赤外分光計測装置１から取得する。安静時光強度取得部２４は、複数の各波長ごとに、安静時光強度情報を取得する。
【００３４】
安静時光強度蓄積部２５は、安静時光強度取得部２４が取得した安静時光強度情報を、複数の各波長に対応付けて安静時光強度情報格納部２３に蓄積する。
【００３５】
脳活動時光強度情報取得部２６は、複数の各波長ごとに脳活動時光強度情報を取得する。脳活動時光強度情報は、脳活動が行われている状態において測定された各受光プローブ１２の光強度に関する情報である。脳活動時光強度情報は、被験者の脳の活動が行われている状態において、近赤外分光計測装置１を用いて、各プローブセットを構成する各送光プローブ１１から、一の光強度の光を発した場合に、各受光プローブ１２が受光した光強度に関する情報である。
【００３６】
光強度変化情報取得部２７は、安静時光強度情報格納部２３の安静時光強度情報と、脳活動時光強度情報取得部２６が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度変化情報を取得する。光強度変化情報は、光強度の変化に関する情報である。光強度変化情報は、例えば、光強度の変化率、光強度の変化率の対数をとった値、光強度の変化量等である。また、光強度変化情報取得部２７は、光強度のゆらぎ（自発脳活動やセンサーの機器ノイズにより光強度に加わるノイズ）の情報Σyも取得しても良い。
【００３７】
推定部２８は、複数の波長ごと、およびプローブセットごとの光強度変化情報を、大脳皮質活動モデル情報および頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する。なお、推定部２８は、推定する値（後述するＸ）の数より少ない観測値（後述するＹ）を用いて、値（後述するＸ）を推定し、取得する。なお、かかる推定部２８が解く問題を逆問題という。推定部２８は、後述する変分ベイズ法を用いて逆問題を解くことは好適であるが、逆問題の解法は問わない。また、推定部２８は、頭部の領域ごとに、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する。頭部の領域ごととは、例えば、頭部のボクセルごと、または頭部を三角錐に区切った領域ごと等である。
【００３８】
感度情報格納手段２８１は、感度情報を格納し得る。感度情報とは、大脳皮質または頭皮の活動の変化（例えば、後述するＸが対応）と検出された光強度の変化の情報（例えば、後述するＹが対応）の関係を示す情報である。感度情報は、例えば、後述する感度行列である。感度情報は、推定部２８が算出し、感度情報格納手段２８１に、少なくとも一時的に格納しても良い。
【００３９】
演算式格納手段２８２は、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式（例えば、演算式１）であり、ベイズ推定の演算式を格納し得る。また、演算式格納手段２８２は、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式（例えば、演算式２）であり、コスト関数である演算式を格納し得る。なお、演算式１は、例えば、後述する数式２８である。また、演算式２は、例えば、後述する数式３１である。
【００４０】
推定手段２８３は、演算式（例えば、演算式１）に、感度情報と光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する。また、推定手段２８３は、演算式（例えば、演算式２）に、感度情報と光強度変化情報を代入し、演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得しても良い。なお、「コストが最小である」とは、コストに「−１」をかけた値が最大であることも含む。つまり、「コストが最小である」とは、実質的に考慮する。
【００４１】
出力部２９は、推定部２８が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。なお、出力部２９は、頭皮の血管のオキシヘモグロビンの濃度変化、および頭皮の血管のデオキシヘモグロビンの濃度変化も出力しても良い。また、出力部２９は、大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を、数値を用いて出力しても良いし、図的に出力しても良い。その出力態様は問わない。
【００４２】
大脳皮質活動モデル情報格納部２１、頭皮血流変化モデル情報格納部２２、安静時光強度情報格納部２３、感度情報格納手段２８１、および演算式格納手段２８２は、不揮発性の記録媒体が好適であるが、揮発性の記録媒体でも実現可能である。
【００４３】
大脳皮質活動モデル情報格納部２１等に大脳皮質活動モデル情報等が記憶される過程は問わない。例えば、記録媒体を介して大脳皮質活動モデル情報等が大脳皮質活動モデル情報格納部２１等で記憶されるようになってもよく、通信回線等を介して送信された大脳皮質活動モデル情報等が大脳皮質活動モデル情報格納部２１等で記憶されるようになってもよく、あるいは、入力デバイスを介して入力された大脳皮質活動モデル情報等が大脳皮質活動モデル情報格納部２１等で記憶されるようになってもよい。
【００４４】
安静時光強度取得部２４は、通常、受信手段等から実現され得る。
【００４５】
安静時光強度蓄積部２５、脳活動時光強度情報取得部２６、光強度変化情報取得部２７、推定部２８、および推定手段２８３は、通常、ＭＰＵやメモリ等から実現され得る。安静時光強度蓄積部２５等の処理手順は、通常、ソフトウェアで実現され、当該ソフトウェアはＲＯＭ等の記録媒体に記録されている。但し、ハードウェア（専用回路）で実現しても良い。
【００４６】
出力部２９は、ディスプレイやスピーカー等の出力デバイスを含むと考えても含まないと考えても良い。出力部２９は、出力デバイスのドライバーソフトまたは、出力デバイスのドライバーソフトと出力デバイス等で実現され得る。
【００４７】
次に、脳活動情報出力システムを構成する脳活動情報出力装置２の動作について、図３のフローチャートを用いて説明する。なお、図３のフローチャートの開始前において、安静時光強度取得部２４は、複数の波長ごとに、安静時光強度情報を取得し、安静時光強度蓄積部２５は、複数の波長ごとに、当該安静時光強度情報を安静時光強度情報格納部２３に蓄積した、とする。そして、安静時光強度情報格納部２３には、複数の各波長ごとに、安静時光強度情報が格納されている、とする。
【００４８】
（ステップＳ３０１）脳活動時光強度情報取得部２６は、複数の各波長ごとに脳活動時光強度情報を取得したか否かを判断する。複数の各波長の脳活動時光強度情報を取得すればステップＳ３０２に行き、取得しなければステップＳ３０１に戻る。
【００４９】
（ステップＳ３０２）光強度変化情報取得部２７は、安静時光強度情報格納部２３に格納されている安静時光強度情報を、複数の各波長に対応付けて読み出す。
【００５０】
（ステップＳ３０３）光強度変化情報取得部２７は、ステップＳ３０１で取得された各波長ごとの複数の脳活動時光強度情報と、ステップＳ３０２で読み出された波長ごとの複数の安静時光強度情報とから、波長ごとに、光強度変化情報を算出する。光強度変化情報取得部２７は、例えば、波長ごとに、安静時光強度情報から脳活動時光強度情報を減算することにより、光強度変化情報を取得する。
【００５１】
（ステップＳ３０４）推定手段２８３は、演算式格納手段２８２から演算式を読み出す。
【００５２】
（ステップＳ３０５）推定手段２８３は、感度情報格納手段２８１から感度情報を読み出す。なお、例えば、推定手段２８３は予め感度情報を算出し、感度情報格納手段２８１に格納している、とする。
【００５３】
（ステップＳ３０６）推定手段２８３は、大脳皮質活動モデル情報格納部２１から大脳皮質活動モデル情報を読み出す。
【００５４】
（ステップＳ３０７）推定手段２８３は、頭皮血流変化モデル情報格納部２２から頭皮血流変化モデル情報を読み出す。
【００５５】
（ステップＳ３０８）推定手段２８３は、演算式に、光強度変化情報、感度情報を代入し、かつ、大脳皮質活動モデル情報と頭皮血流変化モデル情報とを用いて、変分ベイズ推定法により、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する。かかる変分ベイズ推定法について、図４のフローチャートを用いて説明する。
【００５６】
（ステップＳ３０９）出力部２９は、推定部２８が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。ステップＳ３０１に戻る。
【００５７】
なお、図３のフローチャートにおいて、電源オフや処理終了の割り込みにより処理は終了する。
【００５８】
次に、ステップＳ３０８の変分ベイズ推定法について、図４のフローチャートを用いて説明する。
【００５９】
（ステップＳ４０１）推定手段２８３は、予め格納されている初期値を読み出す。読み出す初期値は、は、γα０ｉ、α０ｉ−（−はαの真上に配置）、γβ０、β０−（−はβの真上に配置）である。また、初期値は、ユーザが与えても良いし、他の装置から受信等しても良い。ここで、γα０ｉは大脳皮質の各ボクセルにおける正則化（大きな値にペナルティ）に関するガンマ分布のシェイプパラメータの初期値、α０ｉ−は大脳皮質の各ボクセルにおける正則化に関するガンマ分布の平均値の初期値、γβ０は頭皮の正則化に関するガンマ分布のシェイプパラメータの初期値、β０−は頭皮の正則化に関するガンマ分布の平均値の初期値である。
【００６０】
（ステップＳ４０２）推定手段２８３は、後述する第一ステップの演算を行う。第一ステップの演算は、変分ベイズ法における、「Ｘ，σ最大化ステップ」である。
【００６１】
（ステップＳ４０３）推定手段２８３は、後述する第二ステップの演算を行う。第二ステップの演算は、変分ベイズ法における、「Ａ，β最大化ステップ」である。
【００６２】
（ステップＳ４０４）推定手段２８３は、大脳皮質の各ボクセルにおける正則化に関するガンマ分布のシェイプパラメータ（γαｉ）、大脳皮質の各ボクセルにおける正則化に関するガンマ分布の平均値（αｉ−）、頭皮の正則化に関するガンマ分布のシェイプパラメータ（γβ）、および頭皮の正則化に関するガンマ分布の平均値（β−）の値が収束したか否かを判断する。収束すればステップＳ４０５に行き、収束しなければステップＳ４０２に戻る。なお、収束したか否かは、予め決められた条件を満たすほど、上記の４つの変数の各値の変化が小さくなったか否かにより判断される。
【００６３】
（ステップＳ４０５）推定手段２８３は、第一ステップの演算を行う。
【００６４】
（ステップＳ４０６）推定手段２８３は、推定結果を算出する演算を行う。処理を終了する。
【００６５】
以下、本実施の形態における脳活動情報出力システムの具体的な動作について説明する。脳活動情報出力システムの概念図は、図１である。
【００６６】
まず、本脳活動情報出力システムにおいて、頭部が半無限一層モデルであると仮定し、計測された光強度から脳活動を推定する場合について説明する。
【００６７】
また、プローブは、図５に示すように、１５ｍｍ間隔で頭部モデルの表面に配置されている、とする。図５において、黒丸が送光プローブ１１、白丸が受光プローブ１２を表す。また、図５において、直方体で囲ってある部分はこれから推定値を求めていく領域（頭部）で、一辺４ｍｍからなる立方体のボクセルに分けて考える。図５において、頭部は、縦１１ボクセル、横１１ボクセル、高さ５ボクセルに分割されている。ここで、深さ１を頭皮、深さ２・３を頭蓋骨・脳脊髄液、深さ４以降を大脳皮質とする。
【００６８】
また、今回の頭部モデルの光学特性は、ヒト灰白質と同じであり、換算散乱係数「μｓ'＝１．１１ｍｍ−１」、吸収係数「μａ＝０．０１８６ｍｍ−１」とする。さらに、表面境界における屈折率は「ｎ＝１．４」とする。
【００６９】
上記の条件で、被験者は安静な状態を保っており、かかる状態において、安静時光強度取得部２４は、複数の波長ごとに、安静時光強度情報を取得し、安静時光強度蓄積部２５は、複数の波長ごとに、当該安静時光強度情報を安静時光強度情報格納部２３に蓄積した、とする。そして、例えば、安静時光強度情報格納部２３には、波長＝７８０ｎｍ、８３０ｎｍ（以下、適宜、単位を省略する。）の２つの各波長ごとに、安静時光強度情報が格納された、とする。なお、安静時とは、脳活動や頭皮血流変化があまりない安静状態をいう。
【００７０】
そして、一般に、光強度は機器ノイズや自発脳活動によってゆらぐので、各波長ごと、各チャネルごとに多数計測し、安静時光強度取得部２４は、その平均値と分散共分散行列を算出することは好適である。ここで、光強度の平均値は、数式１のように示される。また、分散共分散行列はΣｙである。
【数１】
【００７１】
なお、数式１において、Φ０は安静時の光強度であり、ｒｓ→（ｒの真上に→が存在）は送光プローブ１１の位置、ｒｄ→（ｒの真上に→が存在）は受光プローブ１２の位置を示す。数式１のバーは、平均値を意味する。
【００７２】
次に、脳活動時光強度情報取得部２６は、各チャネルごと、波長＝７８０、８３０の２つの各波長ごとに、脳活動時光強度情報を取得した、とする。具体的には、脳活動時光強度情報取得部２６は、数式２に示す脳活動時光強度情報を、チャネルごと、波長ごとに取得する。また、脳活動時光強度は複数回計測して取得してもよい。その計測の番号をサンプル番号とする。
【数２】
【００７３】
次に、光強度変化情報取得部２７は、安静時光強度情報格納部２３に格納されている２つの安静時光強度情報を、各波長（波長＝７８０、８３０）に対応付けて読み出す。
【００７４】
次に、光強度変化情報取得部２７は、取得された波長ごとの２つの脳活動時光強度情報と、読み出された波長ごとの２つの安静時光強度情報とから、波長ごとに、光強度変化情報を、以下のように算出する。
【００７５】
つまり、光強度変化情報取得部２７は、ここでは、数式３の演算式を格納しており、それぞれのサンプル番号ごとに、かかる数式３に、安静時光強度情報（ｌｏｇΦ０）と脳活動時光強度情報（ｌｏｇΦ）を代入し、光強度変化情報（Ｙ）を取得する。
【数３】
【００７６】
ここでtはサンプル番号を表す。なお、ＹはＹ（ｔ）を各列に持つ行列であり、その次元は「（チャネル数×波長数）×（サンプル数）」である。今回は、サンプル数は１とする。従って、ここでは、例えば、Ｙは４０次元である、とする。
【００７７】
また、ここで、光強度変化情報取得部２７は、数式４を格納しており、数式４の演算式に、安静時光強度取得部２４が取得した安静時光強度の平均値と、脳活動時光強度情報（数式２）を代入し、光強度変化情報Ｙを算出することは好適である。
【数４】
【００７８】
そして、頭部モデルと、図５に示すようにプローブ配置（送光プローブ１１と受光プローブ１２の配置）が決まっているので、推定手段２８３は、光拡散過程を解くことにより、以下の数式５を用いて、各波長ごとの吸収係数変化に関する感度行列Ｇを計算して、取得する。ここで、感度行列Ｇは、感度情報の一例である。なお、数式５において、定義の仕方によっては定数倍がかかる場合がある。
【数５】
【００７９】
なお、数式５において、jはチャネル（送光プローブ１１と受光プローブ１２のペア）のインデックス、iは頭部ボクセルのインデックスである。また、数式５において、Φ０は安静時の光強度であり、Φ０（第一変数，第二変数）は、第一変数の位置から出発し、第二変数の位置に到達する光の強度を示す。なお、Φ０は頭部モデル内での光子の拡散過程を解くことで得られる。また、ｒｓ， j→（ｒの真上に→が存在）はチャネルjの送光プローブ１１の位置を、ｒｄ，j→（ｒの真上に→が存在）はチャネルjの受光プローブ１２の位置を示す。また、ｒi→（ｒの真上に→が存在）は頭部ボクセルiの位置を示す。
【００８０】
また、ここでは、推定手段２８３は、感度行列を、頭皮（図５の深さ１の領域）と大脳皮質（図５の深さ４または深さ５の領域）のみで求める。図５において、頭蓋骨・脳脊髄液（深さ２または深さ３の領域）にはヘモグロビンとオキシヘモグロビンが存在せず、常にＸ＝０だから、感度行列Ｇ＝０である。また、一般的に、頭部モデルの光学特性は波長ごとに違うので、推定手段２８３は、感度行列を波長毎に算出する。なお、Ｘは、ここでは、オキシヘモグロビンの濃度変化およびデオキシヘモグロビンの濃度変化である。
【００８１】
つまり、推定手段２８３は、Ｇｃ７８０、Ｇｓ７８０、Ｇｃ８３０、Ｇｓ８３０を算出し、感度情報格納手段２８１に蓄積しておく。ここで、上付きのｃは大脳皮質（深さ４または深さ５の領域）、ｓは頭皮（深さ１の領域）を意味する。下付きの７８０は波長７８０ｎｍ、８３０は波長８３０ｎｍを意味する。つまり、Ｇｃ７８０は、波長７８０ｎｍにおける大脳皮質の感度行列である。
【００８２】
次に、吸収係数変化に関する感度行列から、ヘモグロビン濃度変化に関する感度行列を算出する。その変換式は、取り扱う波長が（波長＝７８０、８３０）の２つの場合、以下の数式６に示すように詳細に定義される。
【数６】
【００８３】
数式６において、δμａは脳活動時の吸収係数と安静時の吸収係数の差である。また、δμａは波長毎に値が異なる。つまり、δμａ，７８０は、波長７８０ｎｍにおける、脳活動時の吸収係数と安静時の吸収係数の差である。また、数式６において、εは、モル吸光係数であり、波長や物質（オキシヘモグロビン、デオキシヘモグロビン）により異なる。なお、数式８は、波長ごと、物質ごとのモル吸光係数の値を示す。
【００８４】
また、それぞれの波長における吸収係数の変化は、オキシヘモグロビンの濃度変化Ｘｏｘｙと、デオキシヘモグロビンの濃度変化Ｘｄｅｏｘｙとを用いて、数式７のように表される。そして、数式６の（１）は、数式７を用いて、数式６の（２）に変形される。そして、数式６の（２）は、計算により、数式６の（３）に変形される。そして、ヘモグロビン濃度変化に関する感度行列Ｇは（Ｇｃｏｘｙ，Ｇｃｄｅｏｘｙ，Ｇｓｏｘｙ，Ｇｓｄｅｏｘｙ）である。
【数７】
【数８】
【００８５】
次に、推定手段２８３は、感度情報格納手段２８１から感度行列Ｇ（Ｇｃｏｘｙ，Ｇｃｄｅｏｘｙ，Ｇｓｏｘｙ，Ｇｓｄｅｏｘｙ）を読み出す。なお、Ｇｃ７８０などは吸収係数変化に関する感度行列、Ｇｃｏｘｙなどは、ヘモグロビン濃度変化に関する感度行列である。また、吸収係数変化に関する感度行列とヘモグロビン濃度変化に関する感度行列とは、上述したように、数式６の定義式により簡単に変換できる。なお、数式５を用いて、吸収係数変化に関する感度行列が算出され、次に、当該感度行列にモル吸光係数などをかけて、ヘモグロビン濃度変化に関する感度行列が算出される。そして、ヘモグロビン濃度変化に関する感度行列が上記の推定に用いられる。
【００８６】
次に、推定手段２８３は、大脳皮質活動モデル情報格納部２１から大脳皮質活動モデル情報（例えば、局所的フィルタＷＣ等）を読み出す。また、推定手段２８３は、頭皮血流変化モデル情報格納部２２から頭皮血流変化モデル情報（例えば、広域的フィルタＷＳ等）を読み出す。
【００８７】
次に、推定手段２８３は、以下のように２つのステップ（第一ステップ、第二ステップ）により、変分ベイズ推定を用いて、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を推定する。なお、推定手段２８３は、後述する演算式等を、通常、予め格納している、とする。
【００８８】
まず、第一ステップとは、以下の数式９から数式１５を用いて、Σ、σ―（−はσの上に配置）、γσ、Δ、＜（ΔＸ）２＞を算出する処理である。ここで、Σは、数式９により定義される。また、Σｙは安静時の観測値Ｙの分散共分散行列である。また、σは脳活動時の分散共分散行列の安静時に比べての割合の逆数である。つまり、脳活動時の観測値Yの分散共分散行列はσ−１Σｙであるとする。なお、初期値として、「Ｐ０（σ）∝１／σ」のような分布に従うと仮定するが、アルゴリズム途中ではガンマ分布にしたがう。その平均値がσ―である。また、シェイプパラメータをγσとする。なお、数式１０のγσは、自由エネルギーが最大となるような値である。
【数９】
【００８９】
なお、行列の上付けのＴは、行列の転置をあらわす。また、数式９において、Ｗはスムージングフィルタであり、Ｗｃは局所的フィルタ、Ｗｓは広域的フィルタを示す。具体的には、例えば、Ｗｃは半値幅８ｍｍのガウシアンフィルタを用い、Ｗｓはラプラス行列の逆行列を用いる。
【数１０】
【００９０】
なお、数式１０において、Ｎはチャネル数、Ｔはサンプル数である。また、Ｌは用いる波長の数である。また、Ｔｒは行列のトレースをあらわす。
【数１１】
【数１２】
【００９１】
なお、数式１２において、Ｉは単位行列を示す。
【００９２】
また、数式１１、１２において、Ａはαｉを対角要素に持つ行列（αｉは大脳皮質の各ボクセル毎に在る）である。
【００９３】
つまり、ＡおよびＡ−（−はＡの真上に存在する）は、数式１３のようになる。
【数１３】
【００９４】
数式１３において、αｉはそれぞれのボクセルiにおける、正則化（大きな値にペナルティ）に関するパラメータである。そして、αｉはガンマ分布に従うとする。そして、αｉ−（−はαの真上に存在する）はガンマ分布の平均値である。また、γαｉはガンマ分布のシェイプパラメータである。また、大脳皮質ボクセル数がＭｃ、頭皮ボクセル数がＭｓである。なお、ガンマ分布は、数式１４である。
【数１４】
なお、数式１４においてΓはガンマ関数をあらわす。
【００９５】
また、数式１１、１２において、βは頭皮における正則化に関するパラメータである。βはガンマ分布に従うとする。そしてβ−（−はβの真上に存在する）はガンマ分布の平均値、γβはガンマ分布のシェイプパラメータである。また、ガンマ分布は、数式１５である。
【数１５】
【００９６】
次に、推定手段２８３は、第二ステップを実行する。第二ステップとは、γαｉ、αｉ−（−はαの真上に配置）、γβ、β−（−はβの真上に配置）を、以下の数式１６から数式１９を用いて算出する処理である。なお、γαｉは大脳皮質のボクセルにおけるガンマ分布のシェイプパラメータ、αｉ−は大脳皮質のボクセルにおけるガンマ分布の平均値、γβは頭皮のボクセルにおけるガンマ分布のシェイプパラメータ、β−は頭皮のボクセルにおけるガンマ分布の平均値である。また、上記の初期値は、γα０ｉ、α０ｉ−（−はαの真上に配置）、γβ０、β０−（−はαの真上に配置）であるとし、各初期値は予め与えられている、とする。
【数１６】
【００９７】
数式１６において、Ｔは、サンプル数である。
【数１７】
【数１８】
【数１９】
【００９８】
そして、推定手段２８３は、第一ステップと第二ステップとを繰り返し、第二ステップにおいて、αｉ−（−はαの真上に配置）、γαｉ、β−（−はβの真上に配置）、γβの値の変化が所定の条件を満たすほど小さい場合（ｎ回目の値と（ｎ−１）回目の値との差が閾値以内など）、収束したとして、処理を終了する。なお、最後は、第一ステップを行って終了する。そして、以下の数式２０におけるＸが取得される。
【数２０】
【００９９】
上記において、各行列は以下のような次元数を持つ。Ｙは（Ｎ＊Ｌ）×Ｔ次元行列、Ｇｃは（Ｎ＊Ｌ）×Ｍｃ次元行列、Ｇｓは（Ｎ＊Ｌ）×Ｍｓ次元行列、ＸｃはＭｃ×Ｔ次元行列、ＸｓはＭｓ×Ｔ次元行列である。
（実験１）
【０１００】
次に、本脳活動情報出力システムの第一の実験結果について説明する。本実験において、脳活動情報出力システムを構成する近赤外分光計測装置１は、２０のチャネルを有する。チャネルとは、送光プローブ１１と受光プローブ１２のすべての組み合わせである。つまり、本実験において、図６に示すように、近赤外分光計測装置１は５つの送光プローブ１１と４つの受光プローブ１２を有する。また、図７に示すように、各送光プローブ１１の識別子と各受光プローブ１２の識別子との組み合わせが２０セット存在する。なお、図６において、黒丸が送光プローブ１１、白丸が受光プローブ１２である。また、本実験において、最近傍のプローブの間隔はここでは１５ｍｍとした。また、ボクセルの大きさは４ｍｍ、縦１１ボクセル、横１１ボクセル、深さ５ボクセルとした。さらに、深さ１は頭皮、深さ２・３は頭蓋骨・脳脊髄液、深さ４・５は大脳皮質であるとした。
【０１０１】
そして、以下のようなシミュレーション実験を行った。つまり、シミュレーション実験において、被験者の脳活動があまり表れない安静状態で、近赤外分光計測装置１の送光プローブ１１は、波長７８０ｎｍと８３０ｎｍで発光した、と仮定した。そして、受光プローブ１２が受光する光量の対数値を、各波長において、チャネルごとに、２０回分、取得した。そして、脳活動情報出力装置２の安静時光強度取得部２４は、受光プローブ１２が受光する光量の対数値の平均と分散共分散行列を算出した。なお、図８は、波長７８０ｎｍにおける前記光量の対数値の平均値である。なお、図示しないが、波長８３０ｎｍにおける前記光量の対数値の平均も、同様に、取得した。図９は２波長のデータをまとめた、（チャネル数）×（波長数）＝４０次元を、複数回（今回は２０回）計測したデータから計算した分散共分散行列の一部（波長７８０nmに関する部分）である。そして、安静時光強度蓄積部２５は、安静時光強度情報格納部２３に、波長７８０ｎｍと波長８３０ｎｍにおける光量の対数値の平均および分散共分散行列を蓄積した。
【０１０２】
次に、脳活動情報出力装置２の脳活動時光強度情報取得部２６は、被験者の脳活動がある状態で、波長７８０ｎｍと波長８３０ｎｍにおける光量の対数値を、チャネルごとに取得した。取得したサンプル数は、今回は１回である。この光量の対数値は、図１０の「ｔｒｉａｌ１」の属性値である。なお、被験者の脳活動がある状態では、頭皮の血流変化もある。
【０１０３】
次に、光強度変化情報取得部２７は、波長７８０ｎｍと波長８３０ｎｍの各波長において、チャネルごと光強度変化情報（数式４のＹ）を算出した。つまり、光強度変化情報取得部２７は、数式４に従って、波長ごと、チャネルごと、およびサンプルごとに、光量の対数値の平均値から、脳活動時光強度情報取得部２６が取得した光量の対数値を減算し、光強度変化情報（Ｙ）を算出した。その光強度変化情報（Ｙ）が、図１０の属性値（Ｙ）である。なお、波長780nmにおける観測値Yに関して、Y(1:20)/max(abs(Y))を図示したのが、図１１（ａ）である。また、図１１（ａ）（ｂ）において、線の明るさは、その線の端にある送光プローブと受光プローブのセットからなるチャネルにおける光強度変化情報を表している。図１１（ａ）において、広がった頭皮血流アーチファクトによって、全体的に値が出ていることがわかる。しかし、図１１（ａ）を注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。また、波長830nmにおける観測値Yに関して、Y(21:40)/max(abs(Y))を図示したのが、図１１（ｂ）である。図１１（ｂ）において、広がった頭皮血流アーチファクトによって、全体的に値が出ている。しかし、図１１（ｂ）を注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。なお、780nmよりも強く値が出ているが、これはシミュレーション設定によるものである。
【０１０４】
次に、推定手段２８３は、感度情報格納手段２８１からヘモグロビン濃度変化に関する感度行列（上述したＧｃｏｘｙ，Ｇｃｄｅｏｘｙ，Ｇｓｏｘｙ，Ｇｓｄｅｏｘｙ）を読み出す。次に、推定手段２８３は、大脳皮質活動モデル情報格納部２１から大脳皮質活動モデル情報（例えば、局所的フィルタＷｃ）を読み出す。また、推定手段２８３は、頭皮血流変化モデル情報格納部２２から頭皮血流変化モデル情報（例えば、広域的フィルタＷｓ）を読み出す。
【０１０５】
次に、推定手段２８３は、以下のように変分ベイズアルゴリズムを用いて、Ｘを推定する。つまり、推定手段２８３は、予め格納している初期値を読み出す。初期値は、ここでは、すべてのボクセルに対して「γα０ｉ＝０」、かつ、すべてのボクセルに対して「α０ｉ−＝１００００」、かつ、「γβ０＝０」、かつ、「β０−＝１００００」である。
【０１０６】
そして、推定手段２８３は、上述したステップ１とステップ２の繰り返し処理を１０００回実行し、Ｘを算出した。推定手段２８３が算出したＸの要素であり、頭皮（深さ１）のオキシヘモグロビンの濃度変化（Xoxy）の推定値の集合を図１２に示す。図１２において、各セルの値は、各ボクセルに対応する。したがって、図１２において、１１×１１の値が存在する。また、数値は行（下方向）が頭部モデルの横方向、列（右方向）が頭部モデルの縦方向と対応する。かかることは、図１３、図１４等と同様である。また、推定手段２８３が算出したＸの要素であり、大脳皮質（深さ４）のオキシヘモグロビンの濃度変化（Xoxy）の推定値の集合を図１３に示す。なお、図示しないが、大脳皮質（深さ５）のオキシヘモグロビンの濃度変化（Xoxy）についても、推定手段２８３は、同様に算出した。また、推定手段２８３が算出したＸの要素であり、頭皮（深さ１）のデオキシヘモグロビンの濃度変化（Xdeoxy）の推定値の集合を図１４に示す。さらに、推定手段２８３が算出したＸの要素であり、大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化（Xdeoxy）の推定値の集合を図１５に示す。なお、図示しないが、大脳皮質（深さ５）のデオキシヘモグロビンの濃度変化（Xoxy）についても、推定手段２８３は、同様に算出した。
【０１０７】
次に、出力部２９は、推定部２８が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。出力部２９の出力例を図１６、図１７に示す。図１６（ａ）は、頭皮（深さ１）のオキシヘモグロビンの濃度変化の出力例である。図１６（ａ）等において、縦軸は頭部モデルの縦方向、横軸は頭部モデルの横方向を示し、その数字の単位はmm(ミリメーター)である。また、白黒の濃淡で表されるＸの値の単位はmM(ミリモーラー)である。図１６（ｂ）は、大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の出力例である。また、図１７（ａ）は、頭皮（深さ１）のデオキシヘモグロビンの濃度変化の出力例である。図１７（ｂ）は、大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の出力例である。
【０１０８】
さらに、シミュレーションで仮定した真のＸの値が図１８から図２１である。図１８は、シミュレーションで仮定した真の頭皮（深さ１）のオキシヘモグロビンの濃度変化の値の集合である。また、図１９は、シミュレーションで仮定した真の大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の値の集合である。また、図２０は、シミュレーションで仮定した真の頭皮（深さ１）のデオキシヘモグロビンの濃度変化の値の集合である。図２１は、シミュレーションで仮定した真の大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の値の集合である。
【０１０９】
シミュレーションで仮定した真の値と、推定部２８が取得した推定値とを比較すると、推定部２８は、脳活動をほぼ正確に推定できていることが分かる。
【０１１０】
次に、上記で説明した脳活動情報出力装置２の推定方法の概念について、図２２を用いて説明する。脳活動情報出力装置２における推定は、推定する値（Ｘ）の数より少ない観測値（Ｙ）を用いて、値（Ｘ）を推定しなければならない。これは逆問題の中でも不良設定問題と呼ばれ、正則化と呼ばれるような、なんらかの制限を入れなければ解くことはできない。今回用いる推定法ではベイズ推定を用いる。ベイズ推定ではその正則化と呼ばれるような制限を事前分布において確率的に組み込む。今回の推定法では、その事前分布において生理学的知見を反映させたような、大脳皮質と頭皮において異なる正則化とモデル化を行う。また、今回、ベイズ推定を現実的に行う方法の一例として、変分ベイズ法を用いる。また、値（Ｘ）を推定するために、脳活動情報出力装置２は、大脳皮質活動モデル情報として、スパース化を含む正則化のプログラムおよび局所的フィルタ（Ｗｃ）を保持しており、かつ、頭皮血流変化モデル情報として、スパース化を含まない正則化のプログラムおよび広域的フィルタ（Ｗｓ）を保持している。
【０１１１】
図２２において、Ｙ１からＹＬが観測値である。図２２の（４）における観測値YとXの関係は、光拡散過程においてRytov近似により吸収係数の差δμａが十分小さければ数式２１であらわされる。
【数２１】
数式２１を、数式３、数式５、数式６を用いて書き換えると数式２２となる。しかし、数式２２は観測ノイズを考慮していない関係式である。
【数２２】
【０１１２】
推定部２８は、図２２の（４）において、観測値ＹはＧＸにノイズσ−１Σｙがのって出てくると仮定する。この仮定は、脳活動時のノイズが安静時のノイズの定数倍でうまく近似できるという考えに基づく。
【数２３】
【０１１３】
なお、数２３において、Nは正規分布をあらわす。例えば、N(x|μ,Σｘ)はxが平均μ、分散共分散行列Σｘの正規分布に従うことをあらわす。
【０１１４】
そして、推定部２８は、図２２の（３）において、生理学的知見により、大脳皮質には局所的フィルタ、頭皮は広域的フィルタをかける。また、フィルタをかける前のＸの値をＺとおく。かかる処理の具体的な演算を数式２４に示す。
【数２４】
【０１１５】
Ｚの値に制限をかけないと、逆問題は解くことはできない。脳活動情報出力装置２は、以下のように、確率的に条件をつける。推定部２８は、図２２の（２）において、頭皮においては、頭皮血流変化モデル情報に従って、全部のボクセルで共通な正則化を行う。ここでいう正則化とは、大きなＺの値（つまり大きな濃度変化）についてペナルティをかけることである。これは、図２２の（３）の広域的フィルタ（Ｗｓ）と合わせて、空間的に広がった活動を意味する。かかる処理の具体的な演算を数式２５に示す。
【数２５】
【０１１６】
数式２５において、σは脳活動時の分散共分散行列の安静時に比べての割合の逆数である。また、Ｔはサンプル数である。βは頭皮の正則化に関するパラメータである。Iは単位行列である。
【０１１７】
また、推定部２８は、大脳皮質においては、各ボクセル毎に異なった正則化を行う。ここでは、少数のボクセルのＺの値が大きくなるように選ばれるようになっている（これをスパース化という。）。これは、図２２の（３）の局所的フィルタ（Ｗｃ）と合わせて、局在した活動を意味する。かかる処理の具体的な演算を数式２６に示す。
【数２６】
【０１１８】
数式２６において、Ａはαｉを対角要素に持つ行列であり、パラメータαｉは大脳皮質の各ボクセル毎に存在する正則化に関するパラメータである。
【０１１９】
また、パラメータαｉ、β、σの分布は、数式２７のように仮定する、とする。
【数２７】
【０１２０】
なお、数式２７において、Gammaはガンマ分布を示す。
【０１２１】
以上のように、推定部２８は、確率モデル化した。以下、ベイズ推定により、観測値Yから、最適なσ、Ａ、βの値と、Ｘの推定値を求める。
【数２８】
【０１２２】
数式２８はベイズの定理である。例えば、左辺の事後分布を最大とするようなＸ、σ、Ａ、β、またはＸ、σ、Ａ、βの期待値を推定値とすればよい。しかしながら、直接計算することは難しいため、今回の推定法は変分ベイズ法を用いる。変分ベイズ法では、事後分布を直接計算する代わりに、試験事後分布Qを用意し、真の事後分布に近似させる。この近似は数式２９で定義される自由エネルギーFの最大化によって行われる。
【数２９】
【０１２３】
今回は試験分布におけるＸ、σとＡ、βの間に数３０で示すように独立分解仮定をおく。これにより、計算を簡単に行うことが可能となる。
【数３０】
【０１２４】
この仮定により自由エネルギーの最大化は、Q1に関してFを最大化する「Ｘ、σ最大化ステップ」と、Q2に関してFを最大化する「Ａ、β最大化ステップ」を交互に繰り返すことにより実現できる。「Ｘ、σ最大化ステップ」に対応するアルゴリズムが、前述した第一ステップである。また、「Ａ、β最大化ステップ」に対応するアルゴリズムが、前述した第二ステップである。
【０１２５】
値が収束するまで繰り返した後、今回の推定においては、試験事後分布QにおけるＸ、σ、Ａ、βの期待値を推定値とした。
（実験２）
【０１２６】
次に、本脳活動情報出力システムの第二の実験結果について説明する。本実験において、最近傍のプローブの間隔は２０ｍｍとした。実験２における他の条件は、実験１と同じである。
【０１２７】
そして、以下のようなシミュレーション実験を行った。つまり、シミュレーション実験において、被験者の脳活動があまり表れない安静状態で、近赤外分光計測装置１の送光プローブ１１は、波長７８０ｎｍと８３０ｎｍで発光した、と仮定した。そして、受光プローブ１２が受光する光量の対数値を、各波長において、チャネルごとに、２０回分、取得した。そして、脳活動情報出力装置２の安静時光強度取得部２４は、受光プローブ１２が受光する光量の対数値の平均と分散共分散行列を算出した。なお、図２５は、波長７８０ｎｍにおける前記光量の対数値の平均値である。なお、図示しないが、波長８３０ｎｍにおける前記光量の対数値の平均も、同様に、取得した。図２６は２波長のデータをまとめた、（チャネル数）×（波長数）＝４０次元を、複数回（今回は２０回）計測したデータから計算した分散共分散行列の一部（波長７８０nmに関する部分）である。そして、安静時光強度蓄積部２５は、安静時光強度情報格納部２３に、波長７８０ｎｍと波長８３０ｎｍにおける光量の対数値の平均および分散共分散行列を蓄積した。
【０１２８】
次に、脳活動情報出力装置２の脳活動時光強度情報取得部２６は、被験者の脳活動がある状態で、波長７８０ｎｍと波長８３０ｎｍにおける光量の対数値を、チャネルごとに取得した。取得したサンプル数は、今回は１回である。この光量の対数値は、図２７の「ｔｒｉａｌ１」の属性値である。なお、被験者の脳活動がある状態では、頭皮の血流変化もある。
【０１２９】
次に、光強度変化情報取得部２７は、波長７８０ｎｍと波長８３０ｎｍの各波長において、チャネルごと光強度変化情報（数式４のＹ）を算出した。つまり、光強度変化情報取得部２７は、数式４に従って、波長ごと、チャネルごと、およびサンプルごとに、光量の対数値の平均値から、脳活動時光強度情報取得部２６が取得した光量の対数値を減算し、光強度変化情報（Ｙ）を算出した。その光強度変化情報（Ｙ）が、図２７の属性値（Ｙ）である。なお、波長780nmにおける観測値Yに関して、Y(1:20)/max(abs(Y))を図示したのが、図２８（ａ）である。また、図２８（ａ）（ｂ）において、線の明るさは、その線の端にある送光プローブと受光プローブのセットからなるチャネルにおける光強度変化情報を表している。図２８（ａ）において、広がった頭皮血流アーチファクトによって、全体的に値が出ていることがわかる。しかし、図２８（ａ）を注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。また、波長830nmにおける観測値Yに関して、Y(21:40)/max(abs(Y))を図示したのが、図２８（ｂ）である。図２８（ｂ）において、広がった頭皮血流アーチファクトによって、全体的に値が出ている。しかし、図２８（ｂ）を注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。なお、780nmよりも強く値が出ているが、これはシミュレーション設定によるものである。
【０１３０】
次に、推定手段２８３は、感度情報格納手段２８１からヘモグロビン濃度変化に関する感度行列（上述したＧｃｏｘｙ，Ｇｃｄｅｏｘｙ，Ｇｓｏｘｙ，Ｇｓｄｅｏｘｙ）を読み出す。次に、推定手段２８３は、大脳皮質活動モデル情報格納部２１から大脳皮質活動モデル情報（例えば、局所的フィルタＷｃ）を読み出す。また、推定手段２８３は、頭皮血流変化モデル情報格納部２２から頭皮血流変化モデル情報（例えば、広域的フィルタＷｓ）を読み出す。
【０１３１】
次に、推定手段２８３は、変分ベイズアルゴリズムを用いて、Ｘを推定する。つまり、推定手段２８３は、予め格納している初期値を読み出す。初期値は、ここでは、すべてのボクセルに対して「γα０ｉ＝０」、かつ、すべてのボクセルに対して「α０ｉ−＝１００００」、かつ、「γβ０＝０」、かつ、「β０−＝１００００」である。
【０１３２】
そして、推定手段２８３は、上述したステップ１とステップ２の繰り返し処理を１０００回実行し、Ｘを算出した。推定手段２８３が算出したＸの要素であり、頭皮（深さ１）のオキシヘモグロビンの濃度変化（Xoxy）の推定値の集合を図２９に示す。図２９において、各セルの値は、各ボクセルに対応する。したがって、図２９において、１１×１１の値が存在する。また、数値は行（下方向）が頭部モデルの横方向、列（右方向）が頭部モデルの縦方向と対応する。かかることは、図３０、図３１等と同様である。また、推定手段２８３が算出したＸの要素であり、大脳皮質（深さ４）のオキシヘモグロビンの濃度変化（Xoxy）の推定値の集合を図３０に示す。なお、図示しないが、大脳皮質（深さ５）のオキシヘモグロビンの濃度変化（Xoxy）についても、推定手段２８３は、同様に算出した。また、推定手段２８３が算出したＸの要素であり、頭皮（深さ１）のデオキシヘモグロビンの濃度変化（Xdeoxy）の推定値の集合を図３１に示す。さらに、推定手段２８３が算出したＸの要素であり、大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化（Xdeoxy）の推定値の集合を図３２に示す。なお、図示しないが、大脳皮質（深さ５）のデオキシヘモグロビンの濃度変化（Xoxy）についても、推定手段２８３は、同様に算出した。
【０１３３】
次に、出力部２９は、推定部２８が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。出力部２９の出力例を図３３、図３４に示す。図３３（ａ）は、頭皮（深さ１）のオキシヘモグロビンの濃度変化の出力例である。図３３（ａ）等において、縦軸は頭部モデルの縦方向、横軸は頭部モデルの横方向を示し、その数字の単位はmm(ミリメーター)である。また、白黒の濃淡で表されるＸの値の単位はmM(ミリモーラー)である。図３３（ｂ）は、大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の出力例である。また、図３４（ａ）は、頭皮（深さ１）のデオキシヘモグロビンの濃度変化の出力例である。図３４（ｂ）は、大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の出力例である。
【０１３４】
以上、本実施の形態によれば、精度高く脳活動を推定できる。また、変分ベイズ法を用いることにより、極めて精度高く脳活動を推定できる。
【０１３５】
なお、本実施の形態の実験によれば、７８０ｎｍと８３０ｎｍの２波長を用いた。しかし、異なる波長を２波長以上用いれば良い。
【０１３６】
また、本実施の形態において、頭皮血流変化と脳活動の分離は、高密度プローブ計測によってはじめて可能となる。短い距離のチャネルは浅い頭皮の活動を拾いやすく、長い距離のチャネルは深い脳活動の活動をも拾うからである。つまり、本実施の形態において、プローブの配置密度が１．５ｃｍ以下であることは極めて好適である。また、プローブの配置密度が１．７ｃｍ以下であることは好適である。なお、プローブの配置密度は２ｃｍ以下、または３ｃｍ未満等でも良い。
【０１３７】
また、本実施の形態によれば、頭部が半無限一層モデルであると仮定し、計測データから脳活動を推定した。しかし、脳活動情報出力システムにおいて、ＭＲＩ構造画像などを用いて、頭部が、頭皮・頭蓋骨・脳脊髄液・灰白質などからなる多層モデルであると仮定し、計測データから脳活動を推定しても良い。
【０１３８】
また、本実施の形態によれば、推定部２８が行う推定法は、主として、変分ベイズ法であった。ただし、推定部２８は、逆問題を解くための他の方法を用いても良い。例えば、推定部２８は、数式３１を格納しており、当該数式３１を用いて、Ｘを推定しても良い。
【数３１】
【０１３９】
つまり、推定部２８は、例えば、数式３１のようなコスト関数Ｃを最小にするようなＸを算出すればよい。なお、数式３１において、１項目が誤差項、２項目以降が正則化項と呼ばれる項であり、２項目以降が大きな値にペナルティを与える項になっている。
【０１４０】
数式３１において、大脳皮質Ｘｃには値（ノルム）の１乗、つまり「Ｌ１ノルム」、頭皮Ｘｓには値（ノルム）の２乗、つまり「Ｌ２ノルム」の正則化を与えている。後者のＬ２ノルムでの正則化は、我々が図２２の（２）で頭皮活動に課した正規分布の条件とほぼ同じ働きをする。また、一般に広がった推定値が得られやすいことが知られている。それに対し、前者のＬ１ノルムでの正則化を与えると、スパース（局在した）な推定値が得られることが知られている。
【０１４１】
なお、数式３１において、フィルタＷはかけていないが、上のλのパラメータや、正則化のノルムを０．９乗にする等の細かい値の修正をすれば、フィルタＷをかけた結果と似たようなことが再現できる。
【０１４２】
以上から、数式３１を用いても、変分ベイズ法を用いた場合と似たようなＸの推定は可能である、と言える。
【０１４３】
さらに、本実施の形態における処理は、ソフトウェアで実現しても良い。そして、このソフトウェアをソフトウェアダウンロード等により配布しても良い。また、このソフトウェアをＣＤ−ＲＯＭなどの記録媒体に記録して流布しても良い。なお、このことは、本明細書における他の実施の形態においても該当する。なお、本実施の形態における脳活動情報出力装置を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、記憶媒体に、局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報と、広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報と、送光プローブと受光プローブの組であるプローブセットを１組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納しており、コンピュータを、前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、前記記憶媒体の安静時光強度情報と、前記脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、前記推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部として機能させるためのプログラム、である。
【０１４４】
また、上記プログラムにおいて、前記推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式を格納し得る演算式格納手段と、前記演算式に、前記感度情報と前記光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する推定手段とを具備するものとして、コンピュータを機能させることは好適である。
【０１４５】
また、上記プログラムにおいて、前記推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式であり、コスト関数である演算式を格納し得る演算式格納手段と、前記演算式に、前記感度情報と前記光強度変化情報を代入し、当該演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定手段とを具備するものとして、コンピュータを機能させることは好適である。
【０１４６】
また、図２３は、本明細書で述べたプログラムを実行して、上述した実施の形態の脳活動情報出力装置等を実現するコンピュータの外観を示す。上述の実施の形態は、コンピュータハードウェア及びその上で実行されるコンピュータプログラムで実現され得る。図２３は、このコンピュータシステム３４０の概観図であり、図２４は、コンピュータシステム３４０のブロック図である。
【０１４７】
図２３において、コンピュータシステム３４０は、ＦＤドライブ、ＣＤ−ＲＯＭドライブを含むコンピュータ３４１と、キーボード３４２と、マウス３４３と、モニタ３４４とを含む。
【０１４８】
図２４において、コンピュータ３４１は、ＦＤドライブ３４１１、ＣＤ−ＲＯＭドライブ３４１２に加えて、ＭＰＵ３４１３と、ＣＤ−ＲＯＭドライブ３４１２及びＦＤドライブ３４１１に接続されたバス３４１４と、ブートアッププログラム等のプログラムを記憶するためのＲＯＭ３４１５とに接続され、アプリケーションプログラムの命令を一時的に記憶するとともに一時記憶空間を提供するためのＲＡＭ３４１６と、アプリケーションプログラム、システムプログラム、及びデータを記憶するためのハードディスク３４１７とを含む。ここでは、図示しないが、コンピュータ３４１は、さらに、ＬＡＮへの接続を提供するネットワークカードを含んでも良い。
【０１４９】
コンピュータシステム３４０に、上述した実施の形態の脳活動情報出力装置等の機能を実行させるプログラムは、ＣＤ−ＲＯＭ３５０１、またはＦＤ３５０２に記憶されて、ＣＤ−ＲＯＭドライブ３４１２またはＦＤドライブ３４１１に挿入され、さらにハードディスク３４１７に転送されても良い。これに代えて、プログラムは、図示しないネットワークを介してコンピュータ３４１に送信され、ハードディスク３４１７に記憶されても良い。プログラムは実行の際にＲＡＭ３４１６にロードされる。プログラムは、ＣＤ−ＲＯＭ３５０１、ＦＤ３５０２またはネットワークから直接、ロードされても良い。
【０１５０】
プログラムは、コンピュータ３４１に、上述した実施の形態の脳活動情報出力装置等の機能を実行させるオペレーティングシステム（ＯＳ）、またはサードパーティープログラム等は、必ずしも含まなくても良い。プログラムは、制御された態様で適切な機能（モジュール）を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいれば良い。コンピュータシステム３４０がどのように動作するかは周知であり、詳細な説明は省略する。
【０１５１】
なお、上記プログラムにおいて、ハードウェアによって行われる処理、は含まれない。
【０１５２】
また、上記プログラムを実行するコンピュータは、単数であってもよく、複数であってもよい。すなわち、集中処理を行ってもよく、あるいは分散処理を行ってもよい。
【０１５３】
また、上記各実施の形態において、各処理（各機能）は、単一の装置（システム）によって集中処理されることによって実現されてもよく、あるいは、複数の装置によって分散処理されることによって実現されてもよい。
【０１５４】
本発明は、以上の実施の形態に限定されることなく、種々の変更が可能であり、それらも本発明の範囲内に包含されるものであることは言うまでもない。
【産業上の利用可能性】
【０１５５】
以上のように、本発明にかかる脳活動情報出力装置は、精度高く脳活動を推定できる、という効果を有し、脳活動推定装置等として有用である。
【符号の説明】
【０１５６】
１ 近赤外分光計測装置
２ 脳活動情報出力装置
１１ 送光プローブ
１２ 受光プローブ
１３ 受付部
１４ 制御部
２１ 大脳皮質活動モデル情報格納部
２２ 頭皮血流変化モデル情報格納部
２３ 安静時光強度情報格納部
２４ 安静時光強度取得部
２５ 安静時光強度蓄積部
２６ 脳活動時光強度情報取得部
２７ 光強度変化情報取得部
２８ 推定部
２９ 出力部
２８１ 感度情報格納手段
２８２ 演算式格納手段
２８３ 推定手段
【技術分野】
【０００１】
本発明は、脳の活動に関する情報を取得し、出力する脳活動情報出力装置等に関するものである。
【背景技術】
【０００２】
従来、ＮＩＲＳ（Ｎｅａｒ−Ｉｎｆｒａｒｅｄ Ｓｐｅｃｔｒｏｓｃｏｐｙ）データから推定される脳活動の表示技術があった。かかる技術は、計測した光強度変化をヘモグロビン濃度に変換し、２次元のトポグラフィとして表示する方法であった。
【０００３】
また、近年は、拡散光トモグラフィ（ＤＯＴ）を用いて推定したものを３次元のトモグラフィとして表示するやり方が広まってきている（非特許文献１、非特許文献２参照）。
【０００４】
また、これまでのＤＯＴでは、波長毎に単独で吸収係数の変化を推定し、それを変換することで、ヘモグロビンの濃度を推定していた（非特許文献２参照）。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】D. A. Boas and A. M. Dale, "Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function", Applied Optics 44 (2005) pp. 1957-1968.
【非特許文献２】B. W. Zeff, B. R. White, H. Dehghani, B. L. Schlaggar,and J. P. Culver, "Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography",PNAS 104 (2007) pp. 12169-12174.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、従来の技術においては、大脳皮質における血行動態反応は局所的であり、頭皮における血行動態反応は広域的であるという知識をうまく活用することができず、頭皮血流に起因するアーチファクトを適切に除去できないために、精度高く脳活動を推定することができなかった。また、従来の技術においては、波長毎の観測値からそれぞれ別々に吸収係数の変化を推定し、それを変換することで、ヘモグロビンの濃度を推定していたため、精度高く脳活動を推定することができなかった。なお、ここでの脳活動の推定とは、神経活動に伴う脳血行動態反応を推定することである。
【課題を解決するための手段】
【０００７】
本第一の発明の脳活動情報出力装置は、局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報を格納している大脳皮質活動モデル情報格納部と、広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報を格納している頭皮血流変化モデル情報格納部と、送光プローブと受光プローブの組であるプローブセットを１組以上有する近赤外分光計測装置を用いて、各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納している安静時光強度情報格納部と、近赤外分光計測装置を用いて、各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、安静時光強度情報格納部の安静時光強度情報と、脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、複数の波長ごと、およびプローブセットごとの光強度変化情報を、大脳皮質活動モデル情報および頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部とを具備する脳活動情報出力装置である。
【０００８】
かかる構成により、精度高く脳活動を推定できる。
【０００９】
また、本第二の発明の脳活動情報出力装置は、第一の発明に対して、推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式を格納し得る演算式格納手段と、演算式に、感度情報と光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する推定手段とを具備する脳活動情報出力装置である。
【００１０】
かかる構成により、極めて精度高く脳活動を推定できる。
【００１１】
また、本第三の発明の脳活動情報出力装置は、第一の発明に対して、推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の情報の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式であり、コスト関数である演算式を格納し得る演算式格納手段と、演算式に、感度情報と光強度変化情報を代入し、演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定手段とを具備する脳活動情報出力装置である。
【００１２】
かかる構成により、精度高く脳活動を推定できる。
【発明の効果】
【００１３】
本発明による脳活動情報出力装置によれば、精度高く脳活動を推定できる。
【図面の簡単な説明】
【００１４】
【図１】実施の形態１における脳活動情報出力システムの概念図
【図２】同脳活動情報出力システムのブロック図
【図３】同脳活動情報出力装置の動作について説明するフローチャート
【図４】同変分ベイズ推定法について説明するフローチャート
【図５】同プローブの配置を説明する図
【図６】同実験におけるプローブの配置を説明する図
【図７】同実験におけるチャネルを説明する図
【図８】同波長７８０ｎｍにおける光量の計測値等を示す図
【図９】同共分散行列を示す図
【図１０】同実験における各波長の光強度変化情報を示す図
【図１１】同観測値をビジュアルに示した図
【図１２】同推定手段が算出した頭皮（深さ１）のオキシヘモグロビンの濃度変化の推定値を示す図
【図１３】同推定手段が算出した大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の推定値を示す図
【図１４】同推定手段が算出した頭皮（深さ１）のデオキシヘモグロビンの濃度変化の推定値を示す図
【図１５】同推定手段が算出した大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の推定値を示す図
【図１６】同出力部の出力例を示す図
【図１７】同出力部の出力例を示す図
【図１８】同シミュレーションで仮定した真の頭皮（深さ１）のオキシヘモグロビンの濃度変化の値を示す図
【図１９】同シミュレーションで仮定した真の大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の値を示す図
【図２０】同シミュレーションで仮定した真の頭皮（深さ１）のデオキシヘモグロビンの濃度変化の値を示す図
【図２１】同シミュレーションで仮定した真の大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の値を示す図
【図２２】同脳活動情報出力装置の推定方法の概念について説明図
【図２３】同コンピュータシステムの概観図
【図２４】同コンピュータシステムのブロック図
【図２５】同波長７８０ｎｍにおける光量の計測値等を示す図
【図２６】同共分散行列を示す図
【図２７】同実験における各波長の光強度変化情報を示す図
【図２８】同観測値をビジュアルに示した図
【図２９】同推定手段が算出した頭皮（深さ１）のオキシヘモグロビンの濃度変化の推定値を示す図
【図３０】同推定手段が算出した大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の推定値を示す図
【図３１】同推定手段が算出した頭皮（深さ１）のデオキシヘモグロビンの濃度変化の推定値を示す図
【図３２】同推定手段が算出した大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の推定値を示す図
【図３３】同出力部の出力例を示す図
【図３４】同出力部の出力例を示す図
【発明を実施するための形態】
【００１５】
以下、脳活動情報出力装置等の実施形態について図面を参照して説明する。なお、実施の形態において同じ符号を付した構成要素は同様の動作を行うので、再度の説明を省略する場合がある。
【００１６】
（実施の形態１）
本実施の形態において、脳活動を推定する脳活動情報出力システムについて説明する。また、脳活動情報出力システムは、近赤外分光計測装置（ＮＩＲＳ）を具備する。なお、脳活動の推定とは、神経活動に伴う脳血行動態反応を推定することである。
【００１７】
また、本実施の形態における脳活動情報出力システムにおいて、拡散光トモグラフィ（ＤＯＴ）を用いる。ＤＯＴとは、頭部内での光子移動を拡散過程としてモデル化することにより、脳内活動と観測光強度変化の関係を物理的に正しく求め、それを利用して脳活動の推定の精度を高めようというものである。
【００１８】
また、本実施の形態における脳活動情報出力システムにおいて、特に、局所的な活動である大脳皮質の活動、および広域的な活動である頭皮の血流変化をモデル化しており、かかるモデルを用いて、脳活動を推定する方法を採用する。また、特に、本実施の形態において、頭皮血流をＤＯＴのモデルに組み込んで、大脳皮質における神経活動と頭皮血流に起因するアーチファクトとを同時推定することにより、頭皮血流に起因するアーチファクトを除去し、脳活動を推定する階層変分ベイズ法について説明する。なお、後述する数式に出現する上付きのcは大脳皮質(cortex),sは頭皮(scalp)を表す。
【００１９】
さらに、具体的には、本実施の形態における脳活動情報出力システムにおいて、近赤外分光計測装置における観測結果から、神経活動に伴う大脳皮質内のオキシヘモグロビン（Ｏｘｙ Ｈｂ）、デオキシヘモグロビン（Ｄｅｏｘｙ Ｈｂ）の濃度変化の場所と大きさを推定する。
【００２０】
図１は、本実施の形態における脳活動情報出力システムの概念図である。脳活動情報出力システムは、近赤外分光計測装置１、および脳活動情報出力装置２を具備する。
【００２１】
また、図２は、本実施の形態における脳活動情報出力システムのブロック図である。
【００２２】
近赤外分光計測装置１は、送光プローブ１１、受光プローブ１２、受付部１３、および制御部１４を具備する。近赤外分光計測装置１は、送光プローブ１１と受光プローブ１２の組であるプローブセット（チャネルとも言う）を１組以上有する。
【００２３】
脳活動情報出力装置２は、大脳皮質活動モデル情報格納部２１、頭皮血流変化モデル情報格納部２２、安静時光強度情報格納部２３、安静時光強度取得部２４、安静時光強度蓄積部２５、脳活動時光強度情報取得部２６、光強度変化情報取得部２７、推定部２８、および出力部２９を具備する。
【００２４】
また、推定部２８は、感度情報格納手段２８１、演算式格納手段２８２、および推定手段２８３を具備する。
【００２５】
近赤外分光計測装置１を構成する送光プローブ１１は、近赤外光を生体（ここでは、特に頭部）に照射する。
【００２６】
受光プローブ１２は、生体（ここでは、特に頭部）内部で吸収あるいは散乱され、生体内から出てくる光を受光する。なお、プローブの配置密度が３ｃｍ未満であることは好適である。また、プローブの配置密度が２ｃｍ以下であることはさらに好適である。また、プローブの配置密度が１．７ｃｍ以下であることは好適である。さらに、プローブの配置密度が１．５ｃｍ以下であることは極めて好適である。プローブの配置密度とは、送光プローブ１１と受光プローブ１２間の間隔である。なお、プローブの配置密度は密であればある程良いのは、頭皮血流をとらえるための短い距離のチャネルがたくさんできること、および脳活動をとらえるための３〜４ｃｍ程度の中距離のチャネルがたくさんできること等からである。
【００２７】
受付部１３は、外部からの指示を受け付ける。かかる指示は、近赤外光を生体に照射する指示である照射指示などである。
【００２８】
制御部１４は、受付部１３が照射指示を受け付けた場合に、送光プローブ１１に近赤外光の照射を指示する。なお、送光プローブ１１は、通常、制御部１４から照射指示を受け付けた場合に、近赤外光を生体に照射する。
【００２９】
脳活動情報出力装置２は、近赤外分光計測装置１が計測した光強度に関する情報を用いて、脳活動に関する情報である脳活動情報を取得し、出力する装置である。ここで、出力とは、ディスプレイへの表示、プロジェクターを用いた投影、プリンタへの印字、音出力、外部の装置への送信、記録媒体への蓄積、他の処理装置や他のプログラムなどへの処理結果の引渡しなどを含む概念である。
【００３０】
脳活動情報出力装置２を構成する大脳皮質活動モデル情報格納部２１は、大脳皮質活動モデル情報を格納している。大脳皮質活動モデル情報は、局所的な活動である大脳皮質の活動のモデルを示す情報である。大脳皮質活動モデル情報は、例えば、後述する局所的フィルタＷＣ、後述する正則化の処理を示すプログラムやフラグ、大脳皮質の各ボクセルにおける正則化（大きな値にペナルティ）に関するガンマ分布のシェイプパラメータ（後述するγα０ｉ）、大脳皮質の各ボクセルにおける正則化に関するガンマ分布の平均値（後述するα０ｉ−）等である。なお、大脳皮質活動モデル情報に対応する正則化は、ボクセル毎に異なった正則化である。また、ここでいう正則化とは、大きな正則化対象の値についてペナルティをかけることである。
【００３１】
頭皮血流変化モデル情報格納部２２は、頭皮血流変化モデル情報を格納している。頭皮血流変化モデル情報は、広域的な活動である頭皮の血流変化のモデルを示す情報である。頭皮血流変化モデル情報は、例えば、後述する広域的フィルタＷＳ、後述する正則化の処理を示すプログラムやフラグ、頭皮の正則化に関するガンマ分布のシェイプパラメータ（後述するγβ０）、頭皮の正則化に関するガンマ分布の平均値（後述するβ０−）等である。なお、頭皮血流変化モデル情報に対応する正則化は、全部のボクセルで共通な正則化である。
【００３２】
安静時光強度情報格納部２３は、複数の各波長に対応付けて安静時光強度情報を格納している。安静時光強度情報は、送光プローブ１１と受光プローブ１２の組であるプローブセットを１組以上有する近赤外分光計測装置１を用いて、各プローブセットを構成する各送光プローブ１１から、一の光強度の光を発した場合に、安静状態において測定された各受光プローブ１２の光強度に関する情報である
【００３３】
安静時光強度取得部２４は、安静状態の被験者が装着している近赤外分光計測装置１の受光プローブ１２から取得した光強度に関する情報である安静時光強度情報を、近赤外分光計測装置１から取得する。安静時光強度取得部２４は、複数の各波長ごとに、安静時光強度情報を取得する。
【００３４】
安静時光強度蓄積部２５は、安静時光強度取得部２４が取得した安静時光強度情報を、複数の各波長に対応付けて安静時光強度情報格納部２３に蓄積する。
【００３５】
脳活動時光強度情報取得部２６は、複数の各波長ごとに脳活動時光強度情報を取得する。脳活動時光強度情報は、脳活動が行われている状態において測定された各受光プローブ１２の光強度に関する情報である。脳活動時光強度情報は、被験者の脳の活動が行われている状態において、近赤外分光計測装置１を用いて、各プローブセットを構成する各送光プローブ１１から、一の光強度の光を発した場合に、各受光プローブ１２が受光した光強度に関する情報である。
【００３６】
光強度変化情報取得部２７は、安静時光強度情報格納部２３の安静時光強度情報と、脳活動時光強度情報取得部２６が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度変化情報を取得する。光強度変化情報は、光強度の変化に関する情報である。光強度変化情報は、例えば、光強度の変化率、光強度の変化率の対数をとった値、光強度の変化量等である。また、光強度変化情報取得部２７は、光強度のゆらぎ（自発脳活動やセンサーの機器ノイズにより光強度に加わるノイズ）の情報Σyも取得しても良い。
【００３７】
推定部２８は、複数の波長ごと、およびプローブセットごとの光強度変化情報を、大脳皮質活動モデル情報および頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する。なお、推定部２８は、推定する値（後述するＸ）の数より少ない観測値（後述するＹ）を用いて、値（後述するＸ）を推定し、取得する。なお、かかる推定部２８が解く問題を逆問題という。推定部２８は、後述する変分ベイズ法を用いて逆問題を解くことは好適であるが、逆問題の解法は問わない。また、推定部２８は、頭部の領域ごとに、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する。頭部の領域ごととは、例えば、頭部のボクセルごと、または頭部を三角錐に区切った領域ごと等である。
【００３８】
感度情報格納手段２８１は、感度情報を格納し得る。感度情報とは、大脳皮質または頭皮の活動の変化（例えば、後述するＸが対応）と検出された光強度の変化の情報（例えば、後述するＹが対応）の関係を示す情報である。感度情報は、例えば、後述する感度行列である。感度情報は、推定部２８が算出し、感度情報格納手段２８１に、少なくとも一時的に格納しても良い。
【００３９】
演算式格納手段２８２は、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式（例えば、演算式１）であり、ベイズ推定の演算式を格納し得る。また、演算式格納手段２８２は、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式（例えば、演算式２）であり、コスト関数である演算式を格納し得る。なお、演算式１は、例えば、後述する数式２８である。また、演算式２は、例えば、後述する数式３１である。
【００４０】
推定手段２８３は、演算式（例えば、演算式１）に、感度情報と光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する。また、推定手段２８３は、演算式（例えば、演算式２）に、感度情報と光強度変化情報を代入し、演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得しても良い。なお、「コストが最小である」とは、コストに「−１」をかけた値が最大であることも含む。つまり、「コストが最小である」とは、実質的に考慮する。
【００４１】
出力部２９は、推定部２８が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。なお、出力部２９は、頭皮の血管のオキシヘモグロビンの濃度変化、および頭皮の血管のデオキシヘモグロビンの濃度変化も出力しても良い。また、出力部２９は、大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を、数値を用いて出力しても良いし、図的に出力しても良い。その出力態様は問わない。
【００４２】
大脳皮質活動モデル情報格納部２１、頭皮血流変化モデル情報格納部２２、安静時光強度情報格納部２３、感度情報格納手段２８１、および演算式格納手段２８２は、不揮発性の記録媒体が好適であるが、揮発性の記録媒体でも実現可能である。
【００４３】
大脳皮質活動モデル情報格納部２１等に大脳皮質活動モデル情報等が記憶される過程は問わない。例えば、記録媒体を介して大脳皮質活動モデル情報等が大脳皮質活動モデル情報格納部２１等で記憶されるようになってもよく、通信回線等を介して送信された大脳皮質活動モデル情報等が大脳皮質活動モデル情報格納部２１等で記憶されるようになってもよく、あるいは、入力デバイスを介して入力された大脳皮質活動モデル情報等が大脳皮質活動モデル情報格納部２１等で記憶されるようになってもよい。
【００４４】
安静時光強度取得部２４は、通常、受信手段等から実現され得る。
【００４５】
安静時光強度蓄積部２５、脳活動時光強度情報取得部２６、光強度変化情報取得部２７、推定部２８、および推定手段２８３は、通常、ＭＰＵやメモリ等から実現され得る。安静時光強度蓄積部２５等の処理手順は、通常、ソフトウェアで実現され、当該ソフトウェアはＲＯＭ等の記録媒体に記録されている。但し、ハードウェア（専用回路）で実現しても良い。
【００４６】
出力部２９は、ディスプレイやスピーカー等の出力デバイスを含むと考えても含まないと考えても良い。出力部２９は、出力デバイスのドライバーソフトまたは、出力デバイスのドライバーソフトと出力デバイス等で実現され得る。
【００４７】
次に、脳活動情報出力システムを構成する脳活動情報出力装置２の動作について、図３のフローチャートを用いて説明する。なお、図３のフローチャートの開始前において、安静時光強度取得部２４は、複数の波長ごとに、安静時光強度情報を取得し、安静時光強度蓄積部２５は、複数の波長ごとに、当該安静時光強度情報を安静時光強度情報格納部２３に蓄積した、とする。そして、安静時光強度情報格納部２３には、複数の各波長ごとに、安静時光強度情報が格納されている、とする。
【００４８】
（ステップＳ３０１）脳活動時光強度情報取得部２６は、複数の各波長ごとに脳活動時光強度情報を取得したか否かを判断する。複数の各波長の脳活動時光強度情報を取得すればステップＳ３０２に行き、取得しなければステップＳ３０１に戻る。
【００４９】
（ステップＳ３０２）光強度変化情報取得部２７は、安静時光強度情報格納部２３に格納されている安静時光強度情報を、複数の各波長に対応付けて読み出す。
【００５０】
（ステップＳ３０３）光強度変化情報取得部２７は、ステップＳ３０１で取得された各波長ごとの複数の脳活動時光強度情報と、ステップＳ３０２で読み出された波長ごとの複数の安静時光強度情報とから、波長ごとに、光強度変化情報を算出する。光強度変化情報取得部２７は、例えば、波長ごとに、安静時光強度情報から脳活動時光強度情報を減算することにより、光強度変化情報を取得する。
【００５１】
（ステップＳ３０４）推定手段２８３は、演算式格納手段２８２から演算式を読み出す。
【００５２】
（ステップＳ３０５）推定手段２８３は、感度情報格納手段２８１から感度情報を読み出す。なお、例えば、推定手段２８３は予め感度情報を算出し、感度情報格納手段２８１に格納している、とする。
【００５３】
（ステップＳ３０６）推定手段２８３は、大脳皮質活動モデル情報格納部２１から大脳皮質活動モデル情報を読み出す。
【００５４】
（ステップＳ３０７）推定手段２８３は、頭皮血流変化モデル情報格納部２２から頭皮血流変化モデル情報を読み出す。
【００５５】
（ステップＳ３０８）推定手段２８３は、演算式に、光強度変化情報、感度情報を代入し、かつ、大脳皮質活動モデル情報と頭皮血流変化モデル情報とを用いて、変分ベイズ推定法により、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する。かかる変分ベイズ推定法について、図４のフローチャートを用いて説明する。
【００５６】
（ステップＳ３０９）出力部２９は、推定部２８が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。ステップＳ３０１に戻る。
【００５７】
なお、図３のフローチャートにおいて、電源オフや処理終了の割り込みにより処理は終了する。
【００５８】
次に、ステップＳ３０８の変分ベイズ推定法について、図４のフローチャートを用いて説明する。
【００５９】
（ステップＳ４０１）推定手段２８３は、予め格納されている初期値を読み出す。読み出す初期値は、は、γα０ｉ、α０ｉ−（−はαの真上に配置）、γβ０、β０−（−はβの真上に配置）である。また、初期値は、ユーザが与えても良いし、他の装置から受信等しても良い。ここで、γα０ｉは大脳皮質の各ボクセルにおける正則化（大きな値にペナルティ）に関するガンマ分布のシェイプパラメータの初期値、α０ｉ−は大脳皮質の各ボクセルにおける正則化に関するガンマ分布の平均値の初期値、γβ０は頭皮の正則化に関するガンマ分布のシェイプパラメータの初期値、β０−は頭皮の正則化に関するガンマ分布の平均値の初期値である。
【００６０】
（ステップＳ４０２）推定手段２８３は、後述する第一ステップの演算を行う。第一ステップの演算は、変分ベイズ法における、「Ｘ，σ最大化ステップ」である。
【００６１】
（ステップＳ４０３）推定手段２８３は、後述する第二ステップの演算を行う。第二ステップの演算は、変分ベイズ法における、「Ａ，β最大化ステップ」である。
【００６２】
（ステップＳ４０４）推定手段２８３は、大脳皮質の各ボクセルにおける正則化に関するガンマ分布のシェイプパラメータ（γαｉ）、大脳皮質の各ボクセルにおける正則化に関するガンマ分布の平均値（αｉ−）、頭皮の正則化に関するガンマ分布のシェイプパラメータ（γβ）、および頭皮の正則化に関するガンマ分布の平均値（β−）の値が収束したか否かを判断する。収束すればステップＳ４０５に行き、収束しなければステップＳ４０２に戻る。なお、収束したか否かは、予め決められた条件を満たすほど、上記の４つの変数の各値の変化が小さくなったか否かにより判断される。
【００６３】
（ステップＳ４０５）推定手段２８３は、第一ステップの演算を行う。
【００６４】
（ステップＳ４０６）推定手段２８３は、推定結果を算出する演算を行う。処理を終了する。
【００６５】
以下、本実施の形態における脳活動情報出力システムの具体的な動作について説明する。脳活動情報出力システムの概念図は、図１である。
【００６６】
まず、本脳活動情報出力システムにおいて、頭部が半無限一層モデルであると仮定し、計測された光強度から脳活動を推定する場合について説明する。
【００６７】
また、プローブは、図５に示すように、１５ｍｍ間隔で頭部モデルの表面に配置されている、とする。図５において、黒丸が送光プローブ１１、白丸が受光プローブ１２を表す。また、図５において、直方体で囲ってある部分はこれから推定値を求めていく領域（頭部）で、一辺４ｍｍからなる立方体のボクセルに分けて考える。図５において、頭部は、縦１１ボクセル、横１１ボクセル、高さ５ボクセルに分割されている。ここで、深さ１を頭皮、深さ２・３を頭蓋骨・脳脊髄液、深さ４以降を大脳皮質とする。
【００６８】
また、今回の頭部モデルの光学特性は、ヒト灰白質と同じであり、換算散乱係数「μｓ'＝１．１１ｍｍ−１」、吸収係数「μａ＝０．０１８６ｍｍ−１」とする。さらに、表面境界における屈折率は「ｎ＝１．４」とする。
【００６９】
上記の条件で、被験者は安静な状態を保っており、かかる状態において、安静時光強度取得部２４は、複数の波長ごとに、安静時光強度情報を取得し、安静時光強度蓄積部２５は、複数の波長ごとに、当該安静時光強度情報を安静時光強度情報格納部２３に蓄積した、とする。そして、例えば、安静時光強度情報格納部２３には、波長＝７８０ｎｍ、８３０ｎｍ（以下、適宜、単位を省略する。）の２つの各波長ごとに、安静時光強度情報が格納された、とする。なお、安静時とは、脳活動や頭皮血流変化があまりない安静状態をいう。
【００７０】
そして、一般に、光強度は機器ノイズや自発脳活動によってゆらぐので、各波長ごと、各チャネルごとに多数計測し、安静時光強度取得部２４は、その平均値と分散共分散行列を算出することは好適である。ここで、光強度の平均値は、数式１のように示される。また、分散共分散行列はΣｙである。
【数１】
【００７１】
なお、数式１において、Φ０は安静時の光強度であり、ｒｓ→（ｒの真上に→が存在）は送光プローブ１１の位置、ｒｄ→（ｒの真上に→が存在）は受光プローブ１２の位置を示す。数式１のバーは、平均値を意味する。
【００７２】
次に、脳活動時光強度情報取得部２６は、各チャネルごと、波長＝７８０、８３０の２つの各波長ごとに、脳活動時光強度情報を取得した、とする。具体的には、脳活動時光強度情報取得部２６は、数式２に示す脳活動時光強度情報を、チャネルごと、波長ごとに取得する。また、脳活動時光強度は複数回計測して取得してもよい。その計測の番号をサンプル番号とする。
【数２】
【００７３】
次に、光強度変化情報取得部２７は、安静時光強度情報格納部２３に格納されている２つの安静時光強度情報を、各波長（波長＝７８０、８３０）に対応付けて読み出す。
【００７４】
次に、光強度変化情報取得部２７は、取得された波長ごとの２つの脳活動時光強度情報と、読み出された波長ごとの２つの安静時光強度情報とから、波長ごとに、光強度変化情報を、以下のように算出する。
【００７５】
つまり、光強度変化情報取得部２７は、ここでは、数式３の演算式を格納しており、それぞれのサンプル番号ごとに、かかる数式３に、安静時光強度情報（ｌｏｇΦ０）と脳活動時光強度情報（ｌｏｇΦ）を代入し、光強度変化情報（Ｙ）を取得する。
【数３】
【００７６】
ここでtはサンプル番号を表す。なお、ＹはＹ（ｔ）を各列に持つ行列であり、その次元は「（チャネル数×波長数）×（サンプル数）」である。今回は、サンプル数は１とする。従って、ここでは、例えば、Ｙは４０次元である、とする。
【００７７】
また、ここで、光強度変化情報取得部２７は、数式４を格納しており、数式４の演算式に、安静時光強度取得部２４が取得した安静時光強度の平均値と、脳活動時光強度情報（数式２）を代入し、光強度変化情報Ｙを算出することは好適である。
【数４】
【００７８】
そして、頭部モデルと、図５に示すようにプローブ配置（送光プローブ１１と受光プローブ１２の配置）が決まっているので、推定手段２８３は、光拡散過程を解くことにより、以下の数式５を用いて、各波長ごとの吸収係数変化に関する感度行列Ｇを計算して、取得する。ここで、感度行列Ｇは、感度情報の一例である。なお、数式５において、定義の仕方によっては定数倍がかかる場合がある。
【数５】
【００７９】
なお、数式５において、jはチャネル（送光プローブ１１と受光プローブ１２のペア）のインデックス、iは頭部ボクセルのインデックスである。また、数式５において、Φ０は安静時の光強度であり、Φ０（第一変数，第二変数）は、第一変数の位置から出発し、第二変数の位置に到達する光の強度を示す。なお、Φ０は頭部モデル内での光子の拡散過程を解くことで得られる。また、ｒｓ， j→（ｒの真上に→が存在）はチャネルjの送光プローブ１１の位置を、ｒｄ，j→（ｒの真上に→が存在）はチャネルjの受光プローブ１２の位置を示す。また、ｒi→（ｒの真上に→が存在）は頭部ボクセルiの位置を示す。
【００８０】
また、ここでは、推定手段２８３は、感度行列を、頭皮（図５の深さ１の領域）と大脳皮質（図５の深さ４または深さ５の領域）のみで求める。図５において、頭蓋骨・脳脊髄液（深さ２または深さ３の領域）にはヘモグロビンとオキシヘモグロビンが存在せず、常にＸ＝０だから、感度行列Ｇ＝０である。また、一般的に、頭部モデルの光学特性は波長ごとに違うので、推定手段２８３は、感度行列を波長毎に算出する。なお、Ｘは、ここでは、オキシヘモグロビンの濃度変化およびデオキシヘモグロビンの濃度変化である。
【００８１】
つまり、推定手段２８３は、Ｇｃ７８０、Ｇｓ７８０、Ｇｃ８３０、Ｇｓ８３０を算出し、感度情報格納手段２８１に蓄積しておく。ここで、上付きのｃは大脳皮質（深さ４または深さ５の領域）、ｓは頭皮（深さ１の領域）を意味する。下付きの７８０は波長７８０ｎｍ、８３０は波長８３０ｎｍを意味する。つまり、Ｇｃ７８０は、波長７８０ｎｍにおける大脳皮質の感度行列である。
【００８２】
次に、吸収係数変化に関する感度行列から、ヘモグロビン濃度変化に関する感度行列を算出する。その変換式は、取り扱う波長が（波長＝７８０、８３０）の２つの場合、以下の数式６に示すように詳細に定義される。
【数６】
【００８３】
数式６において、δμａは脳活動時の吸収係数と安静時の吸収係数の差である。また、δμａは波長毎に値が異なる。つまり、δμａ，７８０は、波長７８０ｎｍにおける、脳活動時の吸収係数と安静時の吸収係数の差である。また、数式６において、εは、モル吸光係数であり、波長や物質（オキシヘモグロビン、デオキシヘモグロビン）により異なる。なお、数式８は、波長ごと、物質ごとのモル吸光係数の値を示す。
【００８４】
また、それぞれの波長における吸収係数の変化は、オキシヘモグロビンの濃度変化Ｘｏｘｙと、デオキシヘモグロビンの濃度変化Ｘｄｅｏｘｙとを用いて、数式７のように表される。そして、数式６の（１）は、数式７を用いて、数式６の（２）に変形される。そして、数式６の（２）は、計算により、数式６の（３）に変形される。そして、ヘモグロビン濃度変化に関する感度行列Ｇは（Ｇｃｏｘｙ，Ｇｃｄｅｏｘｙ，Ｇｓｏｘｙ，Ｇｓｄｅｏｘｙ）である。
【数７】
【数８】
【００８５】
次に、推定手段２８３は、感度情報格納手段２８１から感度行列Ｇ（Ｇｃｏｘｙ，Ｇｃｄｅｏｘｙ，Ｇｓｏｘｙ，Ｇｓｄｅｏｘｙ）を読み出す。なお、Ｇｃ７８０などは吸収係数変化に関する感度行列、Ｇｃｏｘｙなどは、ヘモグロビン濃度変化に関する感度行列である。また、吸収係数変化に関する感度行列とヘモグロビン濃度変化に関する感度行列とは、上述したように、数式６の定義式により簡単に変換できる。なお、数式５を用いて、吸収係数変化に関する感度行列が算出され、次に、当該感度行列にモル吸光係数などをかけて、ヘモグロビン濃度変化に関する感度行列が算出される。そして、ヘモグロビン濃度変化に関する感度行列が上記の推定に用いられる。
【００８６】
次に、推定手段２８３は、大脳皮質活動モデル情報格納部２１から大脳皮質活動モデル情報（例えば、局所的フィルタＷＣ等）を読み出す。また、推定手段２８３は、頭皮血流変化モデル情報格納部２２から頭皮血流変化モデル情報（例えば、広域的フィルタＷＳ等）を読み出す。
【００８７】
次に、推定手段２８３は、以下のように２つのステップ（第一ステップ、第二ステップ）により、変分ベイズ推定を用いて、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を推定する。なお、推定手段２８３は、後述する演算式等を、通常、予め格納している、とする。
【００８８】
まず、第一ステップとは、以下の数式９から数式１５を用いて、Σ、σ―（−はσの上に配置）、γσ、Δ、＜（ΔＸ）２＞を算出する処理である。ここで、Σは、数式９により定義される。また、Σｙは安静時の観測値Ｙの分散共分散行列である。また、σは脳活動時の分散共分散行列の安静時に比べての割合の逆数である。つまり、脳活動時の観測値Yの分散共分散行列はσ−１Σｙであるとする。なお、初期値として、「Ｐ０（σ）∝１／σ」のような分布に従うと仮定するが、アルゴリズム途中ではガンマ分布にしたがう。その平均値がσ―である。また、シェイプパラメータをγσとする。なお、数式１０のγσは、自由エネルギーが最大となるような値である。
【数９】
【００８９】
なお、行列の上付けのＴは、行列の転置をあらわす。また、数式９において、Ｗはスムージングフィルタであり、Ｗｃは局所的フィルタ、Ｗｓは広域的フィルタを示す。具体的には、例えば、Ｗｃは半値幅８ｍｍのガウシアンフィルタを用い、Ｗｓはラプラス行列の逆行列を用いる。
【数１０】
【００９０】
なお、数式１０において、Ｎはチャネル数、Ｔはサンプル数である。また、Ｌは用いる波長の数である。また、Ｔｒは行列のトレースをあらわす。
【数１１】
【数１２】
【００９１】
なお、数式１２において、Ｉは単位行列を示す。
【００９２】
また、数式１１、１２において、Ａはαｉを対角要素に持つ行列（αｉは大脳皮質の各ボクセル毎に在る）である。
【００９３】
つまり、ＡおよびＡ−（−はＡの真上に存在する）は、数式１３のようになる。
【数１３】
【００９４】
数式１３において、αｉはそれぞれのボクセルiにおける、正則化（大きな値にペナルティ）に関するパラメータである。そして、αｉはガンマ分布に従うとする。そして、αｉ−（−はαの真上に存在する）はガンマ分布の平均値である。また、γαｉはガンマ分布のシェイプパラメータである。また、大脳皮質ボクセル数がＭｃ、頭皮ボクセル数がＭｓである。なお、ガンマ分布は、数式１４である。
【数１４】
なお、数式１４においてΓはガンマ関数をあらわす。
【００９５】
また、数式１１、１２において、βは頭皮における正則化に関するパラメータである。βはガンマ分布に従うとする。そしてβ−（−はβの真上に存在する）はガンマ分布の平均値、γβはガンマ分布のシェイプパラメータである。また、ガンマ分布は、数式１５である。
【数１５】
【００９６】
次に、推定手段２８３は、第二ステップを実行する。第二ステップとは、γαｉ、αｉ−（−はαの真上に配置）、γβ、β−（−はβの真上に配置）を、以下の数式１６から数式１９を用いて算出する処理である。なお、γαｉは大脳皮質のボクセルにおけるガンマ分布のシェイプパラメータ、αｉ−は大脳皮質のボクセルにおけるガンマ分布の平均値、γβは頭皮のボクセルにおけるガンマ分布のシェイプパラメータ、β−は頭皮のボクセルにおけるガンマ分布の平均値である。また、上記の初期値は、γα０ｉ、α０ｉ−（−はαの真上に配置）、γβ０、β０−（−はαの真上に配置）であるとし、各初期値は予め与えられている、とする。
【数１６】
【００９７】
数式１６において、Ｔは、サンプル数である。
【数１７】
【数１８】
【数１９】
【００９８】
そして、推定手段２８３は、第一ステップと第二ステップとを繰り返し、第二ステップにおいて、αｉ−（−はαの真上に配置）、γαｉ、β−（−はβの真上に配置）、γβの値の変化が所定の条件を満たすほど小さい場合（ｎ回目の値と（ｎ−１）回目の値との差が閾値以内など）、収束したとして、処理を終了する。なお、最後は、第一ステップを行って終了する。そして、以下の数式２０におけるＸが取得される。
【数２０】
【００９９】
上記において、各行列は以下のような次元数を持つ。Ｙは（Ｎ＊Ｌ）×Ｔ次元行列、Ｇｃは（Ｎ＊Ｌ）×Ｍｃ次元行列、Ｇｓは（Ｎ＊Ｌ）×Ｍｓ次元行列、ＸｃはＭｃ×Ｔ次元行列、ＸｓはＭｓ×Ｔ次元行列である。
（実験１）
【０１００】
次に、本脳活動情報出力システムの第一の実験結果について説明する。本実験において、脳活動情報出力システムを構成する近赤外分光計測装置１は、２０のチャネルを有する。チャネルとは、送光プローブ１１と受光プローブ１２のすべての組み合わせである。つまり、本実験において、図６に示すように、近赤外分光計測装置１は５つの送光プローブ１１と４つの受光プローブ１２を有する。また、図７に示すように、各送光プローブ１１の識別子と各受光プローブ１２の識別子との組み合わせが２０セット存在する。なお、図６において、黒丸が送光プローブ１１、白丸が受光プローブ１２である。また、本実験において、最近傍のプローブの間隔はここでは１５ｍｍとした。また、ボクセルの大きさは４ｍｍ、縦１１ボクセル、横１１ボクセル、深さ５ボクセルとした。さらに、深さ１は頭皮、深さ２・３は頭蓋骨・脳脊髄液、深さ４・５は大脳皮質であるとした。
【０１０１】
そして、以下のようなシミュレーション実験を行った。つまり、シミュレーション実験において、被験者の脳活動があまり表れない安静状態で、近赤外分光計測装置１の送光プローブ１１は、波長７８０ｎｍと８３０ｎｍで発光した、と仮定した。そして、受光プローブ１２が受光する光量の対数値を、各波長において、チャネルごとに、２０回分、取得した。そして、脳活動情報出力装置２の安静時光強度取得部２４は、受光プローブ１２が受光する光量の対数値の平均と分散共分散行列を算出した。なお、図８は、波長７８０ｎｍにおける前記光量の対数値の平均値である。なお、図示しないが、波長８３０ｎｍにおける前記光量の対数値の平均も、同様に、取得した。図９は２波長のデータをまとめた、（チャネル数）×（波長数）＝４０次元を、複数回（今回は２０回）計測したデータから計算した分散共分散行列の一部（波長７８０nmに関する部分）である。そして、安静時光強度蓄積部２５は、安静時光強度情報格納部２３に、波長７８０ｎｍと波長８３０ｎｍにおける光量の対数値の平均および分散共分散行列を蓄積した。
【０１０２】
次に、脳活動情報出力装置２の脳活動時光強度情報取得部２６は、被験者の脳活動がある状態で、波長７８０ｎｍと波長８３０ｎｍにおける光量の対数値を、チャネルごとに取得した。取得したサンプル数は、今回は１回である。この光量の対数値は、図１０の「ｔｒｉａｌ１」の属性値である。なお、被験者の脳活動がある状態では、頭皮の血流変化もある。
【０１０３】
次に、光強度変化情報取得部２７は、波長７８０ｎｍと波長８３０ｎｍの各波長において、チャネルごと光強度変化情報（数式４のＹ）を算出した。つまり、光強度変化情報取得部２７は、数式４に従って、波長ごと、チャネルごと、およびサンプルごとに、光量の対数値の平均値から、脳活動時光強度情報取得部２６が取得した光量の対数値を減算し、光強度変化情報（Ｙ）を算出した。その光強度変化情報（Ｙ）が、図１０の属性値（Ｙ）である。なお、波長780nmにおける観測値Yに関して、Y(1:20)/max(abs(Y))を図示したのが、図１１（ａ）である。また、図１１（ａ）（ｂ）において、線の明るさは、その線の端にある送光プローブと受光プローブのセットからなるチャネルにおける光強度変化情報を表している。図１１（ａ）において、広がった頭皮血流アーチファクトによって、全体的に値が出ていることがわかる。しかし、図１１（ａ）を注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。また、波長830nmにおける観測値Yに関して、Y(21:40)/max(abs(Y))を図示したのが、図１１（ｂ）である。図１１（ｂ）において、広がった頭皮血流アーチファクトによって、全体的に値が出ている。しかし、図１１（ｂ）を注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。なお、780nmよりも強く値が出ているが、これはシミュレーション設定によるものである。
【０１０４】
次に、推定手段２８３は、感度情報格納手段２８１からヘモグロビン濃度変化に関する感度行列（上述したＧｃｏｘｙ，Ｇｃｄｅｏｘｙ，Ｇｓｏｘｙ，Ｇｓｄｅｏｘｙ）を読み出す。次に、推定手段２８３は、大脳皮質活動モデル情報格納部２１から大脳皮質活動モデル情報（例えば、局所的フィルタＷｃ）を読み出す。また、推定手段２８３は、頭皮血流変化モデル情報格納部２２から頭皮血流変化モデル情報（例えば、広域的フィルタＷｓ）を読み出す。
【０１０５】
次に、推定手段２８３は、以下のように変分ベイズアルゴリズムを用いて、Ｘを推定する。つまり、推定手段２８３は、予め格納している初期値を読み出す。初期値は、ここでは、すべてのボクセルに対して「γα０ｉ＝０」、かつ、すべてのボクセルに対して「α０ｉ−＝１００００」、かつ、「γβ０＝０」、かつ、「β０−＝１００００」である。
【０１０６】
そして、推定手段２８３は、上述したステップ１とステップ２の繰り返し処理を１０００回実行し、Ｘを算出した。推定手段２８３が算出したＸの要素であり、頭皮（深さ１）のオキシヘモグロビンの濃度変化（Xoxy）の推定値の集合を図１２に示す。図１２において、各セルの値は、各ボクセルに対応する。したがって、図１２において、１１×１１の値が存在する。また、数値は行（下方向）が頭部モデルの横方向、列（右方向）が頭部モデルの縦方向と対応する。かかることは、図１３、図１４等と同様である。また、推定手段２８３が算出したＸの要素であり、大脳皮質（深さ４）のオキシヘモグロビンの濃度変化（Xoxy）の推定値の集合を図１３に示す。なお、図示しないが、大脳皮質（深さ５）のオキシヘモグロビンの濃度変化（Xoxy）についても、推定手段２８３は、同様に算出した。また、推定手段２８３が算出したＸの要素であり、頭皮（深さ１）のデオキシヘモグロビンの濃度変化（Xdeoxy）の推定値の集合を図１４に示す。さらに、推定手段２８３が算出したＸの要素であり、大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化（Xdeoxy）の推定値の集合を図１５に示す。なお、図示しないが、大脳皮質（深さ５）のデオキシヘモグロビンの濃度変化（Xoxy）についても、推定手段２８３は、同様に算出した。
【０１０７】
次に、出力部２９は、推定部２８が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。出力部２９の出力例を図１６、図１７に示す。図１６（ａ）は、頭皮（深さ１）のオキシヘモグロビンの濃度変化の出力例である。図１６（ａ）等において、縦軸は頭部モデルの縦方向、横軸は頭部モデルの横方向を示し、その数字の単位はmm(ミリメーター)である。また、白黒の濃淡で表されるＸの値の単位はmM(ミリモーラー)である。図１６（ｂ）は、大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の出力例である。また、図１７（ａ）は、頭皮（深さ１）のデオキシヘモグロビンの濃度変化の出力例である。図１７（ｂ）は、大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の出力例である。
【０１０８】
さらに、シミュレーションで仮定した真のＸの値が図１８から図２１である。図１８は、シミュレーションで仮定した真の頭皮（深さ１）のオキシヘモグロビンの濃度変化の値の集合である。また、図１９は、シミュレーションで仮定した真の大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の値の集合である。また、図２０は、シミュレーションで仮定した真の頭皮（深さ１）のデオキシヘモグロビンの濃度変化の値の集合である。図２１は、シミュレーションで仮定した真の大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の値の集合である。
【０１０９】
シミュレーションで仮定した真の値と、推定部２８が取得した推定値とを比較すると、推定部２８は、脳活動をほぼ正確に推定できていることが分かる。
【０１１０】
次に、上記で説明した脳活動情報出力装置２の推定方法の概念について、図２２を用いて説明する。脳活動情報出力装置２における推定は、推定する値（Ｘ）の数より少ない観測値（Ｙ）を用いて、値（Ｘ）を推定しなければならない。これは逆問題の中でも不良設定問題と呼ばれ、正則化と呼ばれるような、なんらかの制限を入れなければ解くことはできない。今回用いる推定法ではベイズ推定を用いる。ベイズ推定ではその正則化と呼ばれるような制限を事前分布において確率的に組み込む。今回の推定法では、その事前分布において生理学的知見を反映させたような、大脳皮質と頭皮において異なる正則化とモデル化を行う。また、今回、ベイズ推定を現実的に行う方法の一例として、変分ベイズ法を用いる。また、値（Ｘ）を推定するために、脳活動情報出力装置２は、大脳皮質活動モデル情報として、スパース化を含む正則化のプログラムおよび局所的フィルタ（Ｗｃ）を保持しており、かつ、頭皮血流変化モデル情報として、スパース化を含まない正則化のプログラムおよび広域的フィルタ（Ｗｓ）を保持している。
【０１１１】
図２２において、Ｙ１からＹＬが観測値である。図２２の（４）における観測値YとXの関係は、光拡散過程においてRytov近似により吸収係数の差δμａが十分小さければ数式２１であらわされる。
【数２１】
数式２１を、数式３、数式５、数式６を用いて書き換えると数式２２となる。しかし、数式２２は観測ノイズを考慮していない関係式である。
【数２２】
【０１１２】
推定部２８は、図２２の（４）において、観測値ＹはＧＸにノイズσ−１Σｙがのって出てくると仮定する。この仮定は、脳活動時のノイズが安静時のノイズの定数倍でうまく近似できるという考えに基づく。
【数２３】
【０１１３】
なお、数２３において、Nは正規分布をあらわす。例えば、N(x|μ,Σｘ)はxが平均μ、分散共分散行列Σｘの正規分布に従うことをあらわす。
【０１１４】
そして、推定部２８は、図２２の（３）において、生理学的知見により、大脳皮質には局所的フィルタ、頭皮は広域的フィルタをかける。また、フィルタをかける前のＸの値をＺとおく。かかる処理の具体的な演算を数式２４に示す。
【数２４】
【０１１５】
Ｚの値に制限をかけないと、逆問題は解くことはできない。脳活動情報出力装置２は、以下のように、確率的に条件をつける。推定部２８は、図２２の（２）において、頭皮においては、頭皮血流変化モデル情報に従って、全部のボクセルで共通な正則化を行う。ここでいう正則化とは、大きなＺの値（つまり大きな濃度変化）についてペナルティをかけることである。これは、図２２の（３）の広域的フィルタ（Ｗｓ）と合わせて、空間的に広がった活動を意味する。かかる処理の具体的な演算を数式２５に示す。
【数２５】
【０１１６】
数式２５において、σは脳活動時の分散共分散行列の安静時に比べての割合の逆数である。また、Ｔはサンプル数である。βは頭皮の正則化に関するパラメータである。Iは単位行列である。
【０１１７】
また、推定部２８は、大脳皮質においては、各ボクセル毎に異なった正則化を行う。ここでは、少数のボクセルのＺの値が大きくなるように選ばれるようになっている（これをスパース化という。）。これは、図２２の（３）の局所的フィルタ（Ｗｃ）と合わせて、局在した活動を意味する。かかる処理の具体的な演算を数式２６に示す。
【数２６】
【０１１８】
数式２６において、Ａはαｉを対角要素に持つ行列であり、パラメータαｉは大脳皮質の各ボクセル毎に存在する正則化に関するパラメータである。
【０１１９】
また、パラメータαｉ、β、σの分布は、数式２７のように仮定する、とする。
【数２７】
【０１２０】
なお、数式２７において、Gammaはガンマ分布を示す。
【０１２１】
以上のように、推定部２８は、確率モデル化した。以下、ベイズ推定により、観測値Yから、最適なσ、Ａ、βの値と、Ｘの推定値を求める。
【数２８】
【０１２２】
数式２８はベイズの定理である。例えば、左辺の事後分布を最大とするようなＸ、σ、Ａ、β、またはＸ、σ、Ａ、βの期待値を推定値とすればよい。しかしながら、直接計算することは難しいため、今回の推定法は変分ベイズ法を用いる。変分ベイズ法では、事後分布を直接計算する代わりに、試験事後分布Qを用意し、真の事後分布に近似させる。この近似は数式２９で定義される自由エネルギーFの最大化によって行われる。
【数２９】
【０１２３】
今回は試験分布におけるＸ、σとＡ、βの間に数３０で示すように独立分解仮定をおく。これにより、計算を簡単に行うことが可能となる。
【数３０】
【０１２４】
この仮定により自由エネルギーの最大化は、Q1に関してFを最大化する「Ｘ、σ最大化ステップ」と、Q2に関してFを最大化する「Ａ、β最大化ステップ」を交互に繰り返すことにより実現できる。「Ｘ、σ最大化ステップ」に対応するアルゴリズムが、前述した第一ステップである。また、「Ａ、β最大化ステップ」に対応するアルゴリズムが、前述した第二ステップである。
【０１２５】
値が収束するまで繰り返した後、今回の推定においては、試験事後分布QにおけるＸ、σ、Ａ、βの期待値を推定値とした。
（実験２）
【０１２６】
次に、本脳活動情報出力システムの第二の実験結果について説明する。本実験において、最近傍のプローブの間隔は２０ｍｍとした。実験２における他の条件は、実験１と同じである。
【０１２７】
そして、以下のようなシミュレーション実験を行った。つまり、シミュレーション実験において、被験者の脳活動があまり表れない安静状態で、近赤外分光計測装置１の送光プローブ１１は、波長７８０ｎｍと８３０ｎｍで発光した、と仮定した。そして、受光プローブ１２が受光する光量の対数値を、各波長において、チャネルごとに、２０回分、取得した。そして、脳活動情報出力装置２の安静時光強度取得部２４は、受光プローブ１２が受光する光量の対数値の平均と分散共分散行列を算出した。なお、図２５は、波長７８０ｎｍにおける前記光量の対数値の平均値である。なお、図示しないが、波長８３０ｎｍにおける前記光量の対数値の平均も、同様に、取得した。図２６は２波長のデータをまとめた、（チャネル数）×（波長数）＝４０次元を、複数回（今回は２０回）計測したデータから計算した分散共分散行列の一部（波長７８０nmに関する部分）である。そして、安静時光強度蓄積部２５は、安静時光強度情報格納部２３に、波長７８０ｎｍと波長８３０ｎｍにおける光量の対数値の平均および分散共分散行列を蓄積した。
【０１２８】
次に、脳活動情報出力装置２の脳活動時光強度情報取得部２６は、被験者の脳活動がある状態で、波長７８０ｎｍと波長８３０ｎｍにおける光量の対数値を、チャネルごとに取得した。取得したサンプル数は、今回は１回である。この光量の対数値は、図２７の「ｔｒｉａｌ１」の属性値である。なお、被験者の脳活動がある状態では、頭皮の血流変化もある。
【０１２９】
次に、光強度変化情報取得部２７は、波長７８０ｎｍと波長８３０ｎｍの各波長において、チャネルごと光強度変化情報（数式４のＹ）を算出した。つまり、光強度変化情報取得部２７は、数式４に従って、波長ごと、チャネルごと、およびサンプルごとに、光量の対数値の平均値から、脳活動時光強度情報取得部２６が取得した光量の対数値を減算し、光強度変化情報（Ｙ）を算出した。その光強度変化情報（Ｙ）が、図２７の属性値（Ｙ）である。なお、波長780nmにおける観測値Yに関して、Y(1:20)/max(abs(Y))を図示したのが、図２８（ａ）である。また、図２８（ａ）（ｂ）において、線の明るさは、その線の端にある送光プローブと受光プローブのセットからなるチャネルにおける光強度変化情報を表している。図２８（ａ）において、広がった頭皮血流アーチファクトによって、全体的に値が出ていることがわかる。しかし、図２８（ａ）を注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。また、波長830nmにおける観測値Yに関して、Y(21:40)/max(abs(Y))を図示したのが、図２８（ｂ）である。図２８（ｂ）において、広がった頭皮血流アーチファクトによって、全体的に値が出ている。しかし、図２８（ｂ）を注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。なお、780nmよりも強く値が出ているが、これはシミュレーション設定によるものである。
【０１３０】
次に、推定手段２８３は、感度情報格納手段２８１からヘモグロビン濃度変化に関する感度行列（上述したＧｃｏｘｙ，Ｇｃｄｅｏｘｙ，Ｇｓｏｘｙ，Ｇｓｄｅｏｘｙ）を読み出す。次に、推定手段２８３は、大脳皮質活動モデル情報格納部２１から大脳皮質活動モデル情報（例えば、局所的フィルタＷｃ）を読み出す。また、推定手段２８３は、頭皮血流変化モデル情報格納部２２から頭皮血流変化モデル情報（例えば、広域的フィルタＷｓ）を読み出す。
【０１３１】
次に、推定手段２８３は、変分ベイズアルゴリズムを用いて、Ｘを推定する。つまり、推定手段２８３は、予め格納している初期値を読み出す。初期値は、ここでは、すべてのボクセルに対して「γα０ｉ＝０」、かつ、すべてのボクセルに対して「α０ｉ−＝１００００」、かつ、「γβ０＝０」、かつ、「β０−＝１００００」である。
【０１３２】
そして、推定手段２８３は、上述したステップ１とステップ２の繰り返し処理を１０００回実行し、Ｘを算出した。推定手段２８３が算出したＸの要素であり、頭皮（深さ１）のオキシヘモグロビンの濃度変化（Xoxy）の推定値の集合を図２９に示す。図２９において、各セルの値は、各ボクセルに対応する。したがって、図２９において、１１×１１の値が存在する。また、数値は行（下方向）が頭部モデルの横方向、列（右方向）が頭部モデルの縦方向と対応する。かかることは、図３０、図３１等と同様である。また、推定手段２８３が算出したＸの要素であり、大脳皮質（深さ４）のオキシヘモグロビンの濃度変化（Xoxy）の推定値の集合を図３０に示す。なお、図示しないが、大脳皮質（深さ５）のオキシヘモグロビンの濃度変化（Xoxy）についても、推定手段２８３は、同様に算出した。また、推定手段２８３が算出したＸの要素であり、頭皮（深さ１）のデオキシヘモグロビンの濃度変化（Xdeoxy）の推定値の集合を図３１に示す。さらに、推定手段２８３が算出したＸの要素であり、大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化（Xdeoxy）の推定値の集合を図３２に示す。なお、図示しないが、大脳皮質（深さ５）のデオキシヘモグロビンの濃度変化（Xoxy）についても、推定手段２８３は、同様に算出した。
【０１３３】
次に、出力部２９は、推定部２８が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。出力部２９の出力例を図３３、図３４に示す。図３３（ａ）は、頭皮（深さ１）のオキシヘモグロビンの濃度変化の出力例である。図３３（ａ）等において、縦軸は頭部モデルの縦方向、横軸は頭部モデルの横方向を示し、その数字の単位はmm(ミリメーター)である。また、白黒の濃淡で表されるＸの値の単位はmM(ミリモーラー)である。図３３（ｂ）は、大脳皮質（深さ４）のオキシヘモグロビンの濃度変化の出力例である。また、図３４（ａ）は、頭皮（深さ１）のデオキシヘモグロビンの濃度変化の出力例である。図３４（ｂ）は、大脳皮質（深さ４）のデオキシヘモグロビンの濃度変化の出力例である。
【０１３４】
以上、本実施の形態によれば、精度高く脳活動を推定できる。また、変分ベイズ法を用いることにより、極めて精度高く脳活動を推定できる。
【０１３５】
なお、本実施の形態の実験によれば、７８０ｎｍと８３０ｎｍの２波長を用いた。しかし、異なる波長を２波長以上用いれば良い。
【０１３６】
また、本実施の形態において、頭皮血流変化と脳活動の分離は、高密度プローブ計測によってはじめて可能となる。短い距離のチャネルは浅い頭皮の活動を拾いやすく、長い距離のチャネルは深い脳活動の活動をも拾うからである。つまり、本実施の形態において、プローブの配置密度が１．５ｃｍ以下であることは極めて好適である。また、プローブの配置密度が１．７ｃｍ以下であることは好適である。なお、プローブの配置密度は２ｃｍ以下、または３ｃｍ未満等でも良い。
【０１３７】
また、本実施の形態によれば、頭部が半無限一層モデルであると仮定し、計測データから脳活動を推定した。しかし、脳活動情報出力システムにおいて、ＭＲＩ構造画像などを用いて、頭部が、頭皮・頭蓋骨・脳脊髄液・灰白質などからなる多層モデルであると仮定し、計測データから脳活動を推定しても良い。
【０１３８】
また、本実施の形態によれば、推定部２８が行う推定法は、主として、変分ベイズ法であった。ただし、推定部２８は、逆問題を解くための他の方法を用いても良い。例えば、推定部２８は、数式３１を格納しており、当該数式３１を用いて、Ｘを推定しても良い。
【数３１】
【０１３９】
つまり、推定部２８は、例えば、数式３１のようなコスト関数Ｃを最小にするようなＸを算出すればよい。なお、数式３１において、１項目が誤差項、２項目以降が正則化項と呼ばれる項であり、２項目以降が大きな値にペナルティを与える項になっている。
【０１４０】
数式３１において、大脳皮質Ｘｃには値（ノルム）の１乗、つまり「Ｌ１ノルム」、頭皮Ｘｓには値（ノルム）の２乗、つまり「Ｌ２ノルム」の正則化を与えている。後者のＬ２ノルムでの正則化は、我々が図２２の（２）で頭皮活動に課した正規分布の条件とほぼ同じ働きをする。また、一般に広がった推定値が得られやすいことが知られている。それに対し、前者のＬ１ノルムでの正則化を与えると、スパース（局在した）な推定値が得られることが知られている。
【０１４１】
なお、数式３１において、フィルタＷはかけていないが、上のλのパラメータや、正則化のノルムを０．９乗にする等の細かい値の修正をすれば、フィルタＷをかけた結果と似たようなことが再現できる。
【０１４２】
以上から、数式３１を用いても、変分ベイズ法を用いた場合と似たようなＸの推定は可能である、と言える。
【０１４３】
さらに、本実施の形態における処理は、ソフトウェアで実現しても良い。そして、このソフトウェアをソフトウェアダウンロード等により配布しても良い。また、このソフトウェアをＣＤ−ＲＯＭなどの記録媒体に記録して流布しても良い。なお、このことは、本明細書における他の実施の形態においても該当する。なお、本実施の形態における脳活動情報出力装置を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、記憶媒体に、局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報と、広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報と、送光プローブと受光プローブの組であるプローブセットを１組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納しており、コンピュータを、前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、前記記憶媒体の安静時光強度情報と、前記脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、前記推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部として機能させるためのプログラム、である。
【０１４４】
また、上記プログラムにおいて、前記推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式を格納し得る演算式格納手段と、前記演算式に、前記感度情報と前記光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する推定手段とを具備するものとして、コンピュータを機能させることは好適である。
【０１４５】
また、上記プログラムにおいて、前記推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式であり、コスト関数である演算式を格納し得る演算式格納手段と、前記演算式に、前記感度情報と前記光強度変化情報を代入し、当該演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定手段とを具備するものとして、コンピュータを機能させることは好適である。
【０１４６】
また、図２３は、本明細書で述べたプログラムを実行して、上述した実施の形態の脳活動情報出力装置等を実現するコンピュータの外観を示す。上述の実施の形態は、コンピュータハードウェア及びその上で実行されるコンピュータプログラムで実現され得る。図２３は、このコンピュータシステム３４０の概観図であり、図２４は、コンピュータシステム３４０のブロック図である。
【０１４７】
図２３において、コンピュータシステム３４０は、ＦＤドライブ、ＣＤ−ＲＯＭドライブを含むコンピュータ３４１と、キーボード３４２と、マウス３４３と、モニタ３４４とを含む。
【０１４８】
図２４において、コンピュータ３４１は、ＦＤドライブ３４１１、ＣＤ−ＲＯＭドライブ３４１２に加えて、ＭＰＵ３４１３と、ＣＤ−ＲＯＭドライブ３４１２及びＦＤドライブ３４１１に接続されたバス３４１４と、ブートアッププログラム等のプログラムを記憶するためのＲＯＭ３４１５とに接続され、アプリケーションプログラムの命令を一時的に記憶するとともに一時記憶空間を提供するためのＲＡＭ３４１６と、アプリケーションプログラム、システムプログラム、及びデータを記憶するためのハードディスク３４１７とを含む。ここでは、図示しないが、コンピュータ３４１は、さらに、ＬＡＮへの接続を提供するネットワークカードを含んでも良い。
【０１４９】
コンピュータシステム３４０に、上述した実施の形態の脳活動情報出力装置等の機能を実行させるプログラムは、ＣＤ−ＲＯＭ３５０１、またはＦＤ３５０２に記憶されて、ＣＤ−ＲＯＭドライブ３４１２またはＦＤドライブ３４１１に挿入され、さらにハードディスク３４１７に転送されても良い。これに代えて、プログラムは、図示しないネットワークを介してコンピュータ３４１に送信され、ハードディスク３４１７に記憶されても良い。プログラムは実行の際にＲＡＭ３４１６にロードされる。プログラムは、ＣＤ−ＲＯＭ３５０１、ＦＤ３５０２またはネットワークから直接、ロードされても良い。
【０１５０】
プログラムは、コンピュータ３４１に、上述した実施の形態の脳活動情報出力装置等の機能を実行させるオペレーティングシステム（ＯＳ）、またはサードパーティープログラム等は、必ずしも含まなくても良い。プログラムは、制御された態様で適切な機能（モジュール）を呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいれば良い。コンピュータシステム３４０がどのように動作するかは周知であり、詳細な説明は省略する。
【０１５１】
なお、上記プログラムにおいて、ハードウェアによって行われる処理、は含まれない。
【０１５２】
また、上記プログラムを実行するコンピュータは、単数であってもよく、複数であってもよい。すなわち、集中処理を行ってもよく、あるいは分散処理を行ってもよい。
【０１５３】
また、上記各実施の形態において、各処理（各機能）は、単一の装置（システム）によって集中処理されることによって実現されてもよく、あるいは、複数の装置によって分散処理されることによって実現されてもよい。
【０１５４】
本発明は、以上の実施の形態に限定されることなく、種々の変更が可能であり、それらも本発明の範囲内に包含されるものであることは言うまでもない。
【産業上の利用可能性】
【０１５５】
以上のように、本発明にかかる脳活動情報出力装置は、精度高く脳活動を推定できる、という効果を有し、脳活動推定装置等として有用である。
【符号の説明】
【０１５６】
１ 近赤外分光計測装置
２ 脳活動情報出力装置
１１ 送光プローブ
１２ 受光プローブ
１３ 受付部
１４ 制御部
２１ 大脳皮質活動モデル情報格納部
２２ 頭皮血流変化モデル情報格納部
２３ 安静時光強度情報格納部
２４ 安静時光強度取得部
２５ 安静時光強度蓄積部
２６ 脳活動時光強度情報取得部
２７ 光強度変化情報取得部
２８ 推定部
２９ 出力部
２８１ 感度情報格納手段
２８２ 演算式格納手段
２８３ 推定手段
【特許請求の範囲】
【請求項１】
局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報を格納している大脳皮質活動モデル情報格納部と、
広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報を格納している頭皮血流変化モデル情報格納部と、
送光プローブと受光プローブの組であるプローブセットを１組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納している安静時光強度情報格納部と、
前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、
前記安静時光強度情報格納部の安静時光強度情報と、前記脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、
前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、
前記推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部とを具備する脳活動情報出力装置。
【請求項２】
前記推定部は、
大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、
大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式を格納し得る演算式格納手段と、
前記演算式に、前記感度情報と前記光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する推定手段とを具備する請求項１記載の脳活動情報出力装置。
【請求項３】
前記推定部は、
大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、
大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式であり、コスト関数である演算式を格納し得る演算式格納手段と、
前記演算式に、前記感度情報と前記光強度変化情報を代入し、当該演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定手段とを具備する請求項１記載の脳活動情報出力装置。
【請求項４】
前記プローブの配置密度が１．５ｃｍ以下である請求項１から請求項３いずれか記載の脳活動情報出力装置。
【請求項５】
前記プローブの配置密度が２．０ｃｍ以下である請求項１から請求項３いずれか記載の脳活動情報出力装置。
【請求項６】
記憶媒体に、
局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報と、
広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報と、
送光プローブと受光プローブの組であるプローブセットを１組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納しており、
脳活動時光強度情報取得部、光強度変化情報取得部、推定部、および出力部により実現される脳活動情報出力方法であって、
前記脳活動時光強度情報取得部が、前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得ステップと、
前記光強度変化情報取得部が、前記記憶媒体の安静時光強度情報と、前記脳活動時光強度情報取得ステップで取得された脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得ステップと、
前記推定部が、前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定ステップと、
前記出力部が、前記推定ステップで取得された大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力ステップとを具備する脳活動情報出力方法。
【請求項７】
記憶媒体に、
局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報と、
広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報と、
送光プローブと受光プローブの組であるプローブセットを１組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納しており、
コンピュータを、
前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、
前記記憶媒体の安静時光強度情報と、前記脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、
前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、
前記推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部として機能させるためのプログラム。
【請求項１】
局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報を格納している大脳皮質活動モデル情報格納部と、
広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報を格納している頭皮血流変化モデル情報格納部と、
送光プローブと受光プローブの組であるプローブセットを１組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納している安静時光強度情報格納部と、
前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、
前記安静時光強度情報格納部の安静時光強度情報と、前記脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、
前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、
前記推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部とを具備する脳活動情報出力装置。
【請求項２】
前記推定部は、
大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、
大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式を格納し得る演算式格納手段と、
前記演算式に、前記感度情報と前記光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する推定手段とを具備する請求項１記載の脳活動情報出力装置。
【請求項３】
前記推定部は、
大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、
大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式であり、コスト関数である演算式を格納し得る演算式格納手段と、
前記演算式に、前記感度情報と前記光強度変化情報を代入し、当該演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定手段とを具備する請求項１記載の脳活動情報出力装置。
【請求項４】
前記プローブの配置密度が１．５ｃｍ以下である請求項１から請求項３いずれか記載の脳活動情報出力装置。
【請求項５】
前記プローブの配置密度が２．０ｃｍ以下である請求項１から請求項３いずれか記載の脳活動情報出力装置。
【請求項６】
記憶媒体に、
局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報と、
広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報と、
送光プローブと受光プローブの組であるプローブセットを１組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納しており、
脳活動時光強度情報取得部、光強度変化情報取得部、推定部、および出力部により実現される脳活動情報出力方法であって、
前記脳活動時光強度情報取得部が、前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得ステップと、
前記光強度変化情報取得部が、前記記憶媒体の安静時光強度情報と、前記脳活動時光強度情報取得ステップで取得された脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得ステップと、
前記推定部が、前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定ステップと、
前記出力部が、前記推定ステップで取得された大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力ステップとを具備する脳活動情報出力方法。
【請求項７】
記憶媒体に、
局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報と、
広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報と、
送光プローブと受光プローブの組であるプローブセットを１組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納しており、
コンピュータを、
前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、
前記記憶媒体の安静時光強度情報と、前記脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、
前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、
前記推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部として機能させるためのプログラム。
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１２】
【図１３】
【図１４】
【図１５】
【図１８】
【図１９】
【図２０】
【図２１】
【図２３】
【図２４】
【図２５】
【図２６】
【図２７】
【図２９】
【図３０】
【図３１】
【図３２】
【図１】
【図１１】
【図１６】
【図１７】
【図２２】
【図２８】
【図３３】
【図３４】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１２】
【図１３】
【図１４】
【図１５】
【図１８】
【図１９】
【図２０】
【図２１】
【図２３】
【図２４】
【図２５】
【図２６】
【図２７】
【図２９】
【図３０】
【図３１】
【図３２】
【図１】
【図１１】
【図１６】
【図１７】
【図２２】
【図２８】
【図３３】
【図３４】
【公開番号】特開２０１２−１８７３７６（Ｐ２０１２−１８７３７６Ａ）
【公開日】平成２４年１０月４日（２０１２．１０．４）
【国際特許分類】
【出願番号】特願２０１１−１１０５５３（Ｐ２０１１−１１０５５３）
【出願日】平成２３年５月１７日（２０１１．５．１７）
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成２２年度、独立行政法人情報通信研究機構「高度通信・放送研究開発委託研究／複数モダリティー統合による脳活動計測技術の研究開発」、産業技術力強化法第１９条の適用を受ける特許出願
【出願人】（３９３０３１５８６）株式会社国際電気通信基礎技術研究所 (905)
【Ｆターム（参考）】
【公開日】平成２４年１０月４日（２０１２．１０．４）
【国際特許分類】
【出願日】平成２３年５月１７日（２０１１．５．１７）
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成２２年度、独立行政法人情報通信研究機構「高度通信・放送研究開発委託研究／複数モダリティー統合による脳活動計測技術の研究開発」、産業技術力強化法第１９条の適用を受ける特許出願
【出願人】（３９３０３１５８６）株式会社国際電気通信基礎技術研究所 (905)
【Ｆターム（参考）】
[ Back to top ]