TOP > 国内特許検索 > 脳機能データ解析方法、脳機能解析装置及び脳機能解析プログラム > 明細書

明細書 :脳機能データ解析方法、脳機能解析装置及び脳機能解析プログラム

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第4308871号 (P4308871)
登録日 平成21年5月15日(2009.5.15)
発行日 平成21年8月5日(2009.8.5)
発明の名称または考案の名称 脳機能データ解析方法、脳機能解析装置及び脳機能解析プログラム
国際特許分類 A61B   5/055       (2006.01)
FI A61B 5/05 380
A61B 5/05 382
請求項の数または発明の数 27
全頁数 52
出願番号 特願2007-539917 (P2007-539917)
出願日 平成18年10月6日(2006.10.6)
国際出願番号 PCT/JP2006/320078
国際公開番号 WO2007/043462
国際公開日 平成19年4月19日(2007.4.19)
優先権出願番号 2005298236
優先日 平成17年10月12日(2005.10.12)
優先権主張国 日本国(JP)
審査請求日 平成20年8月20日(2008.8.20)
特許権者または実用新案権者 【識別番号】800000068
【氏名又は名称】学校法人東京電機大学
発明者または考案者 【氏名】月本 洋
早期審査対象出願または早期審理対象出願 早期審査対象出願
個別代理人の代理人 【識別番号】100079108、【弁理士】、【氏名又は名称】稲葉 良幸
【識別番号】100109346、【弁理士】、【氏名又は名称】大貫 敏史
【識別番号】100117189、【弁理士】、【氏名又は名称】江口 昭彦
【識別番号】100134120、【弁理士】、【氏名又は名称】内藤 和彦
【識別番号】100109586、【弁理士】、【氏名又は名称】土屋 徹雄
【識別番号】100083806、【弁理士】、【氏名又は名称】三好 秀和
審査官 【審査官】伊藤 幸仙
参考文献・文献 特開平09-047438(JP,A)
特開2004-081657(JP,A)
Ion-Florin Talos, et al.,'Diffusion Tensor and Functional MRI Fusion with Anatomical MRI for Image-Guided Neurosurgery',MICCAI 2003,ドイツ,Springer,2003年11月15日,LNCS vol. 2878,p407-p415
月本洋 他4名,論理回帰分析法によるfMRI画像の解析,電子情報通信学会技術研究報告,日本,社団法人電子情報通信学会,2004年 5月20日,Vol.104,No.88,p59-p64
調査した分野 A61B 5/055
JSTPlus(JDreamII)
JMEDPlus(JDreamII)
JST7580(JDreamII)
特許請求の範囲 【請求項1】
ボクセル単位で取得された、脳内のプロトンの拡散度を特定可能とする拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成
前記隣接ボクセル間の結合度の評価量に基づき、前記ボクセル単位で取得された、脳の活性部位を特定可能とする脳機能データの解析を行う、
ことを含む、脳機能データ解析方法。
【請求項2】
前記脳機能データの解析を行うことは、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析により、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを含む、請求項1記載の脳機能データ解析方法。
【請求項3】
前記脳機能データの解析を行うことは、前記脳機能データ検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを含む、請求項1記載の脳機能データ解析方法。
【請求項4】
前記脳機能データの解析を行うことは、前記脳機能データ分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを含む、請求項1記載の脳機能データ解析方法。
【請求項5】
前記脳機能データの解析を行うことは、前記脳機能データと所定のモデルとの相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを含む、請求項1記載の脳機能データ解析方法。
【請求項6】
前記脳機能データの解析を行うことは、前記脳機能データから主要な成分を抽出する抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを含む、請求項1記載の脳機能データ解析方法。
【請求項7】
前記脳機能データの解析を行うことは、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化することを含む、請求項1記載の脳機能データ解析方法。
【請求項8】
前記脳機能データの解析を行うことは、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングすることを含む、請求項1記載の脳機能データ解析方法。
【請求項9】
前記拡散テンソルデータに基づき、前記脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施すことをさらに含む、請求項1記載の脳機能データ解析方法。
【請求項10】
コンピュータを、
脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、
前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、
前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、
前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段と、
して機能させることを特徴とする脳機能解析プログラム。
【請求項11】
請求項10に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析する手段であって、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを特徴とする脳機能解析プログラム。
【請求項12】
請求項10に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データを検定する手段であって、前記検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項13】
請求項10に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データを分類する手段であって、前記分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項14】
請求項10に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データと所定のモデルとの相関をとる手段であって、前記相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項15】
請求項10に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記脳機能データから主要な成分を抽出する手段であって、前記抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析プログラム。
【請求項16】
請求項10に記載の脳機能解析プログラムにおいて、前記データ解析は、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化することを特徴とする脳機能解析プログラム。
【請求項17】
請求項10に記載の脳機能解析プログラムにおいて、前記データ解析手段は、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングすることを特徴とする脳機能解析プログラム。
【請求項18】
請求項10乃至17のいずれかに記載の脳機能解析プログラムにおいて、前記コンピュータを、前記拡散テンソルデータ取得手段によって取得した拡散テンソルデータに基づき、前記脳機能データ取得手段によって取得した脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施すデータ前処理手段として機能させることを特徴とする脳機能解析プログラム。
【請求項19】
脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、
前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、
前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、
前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段と、
を備えたことを特徴とする脳機能解析装置。
【請求項20】
請求項19に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データ及びタスクの一方を説明変数とし、前記脳機能データ及び前記タスクの他方を被説明変数とする回帰分析する手段であって、前記脳機能データのボクセル毎の評価量に前記隣接ボクセル間の結合度の評価量を組み込んだ評価量の最小値及び最大値の一方を算出することを特徴とする脳機能解析装置。
【請求項21】
請求項19に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データを検定する手段であって、前記検定の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
【請求項22】
請求項19に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データを分類する手段であって、前記分類の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
【請求項23】
請求項19に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データと所定のモデルとの相関をとる手段であって、前記相関の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
【請求項24】
請求項19に記載の脳機能解析装置において、前記データ解析手段は、前記脳機能データから主要な成分を抽出する手段であって、前記抽出の基準値を前記隣接ボクセル間の結合度の評価量に基づき調整することを特徴とする脳機能解析装置。
【請求項25】
請求項19に記載の脳機能解析装置において、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データを平滑化するデータ平滑化手段をさらに備えたことを特徴とする脳機能解析装置。
【請求項26】
請求項19に記載の脳機能解析装置において、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データをクラスタリングするデータクラスタリング手段をさらに備えたことを特徴とする脳機能解析装置。
【請求項27】
請求項19乃至26のいずれかに記載の脳機能解析装置において、前記拡散テンソルデータ取得手段によって取得した拡散テンソルデータに基づき、前記脳機能データ取得手段によって取得した脳機能データから特定される活性部位ボクセルの脳機能データの値を他のボクセルの脳機能データの値に伝播させる前処理を施すデータ前処理手段をさらに備えたことを特徴とする脳機能解析装置。
発明の詳細な説明 【技術分野】
【0001】
本発明は、fMRI(functional Magnetic Resonance Imaging:機能的磁気共鳴撮影)やPET(Positron Emission Tomography:陽電子断層撮影)等の計測手法により得られる各種のデータを用いて脳機能を解析するための脳機能データ解析方法、脳機能解析装置及び脳機能解析プログラムに関する。
【背景技術】
【0002】
現在、代表的な脳の非侵襲計測方法としては、fMRIやPET、MEG(magnetoencephalography:脳磁図)等があるが、このうちデータの空間的分解能が最も高いとされ、広く使用されているのがfMRIである。
【0003】
fMRIは、脳の活性部位を特定可能とする各種の物理量を測定量としてイメージングするものであり、脳機能の計測に対して有効な手法である(非特許文献1及び2を参照)。fMRIは、脳の構造をイメージングする解剖MRIの原理と同じく、生体内組織のプロトン密度や縦緩和時間T1、横緩和時間T2を反映するものであるが、特に、脳の活性部位における血流量の増加を捉える点に特徴を有する。脳の活性部位においては局所的に血流量が増加することが知られており、血液中のヘモグロビンは、酸素が結合した状態(酸素化ヘモグロビン)と離れた状態(脱酸素化ヘモグロビン)とでその磁気的性質が異なる。増加した動脈血では磁場を乱す脱酸素化ヘモグロビンの量が少なくなるため、活性部位のfMRI信号(BOLD信号)が増加すると考えられる。したがって、fMRIを用いれば、被験者があるタスクを行っているときのBOLD信号の変化を手がかりにして、そのタスクに関係した脳の部位(活性部位)を特定することができる。
【0004】
fMRIにより計測されたBOLD信号の時系列データを解析する際の代表的な手法としては、一般線形モデルに基づいたSPM(Statistical Parametric Mapping)(非特許文献1を参照)や主成分分析、独立成分分析に基づいたデータ解析(非特許文献3を参照)等がある。これらの手法の特徴は、BOLD信号の時系列データを3次元画素(Voxel:ボクセル)毎に個別に統計処理した結果を画像として出力し、脳の活性部位を特定することにある。
【0005】
しかしながら、上記の解析手法においては、データ解析の際に、神経のネットワーク構造が考慮されていないという課題がある。脳は多数の神経細胞がシナプスを介して複雑なネットワークを構成しており、近年の脳科学の知見によれば、脳はこのような神経ネットワークを介して各部位が互いに連携することで全体として一つの高次脳機能を果たしていると考えられている。例えば、被験者があるタスクを行う際に、脳の複数の部位が活性化するといった現象が観測されるが、この現象を上記の解析手法によって解析すると、各活性部位を特定することは可能となるが、活性部位間の結合を特定することは困難である。
【0006】
これにはそもそも以下の理由がある。つまり、上述したようにfMRIにより計測されるBOLD信号は血流量に基づいているため、脳内において血流量が相対的に多い灰白質(神経細胞の細胞体)の活性は捉えられるが、血流量が相対的に少ない白質(神経細胞の軸索、又は神経線維)の活性はfMRIでは捉えにくいからである。
【0007】
一方において、神経ネットワーク構造の基盤となるこのような神経線維群の走行方向を捉えるために、近年、MRIによる新たな観測量として生体組織内のプロトンの拡散の度合いを計測するDTI(Diffusion Tensor Imaging:拡散テンソル撮像)法が注目を集めている(非特許文献4を参照)。通常の解剖MRIを利用する場合、神経線維はT1強調画像で高信号、T2強調画像で低信号となる。神経線維がT1強調画像で高信号化する理由はミエリンの存在にあり、ミエリンは二層構造の脂質や巨大タンパク質から成り、神経線維の走行方向に沿って配列するという形態をとる。したがって、神経線維の走行方向ではプロトンの拡散定数が大きく、それに直交する方向では拡散定数が小さいといった異方性が生じる。DTIは、プロトンの拡散を強調するためのMPG(Motion Probing Gradient:傾斜磁場):
【数1】
JP0004308871B2_000002t.gif

【0008】
を印加することによって拡散の異方性を計測する手法であり、例えば、SE(Spin Echo)パルスの前後に拡散検出用のSTG(Stejskal-Tanner Gradient)パルスを加えたST(Stejskal-Tanner)パルス系列によって得られるBOLD信号の強度S’(l,m,k,i)を計測するものである。ここで、l,m,kはボクセルの位置を表す正の離散変数であり、それぞれ、ボクセルのX方向、Y方向、Z方向の座標を表す。また、iは、測定時間を表す正の離散変数である。該BOLD信号の強度S’(l,m,k,i)は、
【数2】
JP0004308871B2_000003t.gif

【0009】
と書くことができ、拡散強調画像を生成する際には式(1)における拡散テンソル
【数3】
JP0004308871B2_000004t.gif

【0010】
がデータ解析の対象となる。ここで、ρ’(l,m,k,i)はMPGを印加しない場合のBOLD信号(通常の脳機能解析におけるデータ解析の対象)の強度を表し、bはMPGの強さを表すパラメータである。なお、ρ’(l,m,k,i)は、
【数4】
JP0004308871B2_000005t.gif

【0011】
と表される。ここで、f(ν)は流速、TRは繰り返し時間、TEはエコー時間、そしてξ’(l,m,k,i)はプロトン密度を表す。
【0012】
現在、fMRIにより計測されて得られたBOLD信号の時系列データρ’(l,m,k,i)を解析した後に、脳の関心領域(Region of Interest:ROI)を拡散テンソルデータD(l,m,k)によって接続するといった脳機能解析方法が用いられつつある(非特許文献5を参照)。

【非特許文献1】“Human Brain Function:2nd-Ed.”, Richard S. J. Frackowiak, et al, ELSEVIER ACADEMIC PRESS, 2004
【非特許文献2】“Image of Mind”, M. I. Posner and M. E. Raichle, W H Freeman & Co, 1997
【非特許文献3】“Independent Component Analysis: Theory and Applications”, T.-W. Lee, Kluwer Acadmic, 1988
【非特許文献4】“これでわかる拡散MRI”青木 茂樹, 阿部 修,秀潤社,2002
【非特許文献5】“Combined functional MRI and tractography to demonstrate the connectivity of the human primary motor cortex in vivo”, Guye M,et al., Neuroimage, Vol.19, pp.1349-1360, 2003
【発明の開示】
【0013】
しかしながら、上記の脳機能解析方法においては、BOLD信号の時系列データρ’(l,m,k,i)自体の解析には拡散テンソルデータD(l,m,k)が用いられていないため、その有効性には疑問が残る。
【0014】
本発明は、上記の課題を鑑み為されたものであり、その目的とするところは、fMRIやPET等の非侵襲計測方法を用いて計測された脳機能データと神経線維群の走行方向を特定可能とする拡散テンソルデータとに基づき脳の活性部位間の結合構造をも考慮した脳機能データの解析を可能とする脳機能解析方法および脳機能解析プログラムを提供することにある。
【0015】
本発明の一つの側面は、ボクセル単位で取得された、脳内のプロトンの拡散度を特定可能とする拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成
前記隣接ボクセル間の結合度の評価量に基づき、前記ボクセル単位で取得された、脳の活性部位を特定可能とする脳機能データの解析を行う、
ことを含む、脳機能データ解析方法を提供することである。
本発明の他の側面は、脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段とを備えたことを特徴とする脳機能解析装置を提供することである。
【0016】
本発明の他の側面は、コンピュータを、脳の活性部位を特定可能とする脳機能データをボクセル単位で取得する脳機能データ取得手段と、前記脳内のプロトンの拡散度を特定可能とする拡散テンソルデータをボクセル単位で取得する拡散テンソルデータ取得手段と、前記拡散テンソルデータに基づき隣接ボクセル間の結合度の評価量を構成するデータ評価量構成手段と、前記隣接ボクセル間の結合度の評価量に基づき前記脳機能データの解析を行うデータ解析手段と、して機能させることを特徴とする脳機能解析プログラムを提供することである。
【図面の簡単な説明】
【0017】
【図1】図1は、本発明の第1の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。
【図2】図2は、典型的なfMRI計測の方法を例示した図である。
【図3】図3は、典型的なfMRI画像(2次元スライス画像)を例示した図であり、図3Aは寝台上の被験者の頭部を表し、図3Bは頭部の2次元断面(2次元スライス画像)を構成するボクセルを表す。
【図4】図4は、2次元スライス画像を構成するボクセル上の脳機能情報の時系列データと拡散テンソルデータを例示した表である。
【図5】図5は、図4に示した脳機能情報の時系列データの一例を示した表であり、図5Aは実測値を表し、図5Bは実測値を閾値に基づき二値化した値を表す。
【図6】図6は、図4に示した脳機能情報の時系列データと拡散テンソルデータとを立体構成した図であり、図6Aは寝台上の被験者の頭部を表し、図6Bは図4に示した脳機能情報の時系列データと拡散テンソルデータとからk番目の2次元スライス画像Skを生成する手順を表す。
【図7】図7は、図6に示した2次元スライス画像Skから3次元構成された頭部画像を例示した図であり、図7Aは寝台上の被験者の頭部を表し、図7Bは2次元スライス画像Skを集めて3次元的に構成された頭部画像を表す。
【図8】図8は、図1に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図9】図9は、本発明に係るデータ前処理の手法の一例を示した図であり、図9Aはあるボクセルに隣接する2つのボクセルが共に灰白質と見なされる領域におけるデータ前処理を表し、図9Bはこれら2つのボクセルのうちの1つだけが灰白質と見なされる領域におけるデータ前処理を表し、図9Cはこれら2つのボクセルが共に灰白質とは見なされない領域におけるデータ前処理を表す。
【図10】図10は、本発明の第2の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。
【図11】図11は、図10に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図12】図12は、図10に示したデータ平滑化手段によって行われる平滑化の手法を例示した図である。
【図13】図13は、第2の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。
【図14】図14は、図13に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図15】図15は、本発明の第3の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。
【図16】図16は、図15に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図17】図17は、図15に示したクラスタリング手段300によって行われるクラスタリングの手法を例示した図であり、図17Aはクラスタリング処理を行う前の結合度ベクトルの値を表し、図17Bはクラスタリング処理を行った後の結合ベクトルの値を表す。
【図18】図18は、第3の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。
【図19】図19は、図18に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図20】図20は、本発明の第4の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。
【図21】図21は、図20に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【図22】図22は、第4の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。
【図23】図23は、図22に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャ-トである。
【図24】図24は、被験者に聴覚刺激を用い単純復唱を行わせた場合の各種脳機能解析を行った結果を示した図であり、図24A,24Bは共にSPMを用いたT検定(図24AはT検定の閾値を修正した場合,図24BはT検定の閾値を無修正の場合)の結果を表し、図24CはSPMとトラクトグラフィーとを組み合わせた脳機能解析手法を用いた解析結果を表し、図24Dは本発明の第1の実施形態に係る脳機能解析方法を用いた解析結果を表す。

【発明を実施するための最良の形態】
【0018】
本発明は、脳の活性部位を特定するために、MRI等の非侵襲計測装置から脳機能データと拡散テンソルデータとをボクセル単位で取得し、該拡散テンソルデータから隣接ボクセル間の結合度の評価量を構成し、該隣接ボクセル間の結合度の評価量を用いて、該脳機能データの解析を行うことを特徴とするものである。上記脳機能データの隣接ボクセル間の結合度の評価量はデータ解析の手法に依存して様々考えられるし、また該隣接ボクセル間の結合度の評価量をどのように用いるかについてもデータ解析の手法に応じて様々考えられる。
【0019】
そこで、以下に、本発明の実施形態を、図面を用いて詳細に説明する。
【0020】
[第1の実施形態]
図1は、本発明の第1の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第1の実施形態の脳機能解析装置10は、MRI装置50で計測された脳機能情報の原時系列データρ’(l,m,k,i)(例えば、(3)式)を取得する脳機能データ取得手段1と、同じくMRI装置50で計測された拡散テンソルデータD(l,m,k)(例えば、(2)式)を取得する拡散テンソルデータ取得手段2と、脳機能データ取得手段1で取得した脳機能情報の原時系列データρ’(l,m,k,i)に対して前処理を施すデータ前処理手段3と、拡散テンソルデータ取得手段2で取得された拡散テンソルデータD(l,m,k)から隣接ボクセル間の結合度を表す結合度ベクトル
【数5】
JP0004308871B2_000006t.gif

【0021】
を算出するボクセル間結合度算出手段4と、データ前処理手段3で前処理された脳機能情報の時系列データρ(l,m,k,i)とボクセル間結合度算出手段4で算出された結合度ベクトル
【数6】
JP0004308871B2_000007t.gif

【0022】
とから所定の評価量Qを構成するデータ評価量構成手段5と、該データ評価量Qの極値(本実施形態では最小値)を求めるための計算を行うデータ解析手段6と、データ解析手段6で算出されたデータに基づき画像を生成する画像生成手段7と、画像生成手段7によって生成された画像を表示する画像表示手段8と、上記各手段2~8で取得、算出、構成、解析、及び生成された各種データを記録するサブメモリと、以降で詳述する各ステップを実行するためのコンピュータに読み取り可能な脳機能解析用プログラムを記憶するメインメモリとから成る記憶手段9と、記憶手段9から読み出された該脳機能解析用プログラムに従って上記各手段2~9を制御するCPU(Central Processing Unit:中央演算処理装置)等を備える。なお、上記の各記号の意味は後に説明する。
【0023】
ここで、MRI装置50は、核磁気共鳴(Nuclear Magnetic Resonance:NMR)を起こすための静磁場を生成する静磁場コイル、共鳴周波数の高周波を照射し、共鳴信号を検出する高周波コイル、該共鳴信号に位置情報をエンコードするための傾斜磁場を生成する傾斜磁場コイル等から成るマグネットアセンブリと、これらのコイルの通電を制御するシステムコントローラ等から構成されるものであり、操作者の要求に応じて各種fMRIデータ(血流量に関連する脳機能情報の原時系列データρ’(l,m,k,i)や神経繊維の走行方向に関連する拡散テンソルデータD(l,m,k)等)を生成し、これらのデータを脳機能解析装置10に送信するものである。また、本発明においては、脳機能データの原時系列データρ’(l,m,k,i)を必ずしもMRI装置50から取得する必要はなく、その代わりに、同じく各種の脳機能データを生成可能なPET等の非侵襲計測装置から取得する構成としてもよい。
【0024】
また、画像表示手段8は、例えば、CRTディスプレイ、TFTディスプレイ、プラズマディスプレイなどの各種ディスプレイや、インクジェットプリンタ、レーザープリンタなどの各種プリンタなどが使用可能である。また、記憶手段9は、例えば、RAM(Random Access Memory)やROM(Read Only Memory)等から構成されている。さらに、記録手段9のサブメモリとメインメモリとを別体として構成し、メインメモリ部分を磁気ハードディスク、フロッピー(登録商標)ディスク、CD-ROMなどの光ディスク、磁気テープ、メモリチップ等に記憶させてもよい。
【0025】
本実施形態においては、脳機能解析装置10を、本実施形態の主要構成手段10Aと画像表示手段8及び記憶手段9とを一体型としたディスプレイコンソール形式としたが、画像表示手段8、或いは、画像表示手段8及び記憶手段9をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段10Aから切り離した構成としてもよい。いずれの構成においても、脳機能解析装置10はコンピュータによって実現されるものであり、上記各手段2~9は、CPU(Central Processing Unit:中央演算処理装置)100により記憶手段9から読み出された脳機能解析用プログラムに従って制御される。
【0026】
ここで、コンピュータとは、構造化された入力を所定の規則に従って処理し、処理した結果を構造化して出力する装置のことを指し、例えば、汎用コンピュータ、スーパーコンピュータ、メインフレーム、ワークステーション、マイクロコンピュータ、サーバ等が含まれる。また、通信ネットワーク(例えば、イントラネット、ローカルエリアネットワーク(LAN)、ワイドエリアネットワーク(WAN)、及びこれらの組み合わせから成る通信ネットワーク)を介して接続された2つ以上のコンピュータから成る構成(例えば、分散コンピュータシステム)であってもよい。
【0027】
次に、本発明の理解を容易にするために、図2乃至図7を用いて、本発明による脳機能解析の流れの概略を説明する。図2は、典型的なfMRI計測の方法を例示したものであり、図3は、典型的なfMRI画像(2次元スライス画像)を例示したものであり、図3Aは寝台上の被験者の頭部を表し、図3Bは頭部の2次元断面(2次元スライス画像)を構成するボクセルを表す。
【0028】
例えば、被験者にあるタスク(例えば、指のタッピング)を行わせた際に、そのタスクに関連すると考えられる脳の部位をfMRI計測によって特定することを考える。図2に示す例では、1回の計測において、一定時間Tのタスクと同じ一定時間Tの休憩(レスト)とを1セットとして3回繰り返す場合を示している。ここで、横軸は時間軸であり、“ON”はタスクを表し、“OFF”はレストを表す。通常、fMRI計測は1回の実験において数十回行われる。この例では、fMRI計測は24回行われるわけだが、タスクの立上がりと立下りの画像は計測時刻のずれを考慮して通常使用しないので、実際には、有意なfMRI計測は18回(図中の18本の太線に対応)行われたことになる。なお、ti(for i=1,2,…,24)は計測の実時間を離散的に示している。以降では、tiを単に“i(for i=1,2,…,I)”(正の離散整数)と表す。
【0029】
このようなfMRI計測の結果、図3Bに示すような2次元スライス画像Skが得られる。ここで、kは1回の実験によって得られる2次元スライス画像の枚数を表し、通常のfMRI計測においては20~30枚の2次元スライス画像が得られる(図中、K=20~30)。これらの2次元スライス画像Sk(for k=1, 2,…,K)は、図3Bにおいては平面的に記載されているが、実際にはボクセルと呼ばれる3次元の立体画素から構成されている。同図では、一例として64×64個のボクセルから成る2次元スライス画像を示したが、もちろん、その他の個数のボクセルに分割した2次元スライス画像であってもよい。一般に、任意の2次元スライス画像Sk上のボクセルは2つの正の離散変数l,mを用いて(l,m)と表すことができるが、所定の規則で番号付けすることによってこれを単一の正の離散変数jで表すこともできる。つまり、数学的には、jと(l,m)とを該所定の規則によって1対1対応(j≡(l,m)(for Sk))させることが可能である。この例では、lの上限L及びmの上限Mが共に64(L=64,M=64)であることから、j=1,2,…,4096となる。以降においては、これらの表記を適宜使い分けることにする。また、図3Bの丸は頭部の断面形状の輪郭を表している。
【0030】
図4は、脳機能データ取得手段1で取得し、データ前処理部3で前処理(後述)を施した、ある2次元スライス画像Skを構成するボクセルj上の脳機能情報の時系列データρ(j,k,i)(≡ρk(j,i))と拡散テンソルデータ取得手段2で取得されたD (j,k)(≡Dk(j))を例示したものである。背景技術で述べたように、DTI計測によって拡散強調画像を生成する際には、典型的にはSEパルスの前後に拡散検出用のSTGパルスを加えたSTパルス系列でイメージングした結果として、該2次元スライス画像Sk上の脳機能情報の原時系列データρ’(l,m,k,i)(≡ρ’k(l,m,i))と拡散テンソルデータD(l,m,k) (≡Dk(l,m))とが得られる。ρk(j,i)は、この原時系列データρ’k(l,m,i)に対して後述する所定の前処理を施した結果得られる脳機能情報の時系列データを表す。つまり、図中のρk(j,i)は、ある2次元スライス画像Skのボクセルjにおける時刻iの脳機能情報を表す。また、Dk(j)は、ある2次元スライス画像Skのボクセルjにおける拡散テンソル情報を表す。ρk(j,i)は具体的にはスカラーとして表され、また、Dk(j)は式(2)に示したように行列で表現され、例えば、
【数7】
JP0004308871B2_000008t.gif

【0031】
のような値をとる。この例で示すように、拡散テンソルは通常、対称テンソル(非対角成分が対称)なので、実質的には6個の成分を有することになる。また、図中のCLASSとは、タスクを行っているか否かを表したものであり、“ON(=1)”及び“OFF(=0)”はそれぞれ図2の“ON”及び“OFF”に対応する。
【0032】
図5は、図4に示した脳機能情報の時系列データρk(j,i)の一例を示した表である。図5(A)に例示するように、実際の脳機能情報の時系列データρk(j,i)は連続値を取るが、例えば図5(B)に示すように、値“4000”を活性の閾値として設定して、それよりも大きい値を取るボクセル値は活性“A”、小さい値を取るボクセル値は不活性“I”というように前処理の段階でさらに二値化してもよい。
【0033】
図6は、図4に示した拡散テンソルデータDk(j)を用いて、同じく図4で示した脳機能情報の時系列データρk(j,i)とを立体構成した図であり、図6Aは寝台上の被験者の頭部を表し、図6Bは図4に示した脳機能情報の時系列データと拡散テンソルデータとからk番目の2次元スライス画像Skを生成する手順を表す。図6Bのボクセルj1及びj2上の5角形は、図4に示した脳機能情報の時系列データρk(j,i)及び拡散テンソルデータDk(j)が付随した各ボクセル上の内部空間(つまり、ボクセル上のスカラー場及びテンソル場)を例示したものである。
【0034】
例えば、背景技術で述べたSPMにおいては、MRI装置50から取得した脳機能情報の原時系列データρ’k(l,m,i)に対して、
(a)整列(Realignment):先頭の2次元スライス画像(例えば、S1)に、以降の画像をそろえることで、測定中の頭の動きに伴う位置補正をし、動きによる擬似信号を取り除き、
(b)空間的正規化(Spatial Normalization):複数の被験者のデータをまとめたり比較したりするために、各被験者のデータをTalairach標準脳に合わせ、
(c)スムージング(Smoothing):ノイズが多く含まれる原時系列データにガウス型フィルターを適用することで、空間的分解能を下げることなく、解析の感度を向上させる、
といった前処理を施した後に、ボクセル毎に各種の検定等を行い、タスクに伴う有意な活性部位を、被験者の個人データで、さらには複数の被験者から成るグループ間で検定する。さらに、SPMでは、活動推定のモデルを事前に用意し、上記のようにして前処理された機能情報の時系列データρk (j,i)(=ρk (l,m,i))がどれだけモデルに照合するかを、拡散テンソルデータDk(l,m) を用いることなく、一般線形モデルを使いパラメータ推定する。
【0035】
一方、本発明では、脳機能情報の原時系列データρ’k(l,m,i)に対して、拡散テンソルデータ取得手段2によって取得された拡散テンソルデータDk(j)(=Dk(l,m))を用いて、上記の前処理以外にさらに以下で詳述するような前処理(d)を施す。そして、該拡散テンソルデータDk(j)(=Dk(l,m))を用いて、前処理された脳機能情報の時系列データρk(j,i)(=ρk(l,m,i))に以下で詳述するようなデータ処理を施すことによって画像表示の対象となる2次元画像データak(j)(=ak(l,m))を生成する。
【0036】
なお、上記では簡単のために2次元スライス画像を例にとり説明したが、実際には、図7A,7Bに示すような頭部の立体画像データak(j)(=a(j,k)=a(n))を3次元的にデータ解析する。その際にも、以降で詳述するような、拡散テンソルデータDk(j)を用いたデータ処理が施される。ここで、“n”とは、立体構成された全ボクセルに対して所定の規則で番号付けしたボクセルの位置を表す正の離散変数である。図中では、XY座標上の2次元スライス画像の位置を表す変数jと3次元目の方向(z方向)を表す変数kとを分けて記載したが、これを上記した1変数“n”で表すこととは該所定の規則の下で数学的に同値である(n≡(j,k))。なお、以降においては、同一の3次元ボクセルを表すのに1変数“n”を用いた表記とXYZ座標系における3変数(j,l,k)を用いた表記とを適宜使い分ける(n≡(l,m,k)≡(j,k))。
【0037】
以上の前提の下、図1に示した脳機能解析装置10によって行われるデータ解析方法について詳細に説明する。図8は、図1に示した脳機能解析装置10によって行われるデータ解析の処理手順を示したフローチャートである。
【0038】
ステップS10において、脳機能データ取得手段1は、MRI装置50で計測された脳機能情報の原時系列データρ’(l,m,k,i)を取得する。同ステップにおいて、拡散テンソル情報取得手段2は、同じくMRI装置50で計測された拡散テンソルデータD(l,m,k)を取得する。
【0039】
次に、ステップS20において、データ前処理手段3は、上記のようにして取得した脳機能情報の原時系列データρ’(l,m,k,i)に対して、SPMの上記手法(a)~(c)に代表されるような位置補正、ノイズ除去等の前処理を施す。
【0040】
同ステップにおいて、データ前処理手段3は、上記のようにして通常の前処理が施された脳機能情報の時系列データρ’(l,m,k,i)に対してさらに前処理(d)(データの神経信号化)を施すために、上記のようにして取得した拡散テンソルデータD(l,m,k)を用いて、例えば、拡散テンソルデータD(l,m,k)の固有方程式
【数8】
JP0004308871B2_000009t.gif

【0041】
を満たす固有ベクトル
【数9】
JP0004308871B2_000010t.gif

【0042】
の中で最大固有値λMに対応する固有ベクトル
【数10】
JP0004308871B2_000011t.gif

【0043】
を算出する(α=1,2,3)。説明を簡単にするために、以降においては、|λ1|≧|λ2|≧|λ3|とし、最大固有値であるλ1に対応する固有ベクトル
【数11】
JP0004308871B2_000012t.gif

【0044】

【数12】
JP0004308871B2_000013t.gif

【0045】
と書くことにする。また、該固有ベクトル
【数13】
JP0004308871B2_000014t.gif

【0046】
のノルムは1に規格化されているものとする。
【0047】
次いで、同ステップにおいて、データ前処理手段3は、まず、あるボクセル(l,m,k)に対して、該ボクセル(l,m,k)と隣接する26個のボクセル(3次元)の中から該ボクセル(l,m,k)上の固有ベクトル
【数14】
JP0004308871B2_000015t.gif

【0048】
との内積の絶対値が最大となる方向ベクトルに対応する隣接ボクセルを求める。この条件を満たす隣接ボクセルを、例えば、(l+1,m+1,k)とする。そして、このようにして求められたボクセルと該ボクセル(l,m,k)を中心として点対称の関係にある隣接ボクセルを求める。この場合は、ボクセル(l-1,m-1,k)である。この状況を図9に示す。上述したように、ボクセル(l+1,m+1,k)及び(l-1,m-1,k)は、ボクセル(l,m,k)上の固有ベクトル
【数15】
JP0004308871B2_000016t.gif

【0049】
との内積が最大となるボクセルであるから、ボクセル(l,m,k)との結合が最も強いと考えられる。なぜならば、図9(c)の場合であれば、ボクセル(l,m,k)を中心としてボクセル(l+1,m+1,k)及び(l-1,m-1,k)は同じ走行方向に走る神経線維上のボクセルと見なすことができるからである。
【0050】
次に、同ステップにおいて、データ前処理手段3は、以下の3通りの状況に応じて、ボクセル(l,m,k)上の脳機能情報の時系列データρ’(l,m,k,i)の書き換えを行う。
【0051】
(状況1)予め取得された所定の脳の解剖学的データ(MRI装置50により取得されたT1強調画像)から、ボクセル(l+1,m+1,k)及び(l-1,m-1,k)が共に灰白質(神経細胞の細胞体:図中のG)と見なされる場合(図9(A)を参照):
この場合には、ボクセル(l,m,k)上の脳機能情報の時系列データρ’(l,m,k,i)の値を、
【数16】
JP0004308871B2_000017t.gif

【0052】
と書き換える。
【0053】
(状況2)予め取得された所定の脳の解剖学的データ(MRI装置50により取得されたT1強調画像)から、ボクセル(l+1,m+1,k)(又は、ボクセル(l-1,m-1,k))が灰白質と見なされる場合(図9(B)を参照):
この場合には、ボクセル(l,m,k)上の脳機能情報の時系列データρ’(l,m,k,i)の値を、
【数17】
JP0004308871B2_000018t.gif

【0054】
と書き換える。
【0055】
(状況3)予め取得された所定の脳の解剖学的データ(MRI装置50により取得されたT1強調画像)から、ボクセル(l+1,m+1,k)及び(l-1,m-1,k)が共に灰白質とは見なされない(つまり、白質(神経の神経繊維:図中のW)と見なされる)場合(図9(C)を参照):
この場合は、ボクセル(l,m,k)上の脳機能情報の時系列データρ’(l,m,k,i)の値を、
【数18】
JP0004308871B2_000019t.gif

【0056】
と書き換える。
【0057】
このような処理を全ボクセルに対して行う。その結果、図9(A),(B)の下図に示したように、灰白質と見なされるボクセル(灰白質ボクセル)上の脳機能情報の時系列データの値を白質と見なされる隣接ボクセル(白質ボクセル)上の脳機能情報の時系列データの値に反映(伝播)させることができる。また、図9(C)の下図に示したように、白質と見なされるボクセル同士が隣接している場合には、これらのボクセル上の脳機能データの値が平滑化されることになる。このような前処理を施すことによって、以降のステップの効果を際立たせることができる。
【0058】
なお、このような前処理(d)の仕方は、上記した例に限らず、例えば、灰白質ボクセルに隣接し、且つ該灰白質ボクセルの脳機能情報の時系列データの値を受け取った白質ボクセルが、さらに、隣接する白質ボクセルに対して上記と同様の手順によってその値を順次伝播させるようにしてもよい。また、上記のように灰白質ボクセルの脳機能情報の時系列データの値を白質ボクセルに渡す際に、拡散テンソルデータD(l,m,k)に応じて重み付けをした値を渡すような構成にしてもよい。さらに、上記の手順で書き換えられた白質ボクセルの値と元々の白質ボクセルの値のうちの大きい方を白質ボクセルの値としてもよい。
【0059】
次に、ステップS30において、ボクセル間結合度算出手段4は、上記の拡散テンソルデータ
【数19】
JP0004308871B2_000020t.gif

【0060】
に基づき、隣接ボクセル間の結合度ベクトル
【数20】
JP0004308871B2_000021t.gif

【0061】
を、例えば、
【数21】
JP0004308871B2_000022t.gif

【0062】
から算出する。ここで、
【数22】
JP0004308871B2_000023t.gif

【0063】
の各成分は、例えば、
【数23】
JP0004308871B2_000024t.gif

【0064】
と定義されるのものであり、以降において、これを隣接ボクセル間の平均拡散度ベクトルと称する。αは
【数24】
JP0004308871B2_000025t.gif

【0065】
を満たすパラメータであり、
【数25】
JP0004308871B2_000026t.gif

【0066】

【数26】
JP0004308871B2_000027t.gif

【0067】
を満たす定数ベクトルである。
【0068】
式(13)中の
【数27】
JP0004308871B2_000028t.gif

【0069】
は、プロトンの拡散度を表すベクトル(以降、これを拡散度ベクトルと称する)であり、第1成分Dl(l,m,k)は、位置(l,m,k)のボクセル内のX軸方向へのプロトンの拡散度を表し、第2成分Dm(l,m,k)は、位置(l,m,k)のボクセル内のY方向の拡散度を表し、第3成分Dk(l,m,k)は、位置(l,m,k)のボクセル内のZ方向の拡散度を表す。ここで、X,Y,Z方向とは図7に示した方向と同一の方向を想定している。なお、式(9)においては、拡散テンソルが対称テンソルであることを前提としている。
【0070】
このような意味合いを有する拡散度ベクトル
【数28】
JP0004308871B2_000029t.gif

【0071】
を用いて、式(13)を介して式(11)のように構成された隣接ボクセル間の結合度ベクトル
【数29】
JP0004308871B2_000030t.gif

【0072】
の各成分Cl(l,m,k),Cm(l,m,k),Ck(l,m,k)はそれぞれ、位置(l,m,k)のボクセルと位置(l+1,m,k)のボクセルとの結合の度合い、位置(l,m,k)のボクセルと位置(l,m+1,k)のボクセルとの結合の度合い、位置(l,m,k)のボクセルと位置(l,m,k+1)のボクセルとの結合の度合いを表す。なお、本実施形態においては、斜め方向に隣接するボクセル間(たとえば、ボクセル(l,m,k)とボクセル(l+1,m+1,k))の結合度は直接的には考慮していないが、直交する方向のボクセル(この場合、ボクセル(l+1,m,k)やボクセル(l,m+1,k))を介して間接的に考慮されている。もちろん、斜め方向に隣接するボクセル間の結合度を直接考慮するような構成にしてもよい。
【0073】
さらに、拡散度ベクトル
【数30】
JP0004308871B2_000031t.gif

【0074】
としては、例えば、
(I)楕円体モデルにおける
【数31】
JP0004308871B2_000032t.gif

【0075】
を採用する。ここで、|D|は、式(9)で表される行列の行列式を表す。また、例えば、
(II)式(9)で表される行列の対角成分Dll(l,m,k),Dmm(l,m,k),Dkk(l,m,k)を用いて、
【数32】
JP0004308871B2_000033t.gif

【0076】
とすることもできる。また、例えば、
(III)FA(Fractional Anisotropy):
【数33】
JP0004308871B2_000034t.gif

【0077】
やRA(Relational Anisotropy):
【数34】
JP0004308871B2_000035t.gif

【0078】
に代表される拡散の異方性を表す指標を用いて拡散度ベクトル
【数35】
JP0004308871B2_000036t.gif

【0079】
を決めることも可能である。ここで、DAVは、
【数36】
JP0004308871B2_000037t.gif

【0080】
である。その他にも、
(IV)予め取得された所定の脳の解剖学的データ(MRI装置50により取得されたT1強調画像)から脳の白質ボクセルと灰白質ボクセルとを判断し、白質ボクセルでは上記(I)の拡散度ベクトルを用い、灰白質ボクセルでは、その他の拡散度ベクトル(例えば、上記(II)又は(III))を用いることも可能である。
【0081】
このように、本発明における拡散度ベクトル
【数37】
JP0004308871B2_000038t.gif

【0082】
としては、MRI装置50から取得可能な拡散テンソルデータD(l,m,k)から神経線維の走行方向が特定可能な量であればいかなるものを用いてもよい。
【0083】
また、パラメータαに関する不等式(14)の代わりに、例えば、灰白質ボクセルで大きく、白質ボクセルで小さい値といったようにパラメータαの値を灰白質ボクセルと白質ボクセルとで切替えるように設定をしてもよい。
【0084】
また、式(11)の結合度ベクトル
【数38】
JP0004308871B2_000039t.gif

【0085】
の代わりに、
【数39】
JP0004308871B2_000040t.gif

【0086】
を用いてもよい。
【0087】
さらに、神経繊維の走行方向を強調するために、式(11)又は式(22)の結合度ベクトル
【数40】
JP0004308871B2_000041t.gif

【0088】
の各成分cl(l,m,k),cm(l,m,k),ck(l,m,k)をN乗(Nは2以上の自然数)したものを採用してもよい。また、これらの成分の値を白質ボクセルでa倍、灰白質ボクセルでb倍したものを用いてもよい。また、これらの成分の中で最も大きいもののみを例えば上記の方法で強調してもよい。さらに、{cl(l,m,k)+cm(l,m,k)+ck(l,m,k)}/3以上の値を有する成分のみを例えば上記の方法で強調してもよい。
【0089】
さらに、灰白質ボクセルは隣接する6個のボクセル(6方位)とのみ結合させ、白質ボクセルは固有方程式(5)の最大固有値に対応する固有ベクトルとの内積が最大のボクセル及びその点対称のボクセル(2方位)とのみ結合するように結合度ベクトル
【数41】
JP0004308871B2_000042t.gif

【0090】
を構成してもよい。
【0091】
次に、ステップS40において、データ評価量構成手段5は、ステップS20において前処理された脳機能情報の時系列データρ(l,m,k,i)とステップS40において算出された結合度ベクトル
【数42】
JP0004308871B2_000043t.gif

【0092】
とから、図5(A)に例示したようなデータをノンパラメトリック回帰分析するために、例えば、隣接ボクセルを6個とした場合には、データ評価量Q:
【数43】
JP0004308871B2_000044t.gif

【0093】
を構成する。ここで、iは時間、l,m,kはそれぞれボクセルのX方向、Y方向及びZ方向の位置を表す変数であることは上述した通りである。ただし、実際には、図7に示す立体構成されたボクセルの中には脳(神経)を表さないボクセルもあるので、その領域に関するl,m,kは除外するものとする。また、φ(i)は、図5(A)におけるCLASSを表す変数の時系列データであり、この場合、1(=ON)と0(=OFF)のニ値をとる。本実施形態に係るノンパラメトリック回帰分析では、図5(A)に示したような脳機能情報の時系列データρ(l,m,k,i)を説明変数、CLASSを表す変数の時系列データφ(i)を被説明変数とし、これらの変数の間に、
【数44】
JP0004308871B2_000045t.gif

【0094】
を回帰係数とした線形回帰式:
【数45】
JP0004308871B2_000046t.gif

【0095】
を想定する。
【数46】
JP0004308871B2_000047t.gif

【0096】
及び
【数47】
JP0004308871B2_000048t.gif

【0097】
はそれぞれ推定解を表し、κは隣接ボクセル間の連続性を表す重みパラメータである。このように、式(23)は隣接ボクセル間の連続性に基づいた評価量を構成したが、その代わりに、隣接ボクセル間の滑らかさに基づいた評価量を構成してもよい。
【0098】
本実施形態におけるデータ解析の手法としてノンパラメトリック回帰分析を用いる根拠は、被説明変数φ(i)の数が説明変数ρ(l,m,k,i)の数に比べて極端に少ないことにある。図5(A)に示す例では、説明変数ρ(l,m,k,i)はボクセルの数64×64×64=262144個あるのに対して、被説明変数φ(i)の測定値の数は18個である。通常でも、被説明変数の測定値の数は数10であり、100を越えることもある。このような状況においては通常の回帰分析ができないので、本実施形態ではノンパラメトリック回帰分析を用いる(非特許文献6:“Spline Smoothing and Nonparametric Regression”, R. L. Eubank, Marcel Dekker, Newyork, 1988)。
【0099】
通常、ノンパラメトリック回帰分析では、被説明変数φ(l,m,k,i)に連続性を仮定することによって、次の評価量Q’:
【数48】
JP0004308871B2_000049t.gif

【0100】
を最小にするような
【数49】
JP0004308871B2_000050t.gif

【0101】
を求める。ここで、右辺第1項は、各時刻iにおける被説明変数φ(i)の推定解
【数50】
JP0004308871B2_000051t.gif

【0102】
からの誤差を表す項であり、第2項が、被説明変数φ(i)の時間方向に対する連続性を評価するための項である。しかしながら、図5(A)に示す例からわかるように被被説明変数φ(i)の値に連続性を仮定することはできない。そこで、本実施形態では、このような連続性の制約を被説明変数φ(i)に対してではなく、説明変数ρ(l,m,k,i)に対して課すように拡張したノンパラメトリック回帰分析(非特許文献7:“脳機能画像のノンパラメトリック回帰分析”,月本 洋 他,電子情報通信学会誌D-II,Vol. J84, No.12,pp2623-2633,2001年12月)に基づいて、式(23)で表されるデータ評価量Qを構成する。
【0103】
ここで、式(23)の意味について説明する。右辺第1項は、式(25)と同じく、各時刻iにおける被説明変数φ(i)の推定解
【数51】
JP0004308871B2_000052t.gif

【0104】
からの誤差を表す項である。第2項~第4項は、隣接するボクセルが被説明変数φ(i)に与える影響の連続性を表す項であり、この連続性は隣接する説明変数ρ(l,m,k,i)が被説明変数に与える影響、すなわち、式(24)の回帰係数
【数52】
JP0004308871B2_000053t.gif

【0105】
の連続性、で評価できる。
【0106】
さらに、本実施形態においては、第2項~第4項において、ボクセル間結合度算出手段4で算出された結合度ベクトル
【数53】
JP0004308871B2_000054t.gif

【0107】
の各成分Cl(l,m,k),Cm(l,m,k),Ck(l,m,k)をそれぞれX方向、Y方向、Z方向に対する隣接ボクセル間の連続性に対する重み付け因子として導入している。式(11)に示したように、これらの成分はそれぞれ、平均拡散度ベクトル
【数54】
JP0004308871B2_000055t.gif

【0108】
の各成分D’l(l,m,k),D’m(l,m,k),D’k(l,m,k)の関数となっているので、結果として、式(23)のデータ評価量Qは神経線維(神経細胞)の走行方向(プロトンの拡散の異方性)を反映していることになる。
【0109】
次に、ステップS50において、データ解析手段6は、ステップS40において構成されたデータ評価量Qの極値を求めるための計算を行う。本実施形態においては、式(23)のデータ評価量Qを最小にするような回帰係数
【数55】
JP0004308871B2_000056t.gif

【0110】
を算出する。ここでは、簡単のため、式(23)を
【数56】
JP0004308871B2_000057t.gif

【0111】
のように書き換えたデータ評価量Q”について説明する。ここで、nは立体構成された全ボクセルに対して所定の番号付けを施すことによって1次元配列し直したボクセルの位置を表す正の離散変数であり、3変数表示(l,m,k)と1対1対応している。さらに、簡単のため、ここではL=M=Kとしたので、全ボクセル数はL3となる。また、C(n)は、このように1次元に配列し直した際に、ボクセルnと該ボクセルnと隣接するボクセルn’との間の結合度を表す係数であり、この係数は、式(11)の結合ベクトル
【数57】
JP0004308871B2_000058t.gif

【0112】
から一意に定まる。例えば、図3(B)に示した2次元スライス画像Skに話を限定した場合、n=130のボクセルと隣接するボクセルはn=127,129,131,257の4つとなる。ステップS30で算出された結合度ベクトル
【数58】
JP0004308871B2_000059t.gif

【0113】
により、ボクセルn=127とボクセルn=130との結合度、ボクセルn=129とボクセルn=130との結合度、ボクセルn=131とボクセルn=130との結合度、ボクセルn=257とボクセルn=130との結合度はそれぞれ、例えば、2,2,5,5のように決まる。つまり、この場合、
【数59】
JP0004308871B2_000060t.gif

【0114】
となる。
【0115】
このとき、回帰係数
【数60】
JP0004308871B2_000061t.gif

【0116】
は、最小2乗法を用いて、
【数61】
JP0004308871B2_000062t.gif

【0117】
のように求めることができる。ここで、Xはρ(n,i)を成分とするL3×I行列
【数62】
JP0004308871B2_000063t.gif

【0118】
はφ(i)を成分とするI次元ベクトル、
【数63】
JP0004308871B2_000064t.gif

【0119】

【数64】
JP0004308871B2_000065t.gif

【0120】
を成分とするL3次元ベクトル、κOは交差検証法(Cross Validation)で求めたκの最適値、すなわち、
【数65】
JP0004308871B2_000066t.gif

【0121】
を最小にするようなκである。ここで、
【数66】
JP0004308871B2_000067t.gif

【0122】
はi番目の事例を除いて回帰係数を計算し、その回帰係数を用いてi番目の事例の測定値を予測した値である。
【0123】
次に、ステップS60において、画像生成手段7は、上記のようにして算出された回帰係数
【数67】
JP0004308871B2_000068t.gif

【0124】
を、例えば、回帰係数の値が0(すなわち、無相関)のボクセルをグレイスケールの基準として、回帰係数の値が正(すなわち、正の相関)のボクセルを白方向の色調とし、回帰係数の値が負(すなわち、負の相関)のボクセルを黒方向の色調とするような色調を施した画像を生成する。その他にも、グレイスケール上で取り得るボクセルの値の頻度のヒストグラムにおいて、上位5%の値以上を表示させるようにしてもよいし、所定の規則に従ってフルカラー画像を生成してもよい。
【0125】
そして、ステップS70において、画像表示手段8は、上記のようにして生成された画像を立体表示する。
【0126】
このように、第1の実施形態においては、MRI装置から取得した脳機能情報の時系列データに基づく誤差の評価項に加えて、同じくMRI装置から取得した拡散テンソルデータに基づく空間方向の連続性の評価項を加えた評価量を構成した。そして、この評価項をノンパラメトリック線形回帰分析し、その結果得られた回帰係数を画像表示の対象としたことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。
【0127】
なお、本実施形態においては、ステップS20において通常の前処理(a)~(c)以外に前処理(d)のプロセスを設けたが、前処理(d)をはずしてもよい。このことは、以降の実施形態においても同様である。
【0128】
なお、第1の実施形態では、データ評価量として式(23)を構成したわけだが、本発明の目的を達成するにあたり、データ解析の手法に応じて、その他にも各種の評価量を構成することができる。そこで、次に、上記した第1の実施形態の幾つかの変形例について説明するが、重複を避けるため、以降においては第1の実施形態と異なる部分についてのみ説明を加える。
【0129】
<第1の実施形態の変形例1>
第1の実施形態では、ステップS50において、データ解析手段6は、式(19)の線形回帰式を用いることによりノンパラメトリック回帰分析をしたわけだが、より厳密には、回帰式を2次関数、3次関数等の一般のn次関数で行うことも考えられる。このような非線形回帰分析を行う際には、ニューラルネットワーク(Neural Network)モデルを用いることができる。典型的なニューラルネットワークモデルとして、入力層、中間層、出力層から成る階層型ニューラルネットワークモデルがある。
【0130】
本変形例では、ステップS40において、データ評価量構成手段5は、ニューラルネットワークモデルにおける通常の学習の評価量である2乗誤差の項に対して、ボクセル間結合度算出手段4で算出された結合度ベクトル
【数68】
JP0004308871B2_000069t.gif

【0131】
に基づく隣接ボクセル間の連続性を評価する評価量を加えた評価量を構成する。そして、ステップS50において、データ解析手段6は、このようにして構成された評価量の最小値を求めるために、誤差逆伝播法(Error Back Propagation Method)等を用いて解析した後に、ステップS60において、画像生成手段7は、入力(脳機能情報の時系列データ)と出力の間の感度(入力を変化させたときに出力に与える影響)を表示画像として第1の実施形態と同様にして生成する。そして、ステップS70において、画像表示手段8は、上記のようにして生成された画像を立体表示する。
【0132】
これによって、第1の実施形態と同様の効果が得られる。
【0133】
<第1の実施形態の変形例2>
第1の実施形態では、ステップS40において、データ評価量構成手段5は、脳機能の時系列データρ(l,m,k,i)を説明変数、φ(i)を被説明変数として式(23)の評価量Qを構成したわけだが、本変形例では、その逆、つまり、φ(i)を説明変数、ρ(l,m,k,i)を被説明変数として、ボクセル(l,m,k)毎に、例えば、
【数69】
JP0004308871B2_000070t.gif

【0134】
といった評価量Q”’を構成する。ここで、
【数70】
JP0004308871B2_000071t.gif

【0135】
である。
【0136】
そして、ステップS50において、データ解析手段6は、このようにして構成された評価量Q”’の最小値を求めるために、最適なκを交差検証法を用いて決定し、GA(Genetic Algorithm:遺伝的アルゴリズム)等の確率的探索アルゴリズムによって回帰係数
【数71】
JP0004308871B2_000072t.gif

【0137】
を求める。
【0138】
次いで、ステップS60において、画像生成手段7は、上記の解析結果を各種検定した際に有意であると判断されたボクセルの脳機能データρ(l,m,k)に基づく値を表示画像として生成する。そして、ステップS70において、画像表示手段8は、上記のようにして生成された画像を立体表示する。
【0139】
これによって、第1の実施形態と同様の効果が得られる。
【0140】
<第1の実施形態の変形例3>
第1の実施形態の変形例1では、ステップS40において、データ評価量構成手段5は、脳機能の時系列データρ(l,m,k,i)を説明変数、φ(i)を被説明変数として式(23)の評価量を構成したわけだが、本変形例では、その逆、つまり、φ(i)を説明変数、ρ(l,m,k,i)を被説明変数として、階層型ニューラルネットワークモデルを用いて、非線形回帰分析を行うこともできる。ただし、この場合には、ステップS50において極値を求める際に誤差逆伝播法が使えないので、その代わりに、データ解析手段6は、最適なκを交差検証法を用いて決定し、GA等の確率的探索アルゴリズムによってデータ解析を行う。そして、ステップS60において、画像生成手段7は、上記の解析結果を各種検定した際に有意であると判断されたボクセルの脳機能データρ(l,m,k)に基づく値を表示画像として生成する。そして、ステップS70において、画像表示手段8は、上記のようにして生成された画像を立体表示する。
【0141】
これによって、第1の実施形態と同様の効果が得られる。
【0142】
<第1の実施形態の変形例4>
第1の実施形態では、ステップS50において、データ解析手段6は、式(24)の線形回帰式を用いることによりノンパラメトリック回帰分析をしたわけだが、データ解析の手法として、判別分析(非特許文献8:“多変量解析のはなし”,有馬哲,石村貞夫, 東京図書,1987)、数量化II類(非特許文献9:“多変量解析のはなし”,有馬哲,石村貞夫, 東京図書,1987)、決定木(非特許文献10:“AIによるデータ解析”,J.R.キンラン,トッパン,1995)、サポートベクターマシン(非特許文献11:“サポートベクターマシン入門”,Nello Cristianini and John Shawe‐Taylor,共立出版,2005)等の分類的解析手法を用いることもできる。
【0143】
いずれの場合も、例えば、説明変数を脳機能情報の時系列データρ(l,m,k,i)、外的基準をタスクとレストの時系列データとして、各種の判別関数、決定木等を用いてデータ分類する。本変形例では、その際に、ステップS40において、データ評価量構成手段5は、各解析手法における通常の評価量に対して、ステップS30で算出した結合度ベクトル
【数72】
JP0004308871B2_000073t.gif

【0144】
を組み込んだ評価量を構成する。そして、ステップS50において、データ解析手段6は、各解析手法を用いてデータ解析を行う。上記のように評価量を変更すると、通常の解法が使えなくなる場合があるので、その際には、データ解析手段6は、最適なκを交差検証法を用いて決定し、GA等の確率的探索アルゴリズムによってデータ解析を行えばよい。
【0145】
そして、ステップS60において、画像生成手段7は、ステップS50で得られた分類データを各種検定した際に有意であると判断されたボクセルの脳機能データρ(l,m,k)に基づく値を表示画像として生成し、ステップS70において、画像表示手段8は、上記のようにして生成された画像を立体表示する。
【0146】
これによって、第1の実施形態と同様の効果が得られる。
【0147】
<第1の実施形態の変形例5>
第1の実施形態では、ステップS50において、データ解析手段6は、式(24)の線形回帰式を用いることによりノンパラメトリック回帰分析をしたわけだが、データ解析の手法として、ICA(Independent Component Analysis:独立成分分析)を用いることもできる(非特許文献12:“詳解独立成分分析 信号解析の新しい世界”, Aapo Hyvarinen et al.,東京電機大学出版局 2005)。
【0148】
本変形例では、ステップS40において、データ評価量構成手段5は、ボクセル間結合度算出手段4で算出された結合度ベクトル
【数73】
JP0004308871B2_000074t.gif

【0149】
に基づき評価量を構成する。独立成分分析による脳機能情報の時系列データρ(l,m,k,i)の解析での「空間的な独立」という条件を、本変形例では「拡散テンソル情報にもとづく結合度を考慮した空間的な独立」というように変更する。独立成分分析の具体的手法はいくつかあるが、例えば、相互相関を最小化する手法の場合には、本変形例では、相互相関を最小化するのではなく、結合度に基づく評価量に近づけるというようにする。これにより、結合度の大きいボクセル間の相互相関は結合度に基づく評価量に近いような「独立度」になり、結合度の小さい(もしくは、結合度のない)ボクセル間は相互相関を最小にするというようになる。より具体的には、例えば、相関行列の非対角成分の2乗和を最小化するような評価関数中の2乗和の各項に、結合度に基づく評価量による重み係数をつける。この重み係数は、結合度が大きいほど小さくなるように構成する。
【0150】
次に、ステップS50において、データ解析手段6は、上記のようにして構成された評価量を考慮した独立成分分析を行う。その際、最適なκを交差検証法を用いて決定し、GA等の確率的探索アルゴリズムによってデータ解析を行う。
【0151】
次に、ステップS60において、画像生成手段7は、上記のようにして得られた独立成分を表示画像として生成する。ただし、上記の解析の結果得られる独立成分の数は多いので、全行ベクトルの中からタスク/レストと最も相関の高い行ベクトルのみを抽出する処理を行う。
【0152】
そして、ステップS70において、画像表示手段8は、上記のようにして生成された画像を立体表示する。
【0153】
これによって、第1の実施形態と同様の効果が得られる。
【0154】
なお、上記した独立成分分析以外にも、脳機能情報の時系列データから所定の主たる成分を抽出する手法として、主成分分析、因子分析、数量化III類等があるが、これらの解析手法についても、基本的には、上記と同様にして解析することができる。
【0155】
[第2の実施形態]
図10は、本発明の第2の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第2の実施形態の脳機能解析装置20は、第1の実施形態の脳機能解析装置10におけるボクセル間結合度算出手段4、データ評価量構成手段5、及びデータ解析手段6の代わりに、ボクセル間結合度算出手段24、データ評価量構成手段25、データ平滑化手段200、及びデータ解析手段26を設けた構成を成している。第1の実施形態と同一の構成には同一の符号を付し、以降においては、繰り返しの説明を省略する。なお、本実施形態においては、脳機能解析装置20を、本実施形態の主要構成手段20Aと画像表示手段8及び記憶手段9とを一体型としたディスプレイコンソール形式としたが、画像表示手段8、或いは、画像表示手段8及び記憶手段9をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段20Aから切り離した構成としてもよい。
【0156】
ボクセル間結合度算出手段24は、拡散テンソルデータ取得手段2で取得された拡散テンソルデータD(l,m,k)(式(9))から第1の実施形態に詳述したようにして平均拡散度ベクトル
【数74】
JP0004308871B2_000075t.gif

【0157】
(式(12))を算出した後、該平均拡散度ベクトル
【数75】
JP0004308871B2_000076t.gif

【0158】
の各成分D’l(l,m,k),D’m(l,m,k),D’k(l,m,k)の大きさに応じてデータ前処理手段3で前処理された脳機能情報の時系列データρ(l,m,k,i)を平滑化する際の隣接ボクセル間の平滑度を決める結合度ベクトル
【数76】
JP0004308871B2_000077t.gif

【0159】
を算出する。
【0160】
データ平滑化手段200は、ボクセル間結合度算出手段24で算出された結合度ベクトル
【数77】
JP0004308871B2_000078t.gif

【0161】
を用いて、データ前処理手段3で前処理された脳機能情報の時系列データρ(l,m,k,i)を平滑化するための加重平均フィルターを構成し、該加重平均フィルターを用いて、データ前処理手段3で前処理された脳機能情報の時系列データρ(l,m,k,i)を平滑化する。
【0162】
データ解析手段26は、データ平滑化手段200で平滑化されたデータをSPM等の手法を用いてデータ解析を行う。
【0163】
図11は、図10に示した第2の実施形態に係る脳機能解析装置20によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【0164】
ステップS100において、脳機能データ取得手段1は、MRI装置50で計測された脳機能情報の原時系列データρ’(l,m,k,i)を取得する。同ステップにおいて、拡散テンソル情報取得手段2は、同じくMRI装置50で計測された拡散テンソルデータD(l,m,k)を取得する。
【0165】
次に、ステップS110において、データ前処理手段3は、ステップS100において取得した脳機能情報の原時系列データρ’(l,m,k,i)に対して、第1の実施形態で詳述した前処理(a)~(d)を施すことによって、脳機能情報の時系列データρ(l,m,k,i)を生成する。
【0166】
次に、ステップS120において、ボクセル間結合度算出手段24は、まず、ステップS100において取得した拡散テンソルデータD(l,m,k)から第1の実施形態で詳述したようにして平均拡散度ベクトル
【数78】
JP0004308871B2_000079t.gif

【0167】
を算出する。次いで、該平均拡散度ベクトル
【数79】
JP0004308871B2_000080t.gif

【0168】
の各成分D’l(l,m,k),D’m(l,m,k),D’k(l,m,k)の大きさに応じてデータ前処理手段3で前処理された脳機能情報の時系列データρ(l,m,k,i)を平滑化する際の隣接ボクセル間の平滑度を決める結合度ベクトル
【数80】
JP0004308871B2_000081t.gif

【0169】
を、例えば、
【数81】
JP0004308871B2_000082t.gif

【0170】
として算出する。
【0171】
次に、ステップS130において、データ平滑化手段26は、ステップS120において算出された結合度ベクトル
【数82】
JP0004308871B2_000083t.gif

【0172】
を用いて、ステップS110において前処理された脳機能情報の時系列データρ(l,m,k,i)を、例えば、
【数83】
JP0004308871B2_000084t.gif

【0173】
のように異方的に平滑化する。ここで、w及びWは重み因子であり、
【数84】
JP0004308871B2_000085t.gif

【0174】
を満たすように規格化する。
【0175】
図12は、本実施形態に係るデータ平滑化の手法を例示したものであり、ある2次元スライス画像Skの位置(l,m)を中心とした9つのボクセル上の脳機能情報の時系列データρk(l,m,i)を表している。この例では、ρk(l,m,i)=5,ρk(l-1,m,i)=ρk(l-1,m-1,i)=ρk(l,m-1,i)=ρk(l,m+1,i)=ρk(l+1,m,i)=2,ρk(l+1,m-1,i)=ρk(l-1,m+1,i)=ρk(l+1,m+1,i)=3であり、このうち、ρk(l,m,i)の平滑化に用いられるのは、ρk(l,m-1,i),ρk(l,m+1,i),ρk(l-1,m,i),ρk(l+1,m,i)の4つである。さらに、D’k(l,m,k)=0(つまり、Z方向へのプロトンの拡散がないということ)、且つ、D’l(l,m,k)=0.8,D’m(l,m,k)=0.2,w=0.4,W=0.6を満たすような状況を想定する。このとき、結合度ベクトル
【数85】
JP0004308871B2_000086t.gif

【0176】
は、式(31)から
【数86】
JP0004308871B2_000087t.gif

【0177】
となり、位置(l,m,k)のボクセル上の脳機能情報の時系列データρ(l,m,k,i)=5は、式(36)から、
【数87】
JP0004308871B2_000088t.gif

【0178】
のように平滑化されることがわかる。この処理を該2次元スライス画像Skの全ボクセル上の脳機能情報の時系列データρ(l,m,k,i)に施すことによって、神経繊維の走行方向に沿った方向(つまり、平均拡散度ベクトル
【数88】
JP0004308871B2_000089t.gif

【0179】
の成分D’l(l,m,k),D’m(l,m,k),D’k(l,m,k)のうち、最も大きい成分の方向:この例では、X方向)からの寄与が大きいように平滑化することができる。
【0180】
また、これ以外にも、式(5)の箇所で説明したように、最大固有値λMに対応する固有ベクトル
【数89】
JP0004308871B2_000090t.gif

【0181】
の向きと合致する2個のボクセルとのみ平滑化してもよい。
【0182】
次に、ステップS140において、データ解析手段26は、ステップS140において平滑化された脳機能情報の時系列データ
【数90】
JP0004308871B2_000091t.gif

【0183】
に対して、SPM等の解析手法を用いてデータ解析を施す。
【0184】
次に、ステップS150において、画像生成手段7は、データ解析手段26によって解析された結果をグレイスケールやフルカラー画像等の色調画像を生成する。
【0185】
そして、ステップS160において、画像表示手段8は、上記のようにして生成された画像を立体表示する。
【0186】
このように、第2の実施形態においては、MRI装置から取得した脳機能情報の時系列データを、同じくMRI装置から取得した拡散テンソルデータから算出された結合度ベクトルに基づいて平滑化処理したことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。
【0187】
なお、式(31)の結合度ベクトル
【数91】
JP0004308871B2_000092t.gif

【0188】
や式(32)の平滑化処理のアルゴリズムはあくまでも一例であって、その他の実施例も可能である。さらに、本実施形態におけるデータ平滑化処理の手法は、特定のデータ解析手法に依存しないことは上記の説明から明らかである。
【0189】
<第2の実施形態の変形例>
第2の実施形態においては、ステップS130において脳機能情報の時系列データρ(l,m,k,i)に対して拡散テンソルデータD(l,m,k)から算出された結合度ベクトル
【数92】
JP0004308871B2_000093t.gif

【0190】
を用いて平滑化処理を施した後に、ステップS140においてSPM等の解析手法を用いてデータ解析を行ったが、本発明の目的を達成するにあたって、これらのステップを逆にするような手法も考えられる。そこで、次に、上記した第2の実施形態の変形例について説明する。なお、重複を避けるため、以降においては第2の実施形態と異なる部分についてのみ説明を加える。
【0191】
図13は、第2の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。本変形例に係る脳機能解析装置20’の構成は、第2の実施形態の構成と基本的には同じであるが、データ平滑化手段200とデータ解析手段26とが入れ替わった構成になっている。これによって、脳機能解析装置20’によって行われるデータ解析の処理手順が第2の実施形態と異なる。なお、本実施形態においては、脳機能解析装置20’を、本実施形態の主要構成手段20A’と画像表示手段8及び記憶手段9とを一体型としたディスプレイコンソール形式としたが、画像表示手段8、或いは、画像表示手段8及び記憶手段9をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段20A’から切り離した構成としてもよい。
【0192】
図14は、図13に示した脳機能解析装置20’によって行われるデータ解析の処理手順の一例を示したフローチャートである。本変形例では、ステップS121において、まず、データ解析手段26’は、データ前処理手段3によって前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法を用いてデータ解析する。次に、ステップS131において、ボクセル間結合度算出手段24は、第2の実施形態と同様にして、式(31)で示したような結合度ベクトル
【数93】
JP0004308871B2_000094t.gif

【0193】
を算出する。次いで、ステップS141において、データ平滑化手段200’は、データ解析手段26’によって解析された結果である脳機能データρ(l,m,k)に対して、第2の実施形態と同様にして、式(32)を用いて平滑化処理を行う。次いで、ステップS150において、画像生成手段7は、上記のようにして得られた平滑化処理後のデータ
【数94】
JP0004308871B2_000095t.gif

【0194】
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。そして、ステップS160において、画像表示手段8は、上記のようにして生成された画像を立体表示する。
【0195】
このようにして、本変形例においても、第2の実施形態と同様の効果が得られる。
【0196】
[第3の実施形態]
図15は、本発明の第3の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第3の実施形態の脳機能解析装置30は、第1の実施形態の脳機能解析装置10におけるボクセル間結合度算出手段4、データ評価量構成手段5、及びデータ解析手段6の代わりに、ボクセル間結合度算出手段34、クラスタリング手段300、及びデータ解析手段36を設けた構成を成している。第1の実施形態と同一の構成には同一の符号を付し、以降においては、繰り返しの説明を省略する。なお、本実施形態においては、脳機能解析装置30を、本実施形態の主要構成手段30Aと画像表示手段8及び記憶手段9とを一体型としたディスプレイコンソール形式としたが、画像表示手段8、或いは、画像表示手段8及び記憶手段9をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段30Aから切り離した構成としてもよい。
【0197】
ボクセル間結合度算出手段34は、拡散テンソルデータ取得手段2で取得された拡散テンソルデータD(l,m,k)(式(9))から第1の実施形態で詳述したようにして平均拡散度ベクトル
【数95】
JP0004308871B2_000096t.gif

【0198】
(式(12))を算出した後、該平均拡散度ベクトル
【数96】
JP0004308871B2_000097t.gif

【0199】
の各成分D’l(l,m,k),D’m(l,m,k),D’k(l,m,k)の大きさに応じてボクセル同士をクラスタリングするための指標として結合度ベクトル
【数97】
JP0004308871B2_000098t.gif

【0200】
を算出する。
【0201】
クラスタリング手段300は、ボクセル間結合度算出手段34によって算出された結合度ベクトル
【数98】
JP0004308871B2_000099t.gif

【0202】
を用いてボクセルをクラスタリングする。
【0203】
データ解析手段36は、クラスタリング手段300でクラスタリングされたボクセル集団(クラスター)を単位として、データ前処理手段3で前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法を用いてデータ解析を行う。
【0204】
図16は、図15に示した第3の実施形態に係る脳機能解析装置30によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【0205】
ステップS200において、脳機能データ取得手段1は、MRI装置50で計測された脳機能情報の原時系列データρ’(l,m,k,i)を取得する。同ステップにおいて、拡散テンソル情報取得手段2は、同じくMRI装置50で計測された拡散テンソルデータD(l,m,k)を取得する。
【0206】
次に、ステップS210において、データ前処理手段3は、ステップS200において取得した脳機能情報の原時系列データρ’(l,m,k,i)に対して、第1の実施形態で詳述した前処理(a)~(d)を施すことによって、脳機能情報の時系列データρ(l,m,k,i)を生成する。
【0207】
次に、ステップS220において、ボクセル間結合度算出手段34は、まず、ステップS200において取得した拡散テンソルデータD(l,m,k)から第1の実施形態で詳述したようにして平均拡散度ベクトル
【数99】
JP0004308871B2_000100t.gif

【0208】
を算出する。次いで、該平均拡散度ベクトル
【数100】
JP0004308871B2_000101t.gif

【0209】
の各成分D’l(l,m,k),D’m(l,m,k),D’k(l,m,k)の大きさに応じてボクセル同士をクラスタリングする際の指標となる結合度ベクトル
【数101】
JP0004308871B2_000102t.gif

【0210】
を、例えば、
【数102】
JP0004308871B2_000103t.gif

【0211】
として算出する。ここで、
【数103】
JP0004308871B2_000104t.gif

【0212】
はクラスタリングの閾値を表す定数ベクトルであり、
【数104】
JP0004308871B2_000105t.gif

【0213】
と書ける。
【0214】
次に、ステップS230において、クラスタリング手段300は、ステップS220において算出された結合度ベクトル
【数105】
JP0004308871B2_000106t.gif

【0215】
の各成分Sl(l,m,k),Sm(l,m,k),Sk(l,m,k)の値が正になるボクセルを、それぞれX方向、Y方向、Z方向に対してクラスタリングする。図17は、図15に示したクラスタリング手段300によって行われるクラスタリングの手法を例示した図であり、ある2次元スライス画像Skの16個のボクセルを表している。図17(A)における隣接ボクセル間の矢印の値は、X方向又はY方向の平均拡散度ベクトルの値を表している。式(37)の定数ベクトル
【数106】
JP0004308871B2_000107t.gif

【0216】
の各成分の値をそれぞれ、例えば、hl=hm=3,hk=0とすると、結合度ベクトル
【数107】
JP0004308871B2_000108t.gif

【0217】
の成分Sl(l,m,k),Sm(l,m,k)の各値は、図17(B)のようになる。このとき、X方向、Y方向の各方向において、結合度ベクトル
【数108】
JP0004308871B2_000109t.gif

【0218】
の成分Sl(l,m,k),Sm(l,m,k)が共に正となる隣接ボクセルを順次結んでいくことによって、最終的に、点線で囲んだボクセルが一つのクラスターとして形成されることになる。結合度ベクトル
【数109】
JP0004308871B2_000110t.gif

【0219】
が平均拡散度ベクトルD’(l,m,k)の関数となっていることから、このようなクラスタリングは神経線維の走行方向を反映したものとなっている。
【0220】
次に、ステップS240において、データ解析手段36は、ステップS230においてクラスタリングされたボクセル集団(クラスター)毎に、ステップS210において前処理された脳機能情報の時系列データρ(l,m,k,i)を、SPM等の解析方法を用いてデータ解析する。その際には、例えば、一つのクラスター内での脳機能情報の時系列データを平均化する等の処理が必要となる。
【0221】
次に、ステップS250において、画像生成手段7は、ステップS240においてデータ解析された後のデータ
【数110】
JP0004308871B2_000111t.gif

【0222】
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。
【0223】
そして、ステップS260において、画像表示手段8は、ステップS250において生成された画像を立体表示する。
【0224】
このように、第3の実施形態においては、MRI装置から取得した拡散テンソルデータに基づいてボクセルをクラスタリングした後に、同じくMRI装置から取得した脳機能情報の時系列データをクラスター毎に解析したことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。
【0225】
なお、式(36)の結合度ベクトル
【数111】
JP0004308871B2_000112t.gif

【0226】
はあくまでも一例であってその他の実施例も可能である。また、式(5)の箇所で説明したように、最大固有値λMに対応する固有ベクトル
【数112】
JP0004308871B2_000113t.gif

【0227】
の向きと合致する2個のボクセルとのみクラスタリングしてもよい。さらに、本実施形態におけるクラスタリング処理の手法は、特定のデータ解析手法に依存しないことは上記の説明から明らかである。
【0228】
<第3の実施形態の変形例1>
第3の実施形態においては、ステップS230において、拡散テンソルデータD(l,m,k)から算出された結合度ベクトル
【数113】
JP0004308871B2_000114t.gif

【0229】
を用いてボクセルをクラスタリングした後に、ステップS240においてクラスター毎に脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法を用いてデータ解析したが、本発明の目的を達成するにあたって、これらのステップを逆にするような手法も考えられる。そこで、次に、上記した第3の実施形態の変形例について説明する。なお、重複を避けるため、以降においては第3の実施形態と異なる部分についてのみ説明を加える。
【0230】
図18は、第3の実施形態の変形例に係る脳機能解析装置の概略構成を示したブロック図である。本変形例に係る脳機能解析装置30’の構成は、第3の実施形態の構成と基本的には同じであるが、クラスタリング手段300とデータ解析手段36とが入れ替わった構成になっている。これによって、脳機能解析装置30’によって行われるデータ解析の処理手順が第3の実施形態と異なる。なお、本実施形態においては、脳機能解析装置30’を、本実施形態の主要構成手段30A’と画像表示手段8及び記憶手段9とを一体型としたディスプレイコンソール形式としたが、画像表示手段8、或いは、画像表示手段8及び記憶手段9をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段30A’から切り離した構成としてもよい。
【0231】
図19は、図18に示した脳機能解析装置30’によって行われるデータ解析の処理手順の一例を示したフローチャートである。本変形例では、ステップS221において、まず、データ解析手段36’は、データ前処理手段3によって前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法でデータ解析する。次に、ステップS231において、ボクセル間結合度算出手段34は、第3の実施形態と同様にして、式(36)で示したような結合度ベクトル
【数114】
JP0004308871B2_000115t.gif

【0232】
を算出する。次いで、ステップS241において、クラスタリング手段300’は、データ解析手段36’によって解析された脳機能データρ(l,m,k)に対して、第3の実施形態と同様にして、図17に示したようなボクセルのクラスタリングを行う。次いで、ステップS250において、画像生成手段7は、上記のようにしてクラスタリングされた脳機能データ
【数115】
JP0004308871B2_000116t.gif

【0233】
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。そして、ステップS260において、画像表示手段8は、上記のようにして生成された画像を立体表示する。
【0234】
このようにして、本変形例においても、第3の実施形態と同様の効果が得られる。
【0235】
<第3の実施形態の変形例2>
第3の実施形態の変形例1においては、データ前処理手段3によって前処理された脳機能情報の時系列データρ(l,m,k,i)をSPM等の解析手法でデータ解析した後に、ボクセル間結合度算出手段34で算出された結合度ベクトル
【数116】
JP0004308871B2_000117t.gif

【0236】
を用いてデータ解析後の脳機能データρ(l,m,k)をクラスタリングしたが、データ解析の手法として分類的解析手法を用いることも考えられる。そこで、本変形例では、その一例として、データ解析の手法として分類的解析手法の一つである決定木を用いた場合について説明する。例えば、図4に示したように、説明変数(属性)を脳機能情報の時系列データρ(l,m,k)、外的基準(クラス)をタスクとレストの時系列データとして決定木の手法を用いた場合、不要な属性(ボクセル)は除去され、必要な属性(ボクセル)のみが残る。このようにして残ったボクセルを結合度ベクトル
【数117】
JP0004308871B2_000118t.gif

【0237】
を用いてクラスタリングする。例えば、(1)あるボクセルを選び、(2)そのボクセルの周囲で結合度が一定値以上のボクセルをクラスタリングし、(3)(2)の処理を、周囲に結合度が一定値以上のボクセルがなくなるまで続け、(4)ボクセルが残っていれば、あるボクセルを選び、(2)の処理にもどる。
【0238】
次いで、ステップS250において、画像生成手段7は、上記のようにしてクラスタリングされた脳機能データ
【数118】
JP0004308871B2_000119t.gif

【0239】
に対して、所定の処理規則に従って、グレイスケールやフルカラー画像等の色調画像を生成する。そして、ステップS260において、画像表示手段8は、上記のようにして生成された画像を立体表示する。
【0240】
このようにして、本変形例においても、第3の実施形態と同様の効果が得られる。
【0241】
[第4の実施形態]
図20は、第4の実施形態に係る脳機能解析装置の概略構成を示したブロック図である。第4の実施形態の脳機能解析装置40は、第1の実施形態の脳機能解析装置10におけるボクセル間結合度算出手段4、データ評価量構成手段5、及びデータ解析手段6の代わりに、ボクセル間結合度算出手段44、データ検定手段400を設けた構成を成している。第1の実施形態と同一の構成には同一の符号を付し、以降においては、繰り返しの説明を省略する。なお、本実施形態においては、脳機能解析装置40を、本実施形態の主要構成手段40Aと画像表示手段8及び記憶手段9とを一体型としたディスプレイコンソール形式としたが、画像表示手段8、或いは、画像表示手段8及び記憶手段9をそれぞれ独立した画像表示装置及び記憶装置として主要構成手段40Aから切り離した構成としてもよい。
【0242】
ボクセル間結合度算出手段44は、データ前処理手段3で前処理された脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に各種データ検定(例えば、t検定、f検定、又はz検定)する際の基準値(検定値や有意水準など)を調整するための結合度ベクトル
【数119】
JP0004308871B2_000120t.gif

【0243】
を、式(12)の平均拡散度ベクトル
【数120】
JP0004308871B2_000121t.gif

【0244】
の関数として、拡散テンソルデータ取得手段2で取得された拡散テンソルデータD(l,m,k)(式(9))から算出する。
【0245】
データ検定手段400は、ボクセル間結合度算出手段44で算出された結合度ベクトル
【数121】
JP0004308871B2_000122t.gif

【0246】
によって該基準値を修正し、修正した基準値に基づき、前処理後の脳機能データρ(l,m,k)に対してボクセル毎に該データ検定を行う。
【0247】
図21は、図20に示した脳機能解析装置40によって行われるデータ解析の処理手順の一例を示したフローチャートである。
【0248】
ステップS300において、脳機能データ取得手段1は、MRI装置50で計測された脳機能情報の原時系列データρ’(l,m,k,i)を取得する。同ステップにおいて、拡散テンソル情報取得手段2は、同じくMRI装置50で計測された拡散テンソルデータD(l,m,k)を取得する。
【0249】
次に、ステップS310において、データ前処理手段3は、ステップS300において取得した脳機能情報の原時系列データρ’(l,m,k,i)に対して、第1の実施形態で詳述したような前処理(a)~(d)を施すことによって、脳機能情報の時系列データρ(l,m,k,i)を生成する。
【0250】
次に、ステップS320において、ボクセル間結合度算出手段44は、ステップS310において前処理された脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に、例えば、z検定する際の基準値(検定値や有意水準など)を調整するための結合度ベクトル
【数122】
JP0004308871B2_000123t.gif

【0251】
を、式(12)の平均拡散度ベクトル
【数123】
JP0004308871B2_000124t.gif

【0252】
の関数として、ステップS310で取得した拡散テンソルデータD(l,m,k)から算出する。
【0253】
次に、ステップS330において、データ検定手段400は、ステップS310において前処理された脳機能情報の時系列データρ(l,m,k,i)をボクセル毎にz検定を行いzスコアを算出する。そして、同ステップにおいて、算出したzスコアが大きなボクセル(例えば、ボクセル(l,m,k))を基点にして、その周辺(26方位)のボクセル(例えば、ボクセル(l+1,m,k))における有意水準の値a(l+1,m,k)を、ステップS320で算出した結合度ベクトル
【数124】
JP0004308871B2_000125t.gif

【0254】
に応じて、例えば、
【数125】
JP0004308871B2_000126t.gif

【0255】
のように下げる。ここで、αは重み係数である。
【0256】
次に、ステップS340において、画像生成手段7は、ステップS330において修正された有意水準に基づき、有意であると判断されたボクセル上の脳機能データρ(l,m,k)に基づく値を画像データとして生成する。
【0257】
そして、ステップS350において、画像表示手段8は、ステップS340において生成された画像データを立体表示する。
【0258】
なお、上記においては、拡散テンソルデータ取得手段2で取得した拡散テンソルデータD(l,m,k)から算出した式(38)の結合度ベクトル
【数126】
JP0004308871B2_000127t.gif

【0259】
を用いて検定の有意水準の値を調整したが、その代わりに、検定値を調整してもよい。
【0260】
例えば、タスク画像とレスト画像との間でボクセル毎にz検定を行う場合を考えてみる。タスク画像とは、あるタスクを行っているときの画像(例えば、図2の太線で表された1~3のスライス画像)のことをいい、レスト画像とは、該タスクを行っていないときの画像(例えば、図2の太線で表された5~7のスライス画像)のことをいう。このとき、ボクセル(l,m,k)におけるzスコア(z検定の検定値)は、通常、
【数127】
JP0004308871B2_000128t.gif

【0261】
と表される。ここで、Mt(l,m,k)は、タスク画像の脳機能情報の時系列データρ(l,m,k,i)の平均値、Mc(l,m,k)は、レスト画像の脳機能情報の時系列データρ(l,m,k,i)の平均値、σt(l,m,k)はタスク画像の脳機能情報の時系列データρ(l,m,k,i)の標準偏差、σc(l,m,k)は、レスト画像の脳機能情報の時系列データρ(l,m,k,i)の標準偏差を表す。この前提の下、z=0ならば、タスク画像における脳機能情報の時系列データρ(l,m,k,i)の平均値とレスト画像の脳機能情報の時系列データρ(l,m,k,i)の平均値との間に差がないといえる。Z=1ならば、タスク画像における脳機能情報の時系列データρ(l,m,k,i)の平均値とレスト画像の脳機能情報の時系列データρ(l,m,k,i)の平均値との差が時系列データのばらつきの分だけ大きいことを意味する。このような傾向から、zスコアの値が大きいボクセルがタスクに関連している活性部位であると考えることができる。
【0262】
そこで、式(40)を
【数128】
JP0004308871B2_000129t.gif

【0263】
のように変更して、ステップS330において、b(l,m,k)を、ステップS320で算出した式(38)の結合度ベクトル
【数129】
JP0004308871B2_000130t.gif

【0264】
を用いて、例えば、
【数130】
JP0004308871B2_000131t.gif

【0265】
と修正する。ここで、βは重み係数である。
【0266】
次いで、ステップS340において、画像生成手段7は、全ボクセルにおいて共通の有意水準に基づき、有意であると判断されたボクセル上の脳機能データρ(l,m,k)に基づく値を画像データとして生成する。そして、ステップS350において、画像表示手段8は、ステップS340において生成された画像データを立体表示する。
【0267】
このように、第4の実施形態においては、MRI装置から取得した拡散テンソルデータに基づいて、同じくMRI装置から取得した脳機能情報の時系列データを検定する際の基準値(検定値又は有意水準)の値をボクセル毎に調整したことによって、脳の活性部位間の結合構造をも考慮した脳の活性部位の特定が可能となる。
【0268】
<第4の実施形態の変形例>
第4の実施形態では、脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に各種データ検定(例えば、T検定、F検定、又はZ検定)する際に、拡散テンソルデータ取得手段2で取得した拡散テンソルデータD(l,m,k)から算出した結合度ベクトル
【数131】
JP0004308871B2_000132t.gif

【0269】
を用いて該データ検定の基準値(検定値や有意水準など)を調整したわけだが、本発明の目的を達成するにあたって、脳機能情報の時系列データρ(l,m,k,i)をボクセル毎に各種の解析手法でデータ解析したときに、該解析手法が統計的に有意であるか否かをデータ検定(例えば、T検定、F検定、又はZ検定)するような状況も考えられる。そこで、本変形例として、その際の基準値(検定値や有意水準など)を、結合度ベクトル
【数132】
JP0004308871B2_000133t.gif

【0270】
を用いて調整することを考える。なお、重複を避けるため、以降においては第4の実施形態と異なる部分についてのみ説明を加える。
【0271】
図22は、第4の実施形態の変形例に係る脳機能解析装置の概略構成を示した概念図である。本変形例に係る脳機能解析装置40’の構成は、第4の実施形態の構成のデータ検定手段400の代わりに、データ検定手段400’とデータ解析手段46’とを加えた構成を成す。
【0272】
データ解析手段46’は、脳機能データ取得手段1で取得され、データ前処理手段3で前処理された脳機能情報の時系列データρ(l,m,k,i)を、従来の分類、回帰、相関等の解析手法を用いてボクセル毎にデータ解析する。
【0273】
データ検定手段400’は、データ解析手段46’で行われた解析手法が統計的に有意であるか否かを各種データ検定(例えば、T検定、F検定、又はZ検定)に基づき検定する際の基準値(検定値や有意水準など)を、結合度ベクトル
【数133】
JP0004308871B2_000134t.gif

【0274】
を用いて調整する。
【0275】
図23は、図22に示した脳機能解析装置によって行われるデータ解析の処理手順の一例を示したフローチャ-トである。
【0276】
ステップS360において、データ解析手段46’は、脳機能データ取得手段1で取得され、データ前処理手段3で前処理された脳機能情報の時系列データρ(l,m,k,i)を、従来の分類、回帰、相関等の解析手法を用いてボクセル毎にデータ解析する。
【0277】
1.分類的解析手法(判別分析、決定木、サポートベクターマシン等)
1-1:判別分析の場合
説明変数を脳機能情報の時系列データρ(l,m,k,i)、外的基準をタスクとレストの時系列データとして、各種の判別関数を用いてデータ分類する。
【0278】
1-2:決定木の場合
説明変数(属性)をタスクとレストの時系列データ、外的基準として脳機能情報の時系列データρ(l,m,k,i)を離散化して幾つかのクラスを作り、決定木を用いてデータ分類する。
【0279】
2.回帰的解析手法(各種の関数(線形、非線形)を用いた回帰分析)
2-1:線形回帰分析の場合
説明変数をタスクとレスト、被説明変数を脳機能情報の時系列データρ(l,m,k,i)として、線形回帰式を用いてデータ解析する。
【0280】
3.相関的解析手法(推定モデル(関数)と実際のデータとの相関を取る手法)
予め設定した、モデル(例えば、刺激を入力とし、脳機能情報の時系列データを出力とするような関数)と実際の脳機能情報の時系列データρ(l,m,k,i)との相関を取ることによって、データ解析する。
【0281】
ついで、ステップS320(第4の実施形態と同様)を経て、ステップS331において、データ検定手段400’は、ステップS360において得られた結果が統計的に有意であるか否かを検定する。このとき、検定の基準値(検定値や有意水準など)を結合度ベクトル
【数134】
JP0004308871B2_000135t.gif

【0282】
(式(38))を用いて調整する。例えば、分類的解析手法で得られた結果に対してf検定を行う場合には、ステップS360において得られた結果に対してボクセル毎にf検定を行いf値を算出する。そして、同ステップにおいて、算出したf値が大きなボクセル(例えば、ボクセル(l,m,k))を基点にして、その周辺(26方位)のボクセル(例えば、ボクセル(l+1,m,k))における有意水準の値a(l+1,m,k)を、ステップS320で算出した結合度ベクトル
【数135】
JP0004308871B2_000136t.gif

【0283】
に応じて、例えば、式(39)を用いて調整する。回帰的解析手法の場合には、閾値(有意水準に相当)の値を上記のようにして調整すればよいし、相関的解析手法の場合には、相関の閾値を上記のようにして調整すればよい。
【0284】
次に、ステップS340において、画像生成手段7は、ステップS330において修正された有意水準に基づき、有意であると判断されたボクセル上の脳機能データρ(l,m,k)に基づく値を画像データとして生成する。
【0285】
そして、ステップS350において、画像表示手段8は、ステップS340において生成された画像データを立体表示する。
【0286】
このようにして、本変形例においても、第3の実施形態と同様の効果が得られる。
【0287】
以上で説明した実施形態はあくまでも例示的なものであり、これに限定されるものではない。
【0288】
例えば、各実施形態において、隣接ボクセル間の結合度を算出するステップを脳機能データ及び拡散テンソルデータを取得するステップのすぐ後に持ってきてもよい。
【0289】
また、データ解析の手法についても、ウェブレット解析、ロジスティック回帰分析解析等を適用することも可能である。
【0290】
また、第1の実施形態に第2又は第3の実施形態を組み込むことも可能である。まず、第1の実施形態に第2の実施形態を組み込む場合には、図1に示した脳機能解析装置10に、図10又は図13に示したデータ平滑化手段200又は200’と同じ機能を有するデータ平滑化手段を更に備えるような構成とすればよい。そして、第1の実施形態のステップS50において、データ解析手段6がデータ評価量Qを最小にするような回帰係数を求めた後に、該回帰係数の値に対して上記データ平滑化手段を適用することにより該回帰係数の値を異方的に平滑化すればよい。これにより、第1の実施形態で得られる画像と比較した場合、白質(神経繊維)の活性度がより強調された画像を生成することが可能となる。
【0291】
また、第1の実施形態に第3の実施形態を組み込む場合には、図1に示した脳機能解析装置10に、図15又は図18に示したクラスタリング手段300又は300‘と同じ機能を有するクラスタリング手段を更に備えるような構成とすればよい。そして、第1の実施形態のステップS50において、データ解析手段6がデータ評価量Qを最小にするような回帰係数を求めた後に、該回帰係数の値に対して上記クラスタリング手段を適用することによりボクセル同士をクラスタリングさせればよい。これにより、第1の実施形態で得られる画像と比較した場合、白質(神経線維)の活性度がより強調された画像を生成することが可能となる。
【0292】
さらに、上記した各実施形態で得られた結果に対して、さらに、拡散テンソルデータD(l,m,k)を用いて、白質ボクセルの活性化部位を伸ばして行くことも可能である。その手法は基本的にトラクトグラフィー(tractography)と同様の手法である。
【0293】
さらに、上記の手法(回帰分析、相関、分類、独立成分分析、検定等)を白質ボクセルだけに適用することで、白質の活性化を捕らえることも可能である。
【0294】
さらに、各実施形態における実験としてはブロックデザインを想定したが、本発明は、事象関連fMRIにも適用可能である。
【0295】
また、本発明は接続解析(Connectivity Analysis)にも適用可能である。接続解析の手法にはいくつかあるが、たとえば、SEM(Structural Equation Modeling)(あるいは共分散構造分析)の場合には、その計算過程で相関や検定が行われるが、その際に本発明を適用できるし、平滑化やクラスタリングも適用できる。
【0296】
上記では、神経細胞と神経線維の結合に関する情報を拡散テンソルデータから獲得したが、すでに解剖学的に既知となっている結合に関する情報(知見)があれば、それを用いることも可能である。
【0297】
本発明の趣旨を逸脱しない限りにおいて、その他にも様々な実施形態が可能であり、本発明はあくまでも特許請求の範囲によってのみ規定されるものである。
【0298】
最後に、本発明の脳機能解析方法を用いた脳機能解析の効果を示す。図24は、被験者に聴覚刺激を用い単純復唱を行わせた場合の各種脳機能解析の結果を示したものである。実験は、ヘッドホンから流れてくる音声を聞こえた通りに復唱することをタスクとし、レストは出来るだけ何も考えない状態とした。タスクのボリューム数(計測回数)は48,レストのボリューム数(計測回数)は60,解析に使用しないボリューム数(計測回数)は12とした(図2を参照)。また、TR(Time of Repetition:図2のボリューム幅に対応)は5000msとした。被験者の数は8人であるが、SPMと本発明の脳機能解析方法の両方で良好な結果を得た被験者の解析結果を以下に示す。
【0299】
図24A,24Bは共にSPMを用いたT検定の結果を表すものであり、図24AはT検定の閾値に修正を施した場合の解析結果を表し、図24BはT検定の閾値に修正を施さなかった場合の解析結果を表している。図24Cは、SPMとトラクトグラフィーとを組み合わせた解析手法による解析結果を表し、図24Dは本発明の第1の実施形態に係る脳機能解析方法を用いた解析結果を表す。なお、SPMの前処理として、位置補正、正規化、平滑化を行った。平滑化の半値幅は9mmとした。また、トラクトグラフィーには拡散テンソル解析用のソフトウェアとしてdTVを用いた。
【0300】
前提として、単純復唱の時には、ウェルニッケ野(感覚性言語野)とブローカー野(運動性言語野)が活性化し、両野は弓状束という神経繊維の束で結合されていることが解剖学的に知られている。この弓状束は発話の際に用いられている。
【0301】
図24Aに示すように、SPMを用いた場合でも、T検定の閾値を修正した場合にはウェルニッケ野、ブローカー野の活性化は検出できない。一方、図24Bに示すように、SPMを用い、T検定の閾値を修正しなかった場合には、両野の活性化は検出できたが、両野をつなぐ弓状束は検出できなかった。これらの結果に対し、図24Cに示すように、SPMとトラクトグラフィーとを組み合わせた解析手法を用いると、ウェルニッケ野及びブローカー野の活性化の検出に加え、弓状束を検出できた。
【0302】
一方、本発明の第1の実施形態に係る脳機能解析方法を用いた場合には、図24Dに示すように、両野の活性化や両野をつなぐ弓状束の検出に加えて、弓状束の活性化(弓状束を構成する神経線維におけるパルス伝達)まで検出できる。
【0303】
また、SPMとトラクトグラフィーとを組み合わせた解析手法においては、神経繊維(この場合、弓状束)の検出にあたり、解析者が神経繊維の追跡開始点と終了点とを予め指定しなければならない。したがって、神経線維の追跡開始点と終了点とが予め決定できないような場合(例えば、灰白質の活性領域が予めわかっていないようなタスク)にはこの手法は有効でないことになる。それに対し、本発明の第1の実施形態に係る脳機能解析方法によれば、神経線維の追跡開始点と終了点とが予め指定できない場合に対しても、図24Dに示したように、神経線維の活性化の検出が可能である。
【産業上の利用可能性】
【0304】
このように、本発明によれば、脳の活性部位を特定できるのはもちろんのこと、解析者が予め神経線維の開始点及び終了点を指定しなくとも、神経線維(活性部位間の結合構造)を特定することが可能となり、さらには神経繊維の活性化(神経線維におけるパルス伝達)をも検出可能となるので、認知症、失語症、統合失調症等の高次脳機能障害の診断及び治療に対して極めて有効な脳機能解析方法、及び脳機能解析プログラムを提供することができる。
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14
【図16】
15
【図17】
16
【図18】
17
【図19】
18
【図20】
19
【図21】
20
【図22】
21
【図23】
22
【図24】
23