(57)【要約】【課題】局所神経活動に伴う信号変化とその他の起源を有する生理学的変動が混在する場合においても、信号を生理学的起源の異なる構成成分に精度よく分別し、構成成分別の情報を提供する。【解決手段】生体機能データを計測し（１）、必要に応じて解析領域の選択（２）、前処理（３）および入力信号の生成（４）を行い、生体機能計測データに対して画素またはチャンネル毎に入出力関係を記述する数式モデルを構築する（５）。これにより生体機能計測データの信号に含まれる生理学的起源の異なる構成成分の生体機能情報を分別抽出する。抽出結果は、例えば表示装置に表示される。【選択図】　　　図１
2004041670.tif 000002
【特許請求の範囲】
【請求項１】
生体機能計測データに対して画素またはチャンネル毎に入出力関係を記述する数式モデルを構築し、生体機能計測データの信号に含まれる生理学的起源の異なる構成成分の生体機能情報を分別抽出することを特徴とする生体機能情報の分別抽出方法。
【請求項２】
抽出される生体機能が中枢神経系の機能であることを特徴とする請求項１記載の生体機能情報の分別抽出方法。
【請求項３】
抽出される生体機能が局所生理機能であることを特徴とする請求項１記載の生体機能情報の分別抽出方法。
【請求項４】
数式モデルの構築において、生体機能計測データの代表信号値を使用することを特徴とする請求項１〜３のいずれかに記載の生体機能情報の分別抽出方法。
【請求項５】
数式モデルの構築において、全身生理機能の測定値を使用することを特徴とする請求項１〜３のいずれかに記載の生体機能情報の分別抽出方法。
【請求項６】
数式モデルの構築において、確率的雑音を使用することを特徴とする請求項１〜３のいずれかに記載の生体機能情報の分別抽出方法。
【請求項７】
分別抽出された生体機能情報を表示することを特徴とする請求項１〜６のいずれかに記載の生体機能情報の分別抽出方法。
【請求項８】
生体機能情報の表示において、モデルに基づくシミュレーションを使用することを特徴とする請求項７記載の生体機能情報の分別抽出方法。
【請求項９】
生体機能情報の表示において、ノイズ寄与率（パワー寄与率）の分布または空間的情報の表示を使用することを特徴とする請求項７記載の生体機能情報の分別抽出方法。
【請求項１０】
生体機能情報の表示において、モデル固有値の分布または空間的情報の表示を使用することを特徴とする請求項７記載の生体機能情報の分別抽出方法。
【請求項１１】
生体機能情報の表示において、確率的雑音の強度の分布または空間的情報の表示を使用することを特徴とする請求項７記載の生体機能情報の分別抽出方法。
【請求項１２】
生体機能計測データに対して画素またはチャンネル毎に入出力関係を記述する数式モデルを構築し、生体機能計測データの信号に含まれる生理学的起源の異なる構成成分の生体機能情報を分別抽出する手段を備えたことを特徴とする生体機能情報の分別抽出装置。
【請求項１３】
分別抽出された生体機能情報を表示する手段を備えたことを特徴とする請求項１２記載の生体機能情報の分別抽出装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、磁気共鳴映像法や近赤外線分光法などの生体機能計測データから、生理学的起源の異なる生体機能に由来する構成成分を分別抽出し、各成分の情報を表示する生体機能情報の分別抽出方法およびその装置に関する。
【０００２】
【従来の技術】
磁気共鳴映像法、近赤外線分光法などの生体機能測定法は非侵襲的な脳機能計測法として近年、広く利用されている。これらの計測データには局所神経活動に伴う信号変化のほか、全身循環、呼吸、自律神経活動、血管運動などさまざまな生理学的現象に由来する信号変化の影響が混在している。生理学的現象のうち計測が容易な心拍、呼吸の影響については、これらの周波数成分と生体機能計測データに含まれる周波数成分の比較が行われ、生体機能計測データにこれらの生理学的現象の成分が含まれていることが明らかにされている（非特許文献１、２、３）。
【非特許文献１】
Ｏｇａｗａ　Ｓ，　ｅｔ　ａｌ　Ｖｉｓｕａｌｉｚａｔｉｏｎ　ｏｆ　Ｉｎｆｏｒｍａｔｉｏｎ　Ｐｒｏｃｅｓｓｉｎｇ　ｉｎ　ｔｈｅ　Ｈｕｍａｎ　Ｂｒａｉｎ：　Ｒｅｃｅｎｔ　Ａｄｖａｎｃｅｓ　ｉｎ　ＭＥＧ　ａｎｄ　Ｆｕｎｃｔｉｏｎａｌ　ＭＲＩ　（ＥＥＧ　Ｓｕｐｐｌ．　４７）　１９９６，　ｐｐ５−１４．
【非特許文献２】
Ｍｉｔｒａ　ＰＰ，　ｅｔ　ａｌ，　Ｍａｇｎ　Ｒｅｓｏｎ　Ｍｅｄ　１９９７；３７：５１１−８．
【非特許文献３】
Ｌｏｗｅ　ＭＪ，　ｅｔ　ａｌ，　Ｎｅｕｒｏｉｍａｇｅ　１９９８；７：１１９−３２．
脳機能検査法としての応用では局所神経活動とこれら生理学的現象に伴う信号変化の周波数成分の違いを利用して、トレンド除去、ローパスフィルターの適用などの前処理が行われてきた（非特許文献４、５）。
【非特許文献４】
Ｆｒｉｓｔｏｎ　ＫＪ，　ｅｔ　ａｌ，　Ｈｕｍａｎ　Ｂｒａｉｎ　Ｍａｐｐｉｎｇ　１９９５；２；１８９−２１０．
【非特許文献５】
Ｓｔｒｕｐｐ　ＪＰ，　Ｎｅｕｒｏｉｍａｇｅ　１９９６；３；Ｓ６０７．
従来は局所神経活動に伴う変化を抽出するためには神経活動賦活時の計測値と基準となる計測値の間の統計的検定を行い有意差のある部分を抽出する方法や、典型的な信号変化パターンとの相関を求める方法などが用いられてきた（非特許文献５）。さらにこれらの手法を拡張した一般線形モデルも用いられている（非特許文献４）。また、関心領域の反応パターンを推定するためには賦活パラダイムに同期した加算平均や事前に予測された変化パターンに基づくフィッティングなどが用いられてきた。最近では独立成分分析を用いて由来の異なる信号成分の分別が試みられている（非特許文献６）。
【非特許文献６】
Ｇｕ　Ｈ，　ｅｔ　ａｌ，　Ｎｅｕｒｏｉｍａｇｅ　２００１；１４：１４３２−４３．
一方、その他の生理学的現象に由来する信号は外乱として認識されており、これらの生体機能計測データからの局所生理機能抽出の応用は行われていなかった。
【０００３】
【発明が解決しようとする課題】
生体機能計測データに影響を与える生理学的現象としては循環、呼吸、自律神経活動、血管運動などがあげられる。前記で述べた方法では、局所神経活動に伴う信号変化とは大きく周波数が異なる生理学的変動、たとえば心拍、呼吸の影響や数分単位の血圧、動脈血二酸化炭素分圧の変動などの影響は前処理により容易に除去可能であるが、周波数帯域の重なる変化、たとえば圧受容体反射に伴う血圧変動の影響や血管運動の影響は除去できない。したがって、これらの影響を除去するにはデータサンプル数を増やして統計パワーを増加させる必要がある。しかしながらこれらの生理学的変動は賦活パラダイムとは独立した現象とは限らず、データサンプル数の増加のみではバイアスのない検定、推定結果がえられるとは限らない。また、これらの影響が大きい場合には誤った結果を導き出す可能性がある。独立成分分析を用いた場合、分離された信号成分の生理学的な解釈は必ずしも容易ではない。そして、これら生理学的現象が生体機能計測データに及ぼす影響の強度は空間的に一様ではない。また、時間的にも同期的、非同期的な影響が混在している。
【０００４】
本発明は、局所神経活動に伴う信号変化とその他の起源を有する生理学的変動が混在する場合においても、信号を生理学的起源の異なる構成成分に精度よく分別し、構成成分別の情報を提供することを目的とする。また本発明は、生体機能計測データに含まれるこれら生理学的起源を異にする信号成分の時間的あるいは空間的特徴を容易に理解可能な形として提供することを目的とする。
【０００５】
【課題を解決するための手段】
上記目的は、生体機能計測データに対して画素またはチャンネル毎に入出力関係を記述する数式モデルを構築し、生体機能計測データの信号に含まれる生理学的起源の異なる構成成分の生体機能情報を分別抽出する生体機能情報の分別抽出方法により、達成される。
ここで、抽出される生体機能は中枢神経系の機能または局所生理機能とすることができる。また、数式モデルの構築において、生体機能計測データの代表信号値を使用すること、全身生理機能の測定値を使用すること、または確率的雑音を使用することができる。
また、分別抽出された生体機能情報を表示することができる。生体機能情報の表示において、モデルに基づくシミュレーションを使用すること、ノイズ寄与率（パワー寄与率）の分布または空間的情報の表示を使用すること、モデル固有値の分布または空間的情報の表示を使用すること、または確率的雑音の強度の分布または空間的情報の表示を使用することができる。
本発明に係る生体機能情報の分別抽出装置は、生体機能計測データに対して画素またはチャンネル毎に入出力関係を記述する数式モデルを構築し、生体機能計測データの信号に含まれる生理学的起源の異なる構成成分の生体機能情報を分別抽出する手段を備えて構成される。本分別抽出装置は、さらに、分別抽出された生体機能情報を表示する手段を備えることができる。
【０００６】
本発明では、解析に用いる生体機能計測データとしては、磁気共鳴信号、近赤外線スペクトルなどを使用する。これら生体機能計測データにはさまざまな起源を有する生理学的変動が混在している。一般に局所神経活動に伴う変化は賦活パラダイムに対応した応答を示すが、その応答は空間的、時間的に一様なものではない。その他の起源を有する生理学的変動の由来としては心拍、呼吸といった生体の有する基本的な周期的活動、自律神経活動に伴う血圧変動、動脈血二酸化炭素分圧の変動などがあり、これらについては前記生体機能計測データに対して同期的な信号変化をもたらすが、これらによる影響の強度は空間的に一様なものではない。これに対して血管運動、すなわち小血管の周期的収縮・拡張にともなう信号変化は通常、局所ごとに異なる位相を示し、非同期的である。またその強度も空間的に一様ではない。各画素またはチャンネルに対応して入力−出力関係を記述する数式モデルを構築すれば、各種生理学的変動に対する空間的、時間的応答を推定することができる。本発明では、生体機能計測によりえられた画素またはチャンネル毎の計測値を出力とし、機能賦活パラダイム、同時に計測された生理機能の計測値、生体機能計測データから作成された参照信号、確率的雑音などを入力とする数式モデルを構築することで、各入力が生体機能計測値に及ぼす影響を分別抽出する。えられた数式モデルのモデルパラメーターの表示に加え、各入力を独立に変化させたときの出力のシミュレーションを表示することで入力―出力間の時間的、空間的関係を視覚的に容易に理解可能な形で提供し、前記目的を達成する。入力成分間に相互作用がある系においても各成分が他の成分の信号変化に寄与する割合を表示することで各成分間のフィードバック関係を理解容易な形で提供する。また、明らかな入力を特定しえない生理学的変動、たとえば血管運動に起因する信号変化に対しては確率的雑音を入力としてモデルを構築し、そのモデルパラメーターの表示に加え、その生理学的変動成分の特徴を理解容易な形にして提供する。また、本発明では連続した測定データをモデル推定用とモデル評価用に分割して使用することで、構築されたモデルの妥当性を検証することができる。あるいは別個に収集されたデータを結合してモデルを構築することもできるので、データ収集に関する制約が少なくなる。
【０００７】
【発明の実施の形態】
以下、本発明の実施の形態を例とともに説明する。図１は本発明の実施手順を示すフローチャートである。本フローチャートの各ステップ１〜６は次のような処理が行われる。
【０００８】
（１）生体機能データの計測
図２（ａ）、（ｂ）に生体機能および生理学的データの収集方法を示す。磁気共鳴画像（図中ａ１）では高い時間分解能がえられるエコープラナー法を使用することが望ましい。近赤外線スペクトル（図中ａ２）では空間的情報をえるため、多チャンネル計測を使用することが望ましい。脳機能計測では被験者に対し刺激ないし課題遂行を行う機能賦活期間と刺激ないし課題を与えない対照期間のデータを収集する。図中、右方向は時間軸であり、賦活期間と対照期間が交互に付与される。本発明ではデータ収集に際して従来の方法のように一定の機能賦活期間と対照期間の組み合わせを反復して行う必要はなく、上記期間の長さは一回のデータ収集期間の間であっても、図２に示すように、試行毎に変えて差し支えない。局所生理機能の情報抽出のみを目的とする場合には、刺激ないし課題を与える必要はない。生体機能計測データの収集と平行して、生理学的変動の参照データとして脈波、連続血圧（図中ｂ１）、胸腹壁運動、呼気中二酸化炭素濃度（図中ｂ２）などを計測し、生理学的データの収集を行う。なお、あとで示すように生体機能計測データ計測値の代表値を参照データとして使用する場合は、これらの生理学的データの収集は必要とはしない。
【０００９】
（２）解析領域の選択
磁気共鳴画像を生体機能計測データとして使用する場合は解析を行う領域を脳が存在する部分に制限する。その際には、脳が存在する部分と背景の領域の磁気共鳴信号の強度差を利用し、当該信号強度に閾値を設定することで両者を判別することができる。なお、データによってはこの処理を省略しても差し支えない。
近赤外線スペクトルを生体機能計測データとして使用する場合にはこの処理は必要としない。
【００１０】
（３）前処理
（１）、（２）の方法で生成したデータに対し、直流および超低周波数成分を除去する平均値除去、トレンド除去処理を画素またはチャンネル毎に行う。また、高周波数成分を除去するローパスフィルター処理を行う。さらに関心のある周波数帯域に合わせてダウンサンプリングを行う。なお、データによってはこれらの処理を省略しても差し支えない。
【００１１】
（４）入力信号の生成
入力信号として使用する脳機能賦活パラダイムデータと生理学的変動の参照データを生成する。脳機能賦活パラダイムデータは機能賦活期間と対照期間それぞれに任意の値を割り付けて生成するが、期間全体の平均値が０となるように値を設定することが望ましい。生理学的変動の参照データとしては連続血圧、呼気中二酸化炭素濃度などの実測値を使用するが、データ収集時の被験者の負担が少なく、後に述べる数式モデルの構造を簡略化するための方法として（１）、（２）、（３）の処理をへた生体機能計測データの代表値を使用することも可能である。本実施例では生体機能計測データの代表値として全脳平均値を使用する。なお、病態により正常と異なる信号変化を示す部分が生体機能計測データに含まれている場合には、当該の画素またはチャンネルの計測値を除外して処理を行う。多断面の磁気共鳴画像を用いる場合には断面ごとの撮像タイミングのずれを補正するため、断面毎の撮像タイミングのずれを調整した補間処理を行った後、平均値を算出する。これら測定値に対し、（３）と同様にして直流および超低周波数成分を除去する平均値除去、トレンド除去処理を行う。また、高周波数成分を除去するローパスフィルター処理を行う。なお、データによってはこの処理を省略しても差し支えない。生理学的変動の参照データに遅延がある場合、たとえばサイドストリームとして採取した呼気ガスサンプルを用いて測定した二酸化炭素濃度を使用する場合などには、あらかじめ遅延時間を補正した計測値を入力信号として使用する。生理学的変動の参照データのサンプリングレートが生体機能計測データのそれと異なる場合には両者が同じサンプリングレートとなるようにリサンプリング処理により調節する。
【００１２】
（５）数式モデルの生成
利用する数式モデルとしては線形、時不変性モデルに限らず、非線形、時変性モデルを含むあらゆるブラックボックスモデルを使用することができる。生体の活動は一般的には非線形、時変性であるが、安定した生理学的状態、比較的短い時間の範囲では線形、時不変性システムと仮定して取り扱っても大きな誤差は生じない。本実施例ではモデルパラメーターの計算が容易な自己回帰システムモデル（ａｕｔｏ−ｒｅｇｒｅｓｓｉｖｅ　ｅｘｏｇｅｎｏｕｓ　ｍｏｄｅｌ；　ＡＲＸ　ｍｏｄｅｌ）を使用して説明する（中溝高好著，信号解析とシステム同定，コロナ社，足立修一著．ＭＡＴＬＡＢによる制御のためのシステム同定，東京電機大学出版局．Ｌｊｕｎｇ　Ｌ，　Ｓｙｓｔｅｍ　Ｉｄｅｎｔｉｆｉｃａｔｉｏｎ　Ｔｏｏｌｂｏｘ　Ｕｓｅｒ’ｓ　Ｇｕｉｄｅ　Ｖｅｒｓｉｏｎ　５，　Ｔｈｅ　Ｍａｔｈｗｏｒｋｓ，　Ｉｎｃ）。数式を次に示す。
【００１３】
【数１】
2004041670.tif 000003
【００１４】
ただし、ｙ（ｋ）はシステムの出力信号、ｕ（ｋ）はシステムの入力信号、ｎｋはむだ時間、ｗ（ｋ）は平均値０、有限な分散を持つ白色雑音、Ａ（ｑ）とＢ（ｑ）はシフトオペレーターｑの多項式である。
画素またはチャンネル毎の計測値を出力、脳機能賦活パラダイムデータと生理学的変動の参照データを入力としてモデルパラメーターの推定を行う。ＡＲＸモデルパラメーターの推定には最小二乗推定法を使用するが、これ以外のモデル、たとえば自己回帰移動平均システムモデル（ａｕｔｏ−ｒｅｇｒｅｓｓｉｖｅ　ｍｏｖｉｎｇ　ａｖｅｒａｇｅ　ｅｘｏｇｅｎｏｕｓ　ｍｏｄｅｌ；　ＡＲＭＡＸ　ｍｏｄｅｌ）を使用する場合にはモデルパラメーター推定法として最尤推定法などを使用する（中溝高好著，信号解析とシステム同定，コロナ社）。
【００１５】
示適モデル次数はクロスバリデーション、最小実現、情報規範などを用いて選定する（中溝高好著，信号解析とシステム同定，コロナ社，足立修一著，ＭＡＴＬＡＢによる制御のためのシステム同定，東京電機大学出版局）。本実施例では多数の画素またはチャンネルの計測値を処理するにあたって、オペレーターの関与を少なくする方法として情報規範を用いた自動的選定、あるいは自動的に選定されたモデル次数の結果を参考にしたオペレーター選定によった。自動選定に用いる情報規範としては赤池情報量基準などを使用する（赤池弘次、中川東一郎共著，ダイナミックシステムの統計的解析と制御［新訂版］，サイエンス社．尾崎統、北川源四郎編、時系列解析の方法、朝倉書店）。同様にしてむだ時間ｎｋについても相関法を用いたインパルス応答あるいはクロスバリデーション、情報規範などで選定を行う。なお、データによってはむだ時間を０に固定しても差し支えない。
【００１６】
全脳平均値を入力信号として使用した場合には、全脳平均値の変化に先行して信号変化を生ずる脳部位が存在する可能性があるため、全脳平均値の遅延時間の調整を行う。調整量はクロスバリデーションに基づく選定、あるいは既知の通過時間に基づきオペレーターの指示による選定を行う。なお、データによってはこの処理を省略しても差し支えない。
【００１７】
局所生理機能の情報抽出のみを目的とする場合には、生理学的変動の参照データを入力とするモデルを使用する。あるいは複数の生理学的変動の参照データを用い、これらの間のフィードバック関係を含めて情報抽出を行う場合は多変量時系列モデルを使用する（赤池弘次、中川東一郎共著，ダイナミックシステムの統計的解析と制御［新訂版］，サイエンス社．尾崎統、北川源四郎編、時系列解析の方法、朝倉書店）。明らかな入力を特定しえない生理学的変動、たとえば血管運動に起因する信号変化の特徴を抽出する場合は他の生理学的変動の参照データとともに確率的雑音を入力としてモデル化する。あるいは他の生理学的変動の影響が除外できる場合は時系列モデルを使用する。
【００１８】
（６）結果の表示
画素、チャンネル毎に推定されたモデルパラメーターは多項式表現、伝達関数表現、零点−極点−ゲイン表現などに形式を変換して表示することができる。あるいは周波数伝達関数として表示することができる。零点、極点などのモデル固有値は複素平面上での分布表示あるいは周波数帯域別の空間的分布を示すマップ表示として示すことができる。
図３はヒトの視覚刺激時の頭部磁気共鳴画像データに本発明を適用し得られたモデルに基づき、脳機能賦活パラダイムに対応するゲインの空間的分布を示すマップである。ゲインの値は画像の輝度により表示している。
図４はヒトの安静時の頭部磁気共鳴画像データと生理学的変動の参照値（呼気終末二酸化炭素濃度）に本発明を適用し得られたモデルの固有値（極点）の分布を複素平面上に表示したグラフである。
図５および図６は本発明を適用し得られたモデルの固有値（極点および零点）の周波数帯域別の空間的分布を示すマップである。これらの例においては周波数帯域　０．００８〜０．０５　Ｈｚ、０．０５〜０．１０　Ｈｚ、０．１０〜０．１５　Ｈｚ、０．１５〜０．２０　Ｈｚ、０．２０〜０．２５　Ｈｚにおける画素毎の固有値の数を左上から右下に向かって配置された画像によって示す。なお最下段右はブランクである。固有値の数は画像の輝度により表示している。図６では、最上段右の周波数帯域　０．０５〜０．１０　Ｈｚの画像のみが現われている。
確率的雑音の強度は分散あるいは標準偏差、変動係数などに変換して、その分布あるいは空間的分布を示すマップとして表示することができる。
図７はヒトの安静時の頭部磁気共鳴画像データと生理学的変動の参照値（呼気終末二酸化炭素濃度）に本発明を適用し得られたモデルの確率的雑音の強度の空間的分布を示すマップである。確率的雑音の強度は画像の輝度により表示している。この例では表示スケールは画素毎の確率的雑音の分散に１を加えた値の常用対数により示している。
モデルに基づいてシミュレーションを行う場合にはインパルス応答、ステップ応答、その他、オペレーターの指定する任意の入力波形に対する応答を表示する。応答は入力開始時を基準とする経時的変化、あるいはオペレーターの指定する経過時間後の信号値の空間的分布を示すマップとして表示することができる。
【００１９】
生体機能計測の実測値と入力信号に基づいたシミュレーション結果を比較することで、推定されたモデルの妥当性を検証することができる。
図８はヒトの視覚刺激時の頭部磁気共鳴画像データに本発明を適用し得られたモデルに対し、パルス状の入力波形を脳機能賦活パラダイムとして与えたときの脳内の各部位の経時的信号変化のシミュレーション結果を示す図である。この例においては８秒間の強度１のパルス状入力に対する脳内各部位の０．５秒間隔の反応を左上から右下に向かって配置された画像により示す。２段目左端の画像が入力開始時、上から４段目左から４番目の画像が入力終了時に対応する。入力に対応する信号強度の変化量は画像の輝度により表示している。
図９は後頭葉の１画素の信号変化のシミュレーション結果を示すグラフである。横軸は時間経過を入力開始時点を１とする画像番号により表示している。縦軸は信号強度の変化量を示す。参照のため入力波形を重ねて表示している。
図１０は上記モデルに対し、パルス状の入力波形を全脳平均値の変動として与えたときの脳内の各部位の経時的信号変化のシミュレーション結果を示す図である。この例においては８秒間の強度１のパルス状入力に対する脳内各部位の０．５秒間隔の反応を左上から右下に向かって配置された画像により示す。２段目左端の画像が入力開始時、上から４段目左から４番目の画像が入力終了時に対応する。入力に対応する信号強度の変化量は画像の輝度により表示している。
図１１は頭頂葉の１画素の信号変化のシミュレーション結果を示すグラフである。横軸は時間経過を入力開始時点を１とする画像番号により表示している。縦軸は信号強度の変化量を示す。参照のため入力波形を重ねて表示している。
【００２０】
本実施例は、脳機能賦活パラダイムとして１種類の入力を使用する例を示したが、複数の入力を用いる脳機能賦活パラダイムに対しても適用することができる。
本実施例は、１断面の磁気共鳴画像を用いているが、多断面のデータでも各断面に対して同様の手順で取り扱うことができる。多断面の場合には断面毎の撮像タイミングのずれを補正して処理を行う。
本実施例では、全データ収集後の一括処理を例にとって説明したが、逐次計算アルゴリズムを使用することにより（中溝高好著，信号解析とシステム同定，コロナ社，足立修一著，ＭＡＴＬＡＢによる制御のためのシステム同定，東京電機大学出版局）、オンラインあるいはリアルタイム処理が可能である。
本実施例ではヒトを対象とした脳機能賦活情報の抽出を例にとって説明したが、本発明は脳以外の臓器、組織およびヒト以外の動物を対象として適用することができる。
【００２１】
時系列モデルを使用した場合にはモデルパラメーターの表示または形式を変換した表示、固有値の表示、パワースペクトル、クロススペクトル、コヒーレンスなどの周波数領域の情報に変換した表示、あるいは確率的雑音の強度の表示をすることができる。また、モデルに基づくシミュレーションを表示することができる。多変量自己回帰モデルにおいて、各変量の雑音がお互いに無相関とみなせる場合には、変量間にフィードバックが存在する系であっても、閉ループの周波数応答関数、ノイズ寄与率（またはパワー寄与率）を表示することができる（赤池弘次、中川東一郎共著，ダイナミックシステムの統計的解析と制御［新訂版］，サイエンス社．尾崎統、北川源四郎編、時系列解析の方法、朝倉書店）。
図１２はヒトの安静時の頭部磁気共鳴画像データと生理学的変動の参照値（脈波、胸郭運動、呼気終末二酸化炭素濃度）に本発明を適用し得られたモデルに基づき、横軸に周波数、縦軸に磁気共鳴信号に対する各変量のノイズ（パワー）寄与率を表示したグラフである。図中、Ａは呼気終末二酸化炭素濃度、Ｂは胸郭運動、Ｃは脈波、Ｄは自身（磁気共鳴信号）をそれぞれ示す。
また、本発明を適用し得られたモデルに基づき、周波数帯域を選択して生体機能計測データに対する各変量のノイズ（パワー）寄与率の空間的分布をマップとして表示することができる。
図１３は本発明を適用し得られたモデルに基づき、周波数帯域０．９８〜１．１０　Ｈｚにおける磁気共鳴信号に対する脈波のノイズ（パワー）寄与率の空間的分布を示すマップである。脈波のノイズ（パワー）寄与率は画像の輝度により表示している。
図１４は本発明を適用し得られたモデルに基づき、周波数帯域０．２４〜０．３６　Ｈｚにおける磁気共鳴信号に対する胸郭運動のノイズ（パワー）寄与率の空間的分布を示すマップである。胸郭運動のノイズ（パワー）寄与率は画像の輝度により表示している。
図１５は本発明を適用し得られたモデルに基づき、周波数帯域０．０４〜０．１６　Ｈｚにおける磁気共鳴信号に対する呼気終末二酸化炭素濃度のノイズ（パワー）寄与率の空間的分布を示すマップである。呼気終末二酸化炭素濃度のノイズ（パワー）寄与率は画像の輝度により表示している。
本実施例では変量毎にノイズ（パワー）寄与率の空間的分布を示すマップを作成したが、各変量に色を割り付けることにより１枚のカラーマップとして表示することも可能である。
また、任意の変量に対し入力を与えた場合の応答をシミュレーションとして表示することができる（尾崎統、北川源四郎編、時系列解析の方法、朝倉書店．尾崎統、北川源四郎編、時系列解析の実際Ｉ、朝倉書店）。
【００２２】
本発明に係る生体機能情報の分別抽出装置としては、例えば、コンピュータを用いて生体機能計測データに対して画素またはチャンネル毎に入出力関係を記述する数式モデルを構築し、生体機能計測データの信号に含まれる生理学的起源の異なる構成成分の生体機能情報を分別抽出するようにすることができる。分別抽出された生体機能情報はコンピュータに接続された表示装置に表示することができる。
【発明の効果】
本発明では、既知の抽出法である平均値の有意差検定や相関を用いた方法とは異なり、刺激提示あるいは課題遂行に伴う生体機能計測データの信号変化のパターンを仮定する必要がない。また、本発明では統計的検定の前提となる事象の独立性を満足するための平滑化処理を必要としないので生体機能計測データの分解能を劣化させることがない。本発明では独立成分分析とは異なり、分別された信号成分の生理学的由来について推測する必要がない。本発明では、局所神経活動以外の生理学的変動が生体機能計測データに混在する場合でも精度よく局所神経活動の情報を抽出することができる。さらに、局所神経活動以外の種々の生理機能についての情報を抽出することができる。
【図面の簡単な説明】
【図１】生理機能賦活情報の分別抽出手順を示すフローチャートである。
【図２】（ａ）、（ｂ）は生体機能および生理学的データの収集方法を示す図である。
【図３】本発明を適用し得られたモデルに基づき、脳機能賦活パラダイムに対応するゲインの空間的分布を示すマップである。
【図４】本発明を適用し得られたモデルの固有値（極点）の分布を複素平面上に表示したグラフである。
【図５】本発明を適用し得られたモデルの固有値（極点）の周波数帯域別の空間的分布を示すマップである。
【図６】本発明を適用し得られたモデルの固有値（零点）の周波数帯域別の空間的分布を示すマップである。
【図７】本発明を適用し得られたモデルの確率的雑音の強度の空間的分布を示すマップである。
【図８】本発明を適用し得られたモデルに対し、パルス状の入力波形を脳機能賦活パラダイムとして与えたときの脳内の各部位の経時的信号変化のシミュレーション結果を示す図である。
【図９】本発明を適用し得られたモデルに対し、パルス状の入力波形を脳機能賦活パラダイムとして与えたときの１画素の信号変化のシミュレーション結果を示すグラフである。
【図１０】本発明を適用し得られたモデルに対し、パルス状の入力波形を全脳平均値の変動として与えたときの脳内の各部位の経時的信号変化のシミュレーション結果を示す図である。
【図１１】本発明を適用し得られたモデルに対し、パルス状の入力波形を全脳平均値の変動として与えたときの１画素の信号変化のシミュレーション結果を示すグラフである。
【図１２】本発明を適用し得られたモデルに基づき、磁気共鳴信号に対する各変量のノイズ（パワー）寄与率を表示したグラフである。
【図１３】本発明を適用し得られたモデルに基づき、周波数帯域０．９８〜１．１０　Ｈｚにおける磁気共鳴信号に対する脈波のノイズ（パワー）寄与率の空間的分布を示すマップである。
【図１４】本発明を適用し得られたモデルに基づき、周波数帯域０．２４〜０．３６　Ｈｚにおける磁気共鳴信号に対する胸郭運動のノイズ（パワー）寄与率の空間的分布を示すマップである。
【図１５】本発明を適用し得られたモデルに基づき、周波数帯域０．０４〜０．１６　Ｈｚにおける磁気共鳴信号に対する呼気終末二酸化炭素濃度のノイズ（パワー）寄与率の空間的分布を示すマップである。
【符号の説明】
１　生体機能データの計測
２　解析領域の選択
３　前処理
４　入力信号の生成
５　数式モデルの生成
６　結果の表示
【図１】
2004041670.tif 000004
【図２】
2004041670.tif 000005
【図３】
2004041670.tif 000006
【図４】
2004041670.tif 000007
【図５】
2004041670.tif 000008
【図６】
2004041670.tif 000009
【図７】
2004041670.tif 000010
【図８】
2004041670.tif 000011
【図９】
2004041670.tif 000012
【図１０】
2004041670.tif 000013
【図１１】
2004041670.tif 000014
【図１２】
2004041670.tif 000015
【図１３】
2004041670.tif 000016
【図１４】
2004041670.tif 000017
【図１５】
2004041670.tif 000018