JP2009537201A - 灌流イメージングの方法 - Google Patents

Abstract

本発明は、灌流イメージングの方法、ならびに後イメージングデータ操作の方法およびそのためのプログラムに関する。

Description

本発明は、灌流イメージングの方法、ならびに後イメージングデータ操作およびそのためのプログラムに関する。

灌流イメージングは、体内の物理構造ならび組織の機能・生存度の両方を検出するために使用され得る。灌流イメージングは、脳、心臓または肝臓に損傷（発作、腫瘍、梗塞などによって生じる）のある患者を調べるために、特に有用である。灌流イメージング手法の例としては、磁気共鳴イメージング（ＭＲＩ）、医用超音波検査法（音波造影法）、陽電子断層撮影法（ＰＥＴ）およびコンピュータ断層撮影法（ＣＴ）がある。

一般のＭＲＩ灌流手法において、造影剤（例えば、アマルシャム（Ａｍｅｒｓｈａｍ）およびシェーリング（Ｓｃｈｅｒｉｎｇ）によってＯｍｎｉｓｃａｎ（登録商標）またはＭａｇｎｅｖｉｓｔ（登録商標）として市販されるものなどのガドリニウムキレート）が患者の血管系へボーラス投与され、そして造影剤ボーラスが対象領域における組織を通過する時間を含む期間にわたり対象領域からの画像が集められる。例えば、ＭＲＩにおいては、高速画像取り込みシーケンス（例えば、スピンエコー、グラジエントリコール（ＧＲＡＳＳまたはＦＬＡＳＨ）、エコープラナー（ＥＰＩ）、ＲＡＲＥ、ハイブリッド、半励起など）が使用される。灌流イメージングのためのＭＲ造影剤のそのようなシーケンスおよびボーラス投与は、当該分野において周知である（例えば、「Ｂｉｏｍｅｄｉｃａｌ Ｍａｇｎｅｔｉｃ Ｒｅｓｏｎａｎｃｅ Ｉｍａｇｉｎｇ」、Ｅｄ． Ｗｅｈｒｌｉ ｅｔ ａｌ、ＶＣＨ、１９８８を参照）。

臨床の現場においては、一連の灌流画像を精査し、そしてその結果を定性的に評価することが一般に行われる。

しかし、多くの場合、定量的な結果（例えば、局部的な血流量、局部的な血液体積、局部的な平均通過時間、局部的な到達時間、局部的な透過性表面積の積およびピークまでの局部的な時間の絶対測定値）が所望される（例えば、Ｒｅｍｐｐ ｅｔ ａｌ、Ｒａｄｉｏｌｏｇｙ １９３：６３７〜６４１（１９９４）およびＶｏｎｋｅｎ ｅｔ ａｌ、ＭＲＭ ４１：３４３〜３５０（１９９９）を参照）。

薬物動態モデリングを使用して、血漿流量、血漿体積、平均通過時間、血管外細胞外体積、透過性表面積の積およびピークまでの時間、またはより生理学的でない灌流パラメータ（毛細管移動定数Ｋ_{ｔｒａｎｓ}など）に対するボクセル別の値が抽出される。造影剤の再循環および造影剤の漏れが含まれる。

ＭＲマルチパス灌流イメージングのためにこれまでに導入された薬物動態モデルのうちの最も一般的なものは、ジョンソン（Ｊｏｈｎｓｏｎ）およびウィルソン（Ｗｉｌｓｏｎ）の断熱近似モデルである（ａａＪＷ）［Ｈｅｎｄｅｒｓｏｎ ｅｔ ａｌ． ＪＭＲＩ、Ｖｏｌ．１２：９９１〜１００３（２０００）］。しかし、ａａＪＷモデルには複数のパラメータが存在するので、ａａＪＷモデルは、比較不安定であり、かつ元のデータにおける低い信号対ノイズ比の影響を非常に受けやすい。

各ボクセルについて観測されるトレーサ灌流信号は、未知のボクセル別動脈入力関数および未知のボクセル別組織残留関数（または、インパルス応答関数）の畳み込みである。

組織残留関数は、時刻ｔにある組織領域内に造影剤のうちのどれだけがまだ存在するかを記述し、かつしたがってその組織の生理学的パラメータ（例えば、血液体積および平均通過時間）に依存する関数である。動脈入力関数は、どのように造影剤が組織ボクセルに送達されるかを記述し、かつしたがって器官内の血管構造の印象を与える。

薬物動態モデリングでは、未知のボクセル別組織残留関数は、既知のパラメータ形式を有すると仮定されるが、パラメータ値は未知である。トレーサ到達の時間遅延およびトレーサのピークまでの時間は、ボクセル別動脈入力関数因子のみによって指定される。すべての残りの灌流パラメータは、流量によって拡大・縮小される（ｆｌｏｗ ｓｃａｌｅｄ）ボクセル別組織残留関数因子によって指定される。

マルチパス薬物動態モデリングを使用して灌流パラメータを見つけるすべての公知の方法には問題がある。例えば、未知のボクセル別動脈入力関数は、１つの既知の動脈入力関数に置き換えられ得る。この値は、一般に、血管内のボクセルを手動または自動で特定することによって見つけられる。この血管は、対象組織から離れた位置に存在し得る。さらに、使用された血管から対象組織の動脈入力関数のすべての遅延およびばらつきは、無視される。これらのモデリングの単純化は、推定された灌流パラメータに大きな誤差を生じさせ得る。理想的には、ボクセル別動脈入力関数は、遅延およびばらつきを回避するように使用されるべきである。加えて、同じ患者を連続して検査した場合に測定される流速、移動および飽和が大きいと、大きな血管の動脈入力関数に大きなばらつきを招く。

驚くべきことに、本発明者らは、ボクセル別ａａＪＷモデルとブラインド逆畳み込み方法とを組み合わせることによって、特に磁気共鳴イメージング分野において、増強された灌流イメージングが可能であることを見出した。また、本発明の方法を使用することによって、動脈入力関数および組織残留関数の値が著しくより臨床的に有益となる。

また、本発明の方法は、血漿流量、血漿体積、平均通過時間、摂取率、透過性表面積の積、血管外細胞外液体積、およびピークまでの時間の高品質生理学的パラメトリック画像を提供する。これらのパラメトリック画像は、生理学的パラメータの直接的な定量的推定値を与える、すなわち、得られる結果は、生理学的に正確である。そのような高品質画像は、従来生成されなかったし、かつこれらのパラメトリック値のすべてを導出することが可能な公知の単一方法はなかった。さらに、本発明の方法は、実験ノイズに著しくは影響されない。

推定されたボクセル別動脈入力関数は、造影剤ボーラスが繰り返し画像化領域を通過する際に、造影剤ボーラスのばらつきおよび遅延における変化を指定する。組織残留関数は、手動で特定された動脈入力関数のみを使用することによって導入された遅延およびばらつきの効果を劣化することなく見つけられる。結果として、すべての血管パラメータに対する推定方式は、線形、簡単かつノイズに対して頑強となる。

したがって、一局面によると、本発明は、灌流イメージングデータを処理する方法であって、ブラインド逆畳み込みアルゴリズムならびにジョンソン（Ｊｏｈｎｓｏｎ）およびウィルソン（Ｗｉｌｓｏｎ）の断熱近似モデル（ａａＪＷ）の使用を含む、方法を提供する。好ましくは、データは、ボクセル毎に処理される、すなわち、ボクセル別ブラインド逆畳み込みアルゴリズムおよびボクセル別ａａＪＷモデルが使用される。好ましくは、処理対象のデータは、一続きの時間値［ｔ］において記録される。最も好ましくは、ボクセル別ブラインド逆畳み込みアルゴリズムおよびボクセル別ａａＪＷモデルは、一続きの時間値［ｔ］において処理された灌流イメージングデータに適用される。

本発明のさらなる局面は、灌流イメージングの方法であって、ブラインド逆畳み込みアルゴリズムおよびａａＪＷモデルの使用を介して、一続きの時間値［ｔ］において、被験者の対象領域の画像を対象領域のボクセルから生成することを含む、方法である。

本発明のさらなる局面は、灌流画像の画像処理の方法であって、ブラインド逆畳み込みアルゴリズムおよびａａＪＷモデルの使用および画像の生成を含む、方法である。好ましくは、データは、ボクセル毎に処理される、すなわち、ボクセル別ブラインド逆畳み込みアルゴリズムおよびボクセル別ａａＪＷモデルが使用される。好ましくは、処理対象の灌流画像は、一続きの時間値［ｔ］において記録された一続きのデータポイントを含む。

本発明の一好適な実施形態において、ブラインド逆畳み込みアルゴリズムおよびａａＪＷモデルは、パラメータであるボクセル別動脈入力関数ｃ_ｐ［ｔ］およびボクセル別組織残留関数ｒ［ｔ］の生成において使用される。

本発明の一好適な実施形態において、また、ボクセル別生理学的パラメータである血漿流量、血漿体積、平均通過時間、摂取率、透過性表面積の積、血管外細胞外液体積およびピークまでの時間のうちの少なくとも１つが決定される。

本発明の一好適な実施形態において、上記パラメータのうちの少なくとも１つまたはすべての画像が別々に生成される。

本発明の一好適な実施形態において、ブラインド逆畳み込み方法が繰り返し使用される。特に好ましくは、ブラインド逆畳み込み方法を繰り返し使用して、ボクセル別動脈入力関数ｃ_ｐ［ｔ］および拡大・縮小されるボクセル別組織残留関数ｕ［ｔ］の値を得る。

本発明の一好適な実施形態において繰り返しループは、４回以上の繰り返し回数を含み、さらに好ましくは、６回以上の繰り返し回数を含む。

本発明の一好適な実施形態において、ａａＪＷを使用して、拡大・縮小されるボクセル別組織残留関数ｕ［ｔ］からボクセル別組織残留関数ｒ［ｔ］を生成する。

本発明の一好適な実施形態において、ボクセル別動脈入力関数ｃ_ｐ［ｔ］の初期推定値は、観測されるボクセル別トレーサ濃度ｃ［ｔ］である。

本発明の一好適な実施形態において、ボクセル別動脈入力関数ｃ_ｐ［ｔ］の初期推定値は、第１のパス動脈入力関数である。

本発明の一好適な実施形態において、拡大・縮小されるボクセル別組織残留関数ｕ［ｔ］の積分値の初期推定値は、１である。

本発明の一好適な実施形態において、灌流イメージング方法は、磁気共鳴イメージング（ＭＲＩ）、医用超音波検査法（音波造影法）、陽電子断層撮影法（ＰＥＴ）、およびコンピュータ断層撮影法（ＣＴ）である。最も好ましくは、磁気共鳴イメージング（ＭＲＩ）である。

本発明の一好適な実施形態において、ブラインド逆畳み込み方法は、ルーシー−リチャードソン（Ｌｕｃｙ−Ｒｉｃｈａｒｄｓｏｎ）逆畳み込み方法、またはランドウェバー（Ｌａｎｄｗｅｂｅｒ）方法とウィーナーフィルタリングとの組み合わせである。最も好ましくは、ルーシー−リチャードソン逆畳み込み方法である。

本発明の一好適な実施形態において、低速水交換方式モデルがａａＪＷモデルに組み込まれる。

本発明の一好適な実施形態において、平均細胞内滞留時間のボクセル別の値が生成される。

本発明のさらなる局面は、本発明の方法において使用されるデータ処理装置であって、装置は、一続きの時間値（ｔ）に対して対象領域のボクセルについての１セットの灌流イメージング信号値を入力として受け取る手段と、ブラインド逆畳み込みアルゴリズムならびにジョンソンおよびウィルソンの断熱近似モデル（ａａＪＷ）の使用を含む、１セットの灌流イメージング信号値からパラメータ（好ましくは、薬物動態パラメータ）を生成する手段とを有する、装置である。

本発明のさらなる局面は、命令を搭載するコンピュータソフトウェア製品であって、命令は、本発明の方法において使用されるデータ処理装置上で実行される場合に、装置が、一続きの時間値（ｔ）に対して対象領域のボクセルについての１セットの灌流イメージング信号値を入力として受け取り、かつブラインド逆畳み込みアルゴリズムならびにジョンソンおよびウィルソンの断熱近似モデル（ａａＪＷ）の使用を含む、１セットの灌流イメージング信号値からパラメータの値を生成する手段を有するように構成されるようにする、コンピュータソフトウェア製品である。

本発明のさらなる局面においては、灌流イメージングデータを処理して、前方駆出率を生成する方法であって、薬物動態コンセンサスモデルおよびブラインド逆畳み込みアルゴリズムの使用を含む、方法が提供される。好ましくは、左心室の前方駆出率が生成される。
発明の詳細な説明
ａａＪＷの２コンパートメントモデルは、毛細管通過時間および血管内トレーサの組織増強に対する寄与の両方を考慮する。このモデルにおいて、トレーサが分配される空間は、２つの円柱に分割され、トレーサはその間で交換される。

内側の円柱は、血漿空間である。造影剤濃度は、時間およびこの円柱に沿った距離の両方とともに変化する。外側の円柱は、血管外細胞外空間であり、かつ十分に混合されたコンパートメントであると仮定される。この円柱内の造影剤濃度は、時間のみとともに変化する。造影剤は、２つのコンパートメントの間を流れ得るが、血漿空間と血管外細胞外空間との間の造影剤の交換に対する２つの速度定数は、同じである必要はない。

測定応答の線形性、および血漿流量が測定期間にわたって一定のままであることを仮定すると、ある体積の組織におけるトレーサの濃度は、以下のとおりである。

ここで、記号＊は、１次元の畳み込み演算子であり、Ｆは、血漿流量であり、およびｔ＝１、２、．．．、Ｔは、離散時間サンプルである。ｃ［ｔ］は、組織トレーサ濃度であり、ｃ_ｐ［ｔ］は、動脈血漿のトレーサ濃度であり、およびｒ［ｔ］は、組織残留関数である。

ａａＪＷモデルに対する組織残留関数（インパルス応答関数）は、以下のとおりである。

ここで、Ｔ_ｃは、平均毛細管通過時間であり、およびｋ_２は、血管外細胞外空間から血管内空間へのトレーサの通過に対する速度定数であり、ν_ｅは、組織（漏れ空間）の単位体積当たりの血管外細胞外体積であり、およびＥは、摂取率（ｅｘｔｒａｃｔｉｏｎ ｆｒａｃｔｉｏｎ）である。中央体積定理（ｃｅｎｔｒａｌ ｖｏｌｕｍｅ ｔｈｅｏｒｅｍ）を使用すると、ν_ｐ＝Ｆ・Ｔ_ｃであり、ここで、ν_ｐは、組織の単位体積当たりの血漿体積である。

平均移送時間に関するいくつかの仮定がａａＪＷモデルに適用され得る。平均毛細管通過時間Ｔ_ｃは、信号のサンプリング時間ｓ［ｔ］未満であり得る。そして、ａａＪＷモデルは、以下のとおりである。

Ｔｃ＝Δｓ、ここで、Δｓは、サンプリング間隔である（単位：秒）。再度中央体積定理を使用すると、この場合、ν_ｐ＝Ｆ・Δｓである。しかし、Ｔ_ｃの組織依存性分散についてのすべての情報は、失われる。

血管内血漿体積をすべてまとめて無視すると（すなわち、Ｔ_ｃ＝０）、以下が得られる。

しかし、ν_ｐおよびＴ_ｃについてのすべての情報が失われる。血漿流量または透過性表面積の積を推定するために、灌流パターン（流量制限型、拡散制限型）についてのさらなる仮定が必要とされる。

本発明の一好適な実施形態において、毛細管移動定数Ｋ_{ｔｒａｎｓ}＝Ｆ・Ｅは、ｋ_２と同じであると仮定され、ｔ≧Ｔ_ｃに対するａａＪＷモデル残留関数は、以下のとおりである。

摂取率Ｅは、毛細管床を通る１回のトレーサ通過の間に血漿から周辺の組織へ全方向に拡散するトレーサの割合である。摂取率Ｅは、Ｆおよび透過性表面積の積ＰＳに以下のとおり関係付けられる。

すべてのこれらの関数およびパラメータは、ボクセル毎に推定される。ボクセル別表記を導入すると、式（１）は、以下のとおり記載され得る。

ここで、ｉ＝１、２、．．．、Ｉおよびｊ＝１、２、．．．、Ｊである。ＩおよびＪは、それぞれｘおよびｙ方向のボクセルの数である。この形式の式（７）を後述のように使用して、２つの因子ｃ^ｉ，ｊ _ｐ［ｔ］およびｕ^ｉ，ｊ［ｔ］を繰り返しブラインド逆畳み込み手順によって別々に見つける。このように、各ボクセル（ｉ，ｊ）に対して、動脈血漿のトレーサ濃度関数ｃ^ｉ，ｊ _ｐ［ｔ］、およびそれに対応する拡大・縮小される組織残留関数ｕ^ｉ，ｊ［ｔ］が見つけられる。

ｃ^ｉ，ｊ _ｐ［ｔ］が対象のボクセルにおける組織トレーサ濃度関数から得られるので、トレーサボーラスが血管ツリーを通って対象のボクセルへ移動することによって生じるすべての遅延およびばらつきは、ｃ^ｉ，ｊ _ｐ［ｔ］に含まれる。これは、完全に異なる解剖学的位置からのボクセル別でないｃ_ｐ［ｔ］を使用することと対照的である。

任意の動脈内に完全に位置する任意のボクセルを通過するｃ^ｉ，ｊ _ｐ［ｔ］＝１、２、．．．、Ｔの積分値Ω^ｉ，ｊは、第１パスモデルにおけるすべてのボクセルに対して同じでなければならない。第１パスモデルは、以下のとおりである。

ここで、ｉ＝１、２、．．．、Ｉ、ｊ＝１、２、．．．、Ｊである。式（７）から分かるように、Ωは、Ｆ^ｉ，ｊおよびその他の血管体積灌流パラメータを絶対的に計算するためには既知でなければならない。反対に、時間依存性パラメータである平均通過時間Ｔｃ^ｉ，ｊおよびピークまでの時間τ^ｉ，ｊは、Ωを知らずに、逆畳み込みの前に、常に絶対値で与えられる。本発明の一好適な実施形態において、積分値Ωは、例えば、動脈、別の大きな血管内または左側の心室内に対象領域マスクを描くことによって手動で見つけられ得る。よりよい推定は、前ボーラス注入を使用することによって見つけられ得る。あるいは、積分値Ωは、自動的に見つけられ得る、例えば、集団統計を使用して、適用したトレーサ用量、体重および体の体積から推定され得る。

本発明の一実施形態において、ｃ^ｉ，ｊ _ｐ［ｔ］は、Ω＝Ω^ｉ，ｊ＝１、ｉ＝１、２、．．．、Ｉ、ｊ＝１、２、．．．、Ｊを設定することによって、対象のボクセルを通る全トレーサ量に対して、完全に自動的に拡大・縮小され得る。このようにスケーリングを制限することで、Ｆ^ｉ，ｊ、ｖ_ｐ ^ｉ，ｊおよびｖ_ｅ ^ｉ，ｊの推定値は、絶対値ではなく、相対値に限定される。

ｃ^ｉ，ｊ _ｐ［ｔ］およびＦ^ｉ，ｊ・ｒ^ｉ，ｊ［ｔ］が別々に見つけられるので、薬物動態パラメータの推定は、簡単になる。ｒ^ｉ，ｊ［ｔ］に遅延はなく、遅延は、ｃ^ｉ，ｊ _ｐ［ｔ］に存在する。したがって、ボクセルに対するピークまでの時間τ^ｉ，ｊは、単に、すべてのボクセルに共通の固定開始点からｃ^ｉ，ｊ _ｐ［ｔ］の最大点までの時間である。

本発明の一実施形態において、同時刻におけるすべてのボクセルに対するｃ^ｉ，ｊ _ｐ［ｔ］の時間展開は、従来の方法で表示され得る（例えば、ビデオを介して）。これにより、トレーサの遅延およびばらつきにおける局部的な差を視覚的に評価することが可能になる。

パラメータＴ_ｃ ^ｉ，ｊ、Ｆ^ｉ，ｊおよびＫ_{ｔｒａｎｓ} ^ｉ，ｊは、ｌｏｇ（Ｆ^ｉ，ｊ・ｒ^ｉ，ｊ［ｔ］）に対して二乗平均フィッティング手順を行うことによって同時に見つけられ得る。対数関数ｌｏｇは、Ｆ^ｉ，ｊ・ｒ^ｉ，ｊ［ｔ］の最後の因子を予想される減衰指数関数から線形関数へ変換する。この手順を以下の好適な例を用いて例示する。Ｔ_ｃの後の通過時間に対する式（５）の関連項の対数は、以下のとおり与えられる。

ここで、ｋ_ｅｐ ^ｉ，ｊは、ｋ_２ ^ｉ，ｊ／ｖ_ｅ ^ｉ，ｊに等しい。

Ｋ_{ｔｒａｎｓ} ^ｉ，ｊは、Ｆ^ｉ，ｊ・Ｅから推定される。Ｋ_{ｔｒａｎｓ} ^ｉ，ｊ＝ｋ_２ ^ｉ，ｊの場合、血管外細胞外体積割合ｖ_ｅ ^ｉ，ｊも推定され得る。

本発明の方法において、ｃ^ｉ，ｊ _ｐ［ｔ］の第１パス部分を特定することによって改善された推定値が見つけられ得る。簡単な検索手順を使用して、ｃ^ｉ，ｊ _ｐ［ｔ］の第１パス最大点の後のｃ^ｉ，ｊ _ｐ［ｔ］の第１の局所最小値が見つけられる。この局所最小値は、ｃ^ｉ，ｊ _ｐ［ｔ］の第１パスセグメントに対する時間の上限とされる。そこで、ｖ_ｐ ^ｉ，ｊは、単に、開始時点から局所最小時点までのｃ^ｉ，ｊ _ｐ［ｔ］の時間積分値である。

この手順の好適な実施形態は、以下のように例示される。ＭＲ記録の構成に基づいて、第１のパスパルスの持続時間の推定値が選択され得る（例えば、２５個の時間サンプル）。次いで、これらの（２５−Δ）個のサンプル内のボクセル別のｃ^ｉ，ｊ _ｐ［ｔ］に対する絶対最大値が見つけられる。Δは、通常２〜８の範囲にある。これは、サンプル時間ｔ_ｍａｘである。次に、範囲［ｔ_ｍａｘ，２５］内のボクセル別のｃ^ｉ，ｊ _ｐ［ｔ］に対する絶対最小値が見つけられる。これは、サンプル時間ｔ_ｍｉｎである。この最小値は、ボクセル内の第１のパスパルスの終点として得られる。そこで、血漿体積は、単にｔ＝１からｔ＝ｔ_ｍｉｎまでのｃ^ｉ，ｊ _ｐ［ｔ］の時間積分値である。

最後に、中央体積定理を使用して、Ｔ_ｃ ^ｉ，ｊ＝ｖ_ｐ ^ｉ，ｊ／Ｆ^ｉ，ｊからＴ_ｃ ^ｉ，ｊが推定され得る。

一旦これらのパラメータのうちのいずれか１つの値が決定されると、そのパラメトリック画像（または、「マップ」）が当該技術分野において公知の標準方法のいずれかによって生成され得る。個々のパラメトリック画像は、例えば本発明の装置によって、別々にかつ自動的に生成され得る。本方法は、実験ノイズによって著しく影響されず、かつこれらのパラメータを正確に決定することが可能なので、生成される画像は、従来技術と違って、高品質である。

本発明の一実施形態において、低速水交換方式がａａＪＷモデルに組み込まれ得る（Ｔ．Ｙａｎｋｅｅｌｅｖ ｅｔ ａｌ、ＩＳＭＴＭ ｃｏｎｆｅｒｅｎｃｅ、Ｓｅａｔｔｌｅ、２００６、Ｉｎｃｏｒｐｏｒａｔｉｎｇ ｔｈｅ ｅｆｆｅｃｔｓ ｏｆ ｔｒａｎｓｃｙｔｏｌｅｍｍａｌ ｗａｔｅｒｅｘｃｈａｎｇｅ ｉｎ ｔｈｅ ｒｅｆｅｒｅｎｃｅｒｅｇｉｏｎ ｍｏｄｅｌ ｆｏｒ ＤＣＥ−ＭＲＩ ａｎａｌｙｓｉｓ）。この実施形態において、さらに平均細胞内滞留時間のボクセル別の値およびそのパラメトリック画像を生成することが可能である。また、この拡張によって、特に対象領域が腫瘍である場合、ボクセル別の値およびその他のパラメータに対するパラメータマップの精度が改善される。

本発明の方法において、ブラインド逆畳み込み方法を使用して、各ボクセルの動脈入力関数ｃ^ｉ，ｊ _ｐ［ｔ］およびそれに対応する拡大・縮小される組織残留関数ｕ^ｉ，ｊ［ｔ］が分離される。本発明の好適なブラインド畳み込み方法は、ルーシー−リチャードソン（Ｌｕｃｙ−Ｒｉｃｈａｒｄｓｏｎ）逆畳み込み方法、またはランドウェバー（Ｌａｎｄｗｅｂｅｒ）方法とウィーナー（Ｗｉｅｎｅｒ）フィルタリングとの組み合わせを含む。最も好適な実施形態において、ルーシー−リチャードソン逆畳み込み方法が適用される。

ルーシー−リチャードソン逆畳み込み方法は、灌流イメージング（特に、ＭＲ灌流イメージング）に適用される場合、他の逆畳み込み方法と比較して、いくつかの利点を有する。ルーシー−リチャードソン逆畳み込み方法は、正値性、エネルギー保存およびサポートの画像拘束条件を自動的に満たす。重要なのは、動脈入力関数および組織残留関数の両方に固有の正値性拘束条件によって、実質的に解空間が低減され、かつしたがってこの方法がノイズに対して頑強にされる。

ルーシー−リチャードソン逆畳み込みアルゴリズムは、ポアッソン分布ノイズの付加した画像に対して、空間ドメインにおいて機能する。ｕ^ｉ，ｊ［ｔ］およびｃ^ｉ，ｊ _ｐ［ｔ］に対する期待値最大化（ＥＭ）アルゴリズムは、以下のとおり与えられる。

ここで、ｓは、繰り返し回数である。ボクセルを表示する添え字（ｉ，ｊ）は、省略した。ＥＭアルゴリズムであるので、単独で使用される場合、最適値に収束することが保証される。最適値は、入力関数が与えられると、最大事後確率または尤度である。

大域最適値からのその距離は、動脈入力関数の初期推定値によって決定される。アルゴリズムは、繰り返しループにおいて２回使用され得る。最適値は、通常局所値になる。まず、拡大・縮小される組織残留関数が推定される。次に、この新しい推定値を使用して、同じアルゴリズムを用いて動脈入力関数の新しい推定値が見つけられる。繰り返しを開始するために、ｕ^０［ｔ］の初期値は、最初の４０〜１００個の時間フレームに対して一定に設定され、かつ後のフレームに対してはゼロに設定される。ｕ^０［ｔ］の積分値は、１．０に設定される。動脈入力関数に対して、すべてのｔでｃ_ｐ ^０［ｔ］＝ｃ_ｔ［ｔ］である。

流量によって拡大・縮小される残留関数の血管パラメータは、畳み込みの２つの因子の繰り返しノンパラメトリック推定が完了した後で推定され得る。しかし、また、これらのパラメータは、拡大・縮小される残留関数の各繰り返し更新の直後に推定され得る。次いで、パラメータ値を使用して、拡大・縮小される残留関数の更新されたパラメトリック推定値が構築され得る。このアルゴリズムは、ノンパラメトリックアルゴリズムよりもｌｏｇ尤度関数を低減する（おそらく、増加されたノイズ頑強性による）。

外側の繰り返しループは、４回以上の繰り返し、好ましくは６回以上の繰り返しを含み得る。例えば、４〜１００回、好ましくは４〜１０回、より好ましくは６〜１０回、および最も好ましくは６〜８回である。

２つの内側ループについての好適な繰り返し回数は、１回であるが、１〜１０回であり得る。

本発明の方法は、シングルパスデータおよびマルチパスデータに対して使用され得るが、マルチパスデータが好ましくは使用される。マルチパスとは、造影剤を対象組織に２回以上（例えば、２〜１０回）通すことを意味する。

本発明の方法は、灌流イメージングデータを処理して前方駆出率を生成する（薬物動態コンセンサスモデルおよびブラインド逆畳み込みアルゴリズムの使用を含む）際に使用され得る。これは、後述される。

本発明の一実施形態において、左心室のトレーサ濃度関数は、１つの十分に混合されたコンパートメントに対して、マルチパス動脈入力関数および左心室インパルス応答の畳み込みとしてモデル化され得る。このインパルス応答は、パルス周波数が既知の場合、前方駆出率評価に対する単一コンパートメント指数関数モデルによって補間され得る。これは、式（４）の薬物動態モデル（薬物動態コンセンサスモデル）が前方駆出率を推定するための単一コンパートメント指数関数モデルに対して心室内で有効であるというのに等しい。

したがって、左心室内では、二重ブラインド逆畳み込み方式を用いれば、マルチパス左心室入力関数および左心室指数関数インパルス応答関数を推定できることになる。唯一必要な変更は、式（４）において合成定数−ｋ_２／ｖ_ｅを単一定数１／τに置き換えることである。すなわち、以下に与えられる。

ここで、ｔは、サンプル番号である。パルス間隔がΔｔであり、かつサンプリング間隔がΔｓであるならば、駆出率ｅｆは、以下のとおりである。

ここで、τ_１＝τ・Δｓである。

対応するモデリングは、左心房に対して行われ得る。左心房も単一の十分に混合されたコンパートメントと考えられ得る。末梢静脈注入すると、右心房および右心室は、十分に混合されないことがあり得る。

上記または当該分野で公知のブラインド逆畳み込み方法のいずれかを薬物動態コンセンサスモデルとの組み合わせて使用して、駆出率を決定し得る。また、駆出率のパラメトリック画像は、本明細書中に記載されるような方法または当該分野で公知の方法のいずれかによって生成され得る。

本発明の方法を行うために、本明細書中のアルゴリズムにおいて使用されるパラメータの値がイメージングデバイスおよびイメージングシーケンスの特性ならびに必要に応じて造影剤にしたがって設定される必要がある。したがって、操作者がコンピュータに使用される造影剤の性質を指定するだけでいいように、適切なパラメータ設定がパラメトリック灌流画像生成を行うために使用されるコンピュータプログラム中へ予めロードされ得る。実際には、計算に必要なパラメータの適切な値が自動的に選択され得るように、機械の製造者によって供給されるプログラムは参照テーブルを含み得る。したがって、本発明の方法が与え得るパラメータの高品質な個々のパラメトリック画像が容易に生成され得る。

本発明の方法における使用に適用可能な灌流イメージング手法の例としては、磁気共鳴イメージング（ＭＲＩ）、医用超音波検査法（音波造影法）、陽電子断層撮影法（ＰＥＴ）およびコンピュータ断層撮影法（ＣＴ）がある。特に、好適なのは、ＭＲＩを使用することである。

ＭＲＩが適用されるならば、スピンエコー（Ｔ_１またはＴ_２イメージング）またはグラジエントエコー（Ｔ_２ ^＊イメージング）信号生成に基づいたエコープラナーシーケンスの使用が好適である。

当該技術水準において公知のいずれかの従来の手法を本発明の方法と組み合わせて、観測データの品質を改善し得る。例えば、解像度は、ゼロパディング（シンク補間およびゼロフィリングとしても知られる）によって改善され得る。ゼロパディングは、フーリエ変換する前にゼロを生データに付加することを単に含む。ゼロパディングを適用することによって、動脈入力関数と組織残留関数との間の畳み込みのラップアラウンド効果（ｗｒａｐ ａｒｏｕｎｄ ｅｆｆｅｃｔｓ）が起きることが回避される。そうしなければ、ラップアラウンド効果がフーリエドメインに生じる。

ＭＲ画像におけるエッジアーチファクトは、観測シーケンスの最後の部分をガンマ関数で近似することによって観測時間シーケンスデータセットをゼロに近づくように外挿することによって、最小化され得る。ゼロパディングおよびガンマ関数の使用の両方は、周知の方法であるが、本明細書中に記載の方法と組み合わせる場合は、特定の値を取る。

本発明の方法は、２次元および３次元空間記録に適用され得る。３次元データセットは、２次元データセット（「スライス」）に分割され得、この場合、時間が第３の次元である。一好適な実施形態において、計算中の種々の位置間の相互作用はない、すなわち、各「スライス」および各ボクセルに対する計算は、他のスライスまたは他のボクセルに対して独立である。

必要なサンプリングレートおよび時間サンプルの数が一連の因子に依存することは、当業者の知るところである。例えば、イメージング方法および使用される造影剤の性質、単独の造影剤「ボーラス」が投与できる程度、および対象組織の血管特性は、すべて影響を有する。公知のサンプリング定理（ナイキストーシャノンサンプリング定理など）を適用して、適切なサンプリングレートが確保され得る。したがって、当該技術水準におけるすでに利用可能な知見から、当業者は、対象領域および使用される灌流イメージング手法に基づいて、過度の負担なく、適切なサンプリングレートで取得される適切な数の時間サンプルを選択することができる。

一般に、好ましくは０．１〜５秒当たり１時間サンプルのサンプリングレート、より好ましくは１〜３秒当たり、より好ましくは１．５〜３秒当たり１時間サンプルのレートが適用され得る。利用可能な装置のイメージングおよび処理能力によるが、一般に、サンプリングレートが高くなるほど、得られる推定値はよりよくなる。

取られるサンプルの数は、１００〜１０００、より好ましくは１００〜５００、例えば、１５０、２００、３００などであり得る。一般に、あるサンプリングレートにおいて、血管特性が遅いほど、より多くのサンプルが必要とされる。例えば、血管腫瘍特性が低速の組織に対しては、１秒当たり１．５のサンプリングレートで３００個のサンプルを取れば、通常、ガドリニウムキレート造影剤を用いたＭＲＩイメージングに対しては非常によい結果が得られる。同様の状況下における血管特性が中高速の組織に対しては、１秒当たり１．５のサンプリングレートおよび１５０〜２００個のサンプルが適用可能である。これらの特性は、他の灌流イメージング手法に必要なパラメータの目安として使用され得る。当業者は、過度の負担なく、必要なパラメータを選択することが容易にできる。

本発明の一実施形態において、本発明の方法は、造影剤（好ましくは、ポーラスの形態）を予め投与された被験者に対して実施される。当業者は、適用される灌流イメージングの方法に基づいて適切な公知の造影剤を選択することができる。

細胞外液に分散する従来の造影剤（Ｇｄ ＤＴＰＡ、Ｇｄ ＤＴＰＡ−ＢＭＡおよびＧｄ ＨＰＤＯ３Ａなど）を本発明の方法において使用して、透過性表面積の積および細胞外血管外液体積ならびに上記のその他の生理学的パラメータの推定を可能にし得る。そのような剤は、ＭＲ造影剤製造者（例えば、シェーリング、アマルシャム、ナイコメド（Ｎｙｃｏｍｅｄ）など）の特許文献において十分に記載される。

血管内貯留型（ｐｏｏｌｉｎｇ）造影剤（すなわち、脈管構造に滞留する造影剤）は、透過性表面積の積および細胞外血管外液体積のいずれの推定も可能でない。なぜなら、血液プール剤は、毛細管壁を通過しないからである。しかし、上記のすべての他の生理学的パラメータは、血液プール造影剤（例えば、シェーリングのＶａｓｏｖｉｓｔ）を使用して推定され得る。

本発明を実現するために使用されるコンピュータソフトウェア製品は、物理的担体（ディスク、テープまたは他の記憶デバイスなど）の形態で提供されるか、または遠隔地から（例えば、インターネットなどのネットワークを介して）送信される信号として提供され得る。

ここで、本発明を以下の実施例および添付の図面を参照して説明する。

図１は、灌流記録からの解剖学的画像を示す。ａ）は、左後部に上皮内癌を有する左乳房の軸方向断面を示す。ｂ）は、左心室壁内に陳旧性梗塞を有する心臓の軸方向断面を示す。ｃ）は、左正常腎臓の冠状断面を示す。 図２は、乳管内上皮内癌（一様目盛）を示す。横座標：０〜１。縦座標：０〜２７６ｓ。ａ）は、腫瘍からの４８個のボクセル内の観測トレーサ濃度の平均を示す。ｂ）は、同じボクセルの推定トレーサ濃度を示す。ｃ）は、同じボクセルの推定動脈入力関数を示す。ｄ）は、同じボクセルの組織残留関数を示す。 図３は、左後部における上皮内癌を有する乳房を示す。ａ）は、組織１００ｍｌかつ１分当たりの血漿流量（単位：ｍｌ）（ウィンドウ：０〜２００ｍｌ）を示す。ｂ）組織１００ｍｌ当たりの血漿体積（単位：ｍｌ）（ウィンドウ：０〜１５ｍｌ）を示す。ｃ）は、組織１００ｍｌ当たりの血管外細胞外液体積（単位：ｍｌ）（ウィンドウ：０〜７ｍｌ）を示す。ｄ）１ｍｌ当たりの摂取率（単位：ｍｌ）（ウィンドウ：０〜０．４）を示す。ｅ）平均通過時間（単位：秒）（範囲：０〜５．４ｓ）を示す。ｆ）ピークまでの時間（単位：秒）（範囲：０〜５２ｓ）を示す。 図４は、陳旧性梗塞を有する心臓を示す（一様目盛）。横座標：０〜１。縦座標：０〜６５．５ｓ。ａ）は、正常な心筋層の６２個のボクセルの観測トレーサ濃度の平均を示す。ｂ）同じボクセルの推定トレーサ濃度を示す。ｃ）同じボクセルの推定動脈入力関数を示す。ｄ）同じボクセルの組織残留関数を示す。 図５は、側壁梗塞を有する左心室を示す。ａ）は、組織１００ｍｌかつ１分当たりの血漿流量（単位：ｍｌ）（ウィンドウ：０〜４００ｍｌ）を示す。ｂ）組織１００ｍｌ当たりの血漿体積（単位：ｍｌ）（ウィンドウ：０〜１２ｍｌ）を示す。ｃ）は、組織１００ｍｌ当たりの血管外細胞外液体積（単位：ｍｌ）（ウィンドウ：０〜８ｍｌ）を示す。ｄ）１ｍｌ当たりの摂取率（単位：ｍｌ）（ウィンドウ：０〜０．４）を示す。ｅ）平均通過時間（単位：秒）（範囲：０〜２．０８ｓ）を示す。ｆ）ピークまでの時間（単位：秒）（範囲：０〜２０ｓ）を示す。 図６は、正常な腎臓皮質を示す（一様目盛）。横座標：０〜１。縦座標：０〜１８６ｓ。ａ）は、皮質からの７３個のボクセルの観測トレーサ濃度の平均を示す。ｂ）同じボクセルの推定トレーサ濃度を示す。ｃ）同じボクセルの推定動脈入力関数を示す。ｄ）同じボクセルの組織残留関数を示す。 図７は、正常腎臓皮質を示す（なお、髄質および腎盤に対する結果は無効である）。ａ）は、組織１００ｍｌかつ１分当たりの血漿流量（単位：ｍｌ）（ウィンドウ：０〜４００ｍｌ）を示す。ｂ）組織１００ｍｌ当たりの血漿体積（単位：ｍｌ）（ウィンドウ：０〜２５ｍｌ）を示す。ｃ）は、組織１００ｍｌ当たりの（皮質糸球体ろ過速度（ＧＦＲ）中の）毛細管移動定数Ｆ・Ｅ（単位：ｍｌ）（範囲：０〜１２７ｍｌ）を示す。ｄ）１ｍｌ当たりの摂取率（単位：ｍｌ）（ウィンドウ：０〜０．４）を示す。ｅ）平均通過時間（単位：秒）（範囲：０〜４．０ｓ）を示す。ｆ）ピークまでの時間（単位：秒）（範囲：７．０〜１５．０ｓ）を示す。

実施例
患者に対する検査は、地元の倫理委員会によって承認された。左乳房に乳管内上皮内癌腫瘍を有する１人の女性（５７歳）の乳房を、うつぶせの姿勢で、１．５テスラ全身装置（インテラ（Ｉｎｔｅｒａ）、フィリップス（Ｐｈｉｌｉｐｓ））における乳房送信コイルおよび並列処理（センス）受信ボディコイルを用いて検査した。陳旧性心臓梗塞を有する１人の女性（４６歳）の心臓および１人の志願者（５４歳男性）の腎臓を、あおむけの姿勢で、３．０テスラ全身装置（ＧＥシグマエキサイト（ＧＥ Ｓｉｇｎａ Ｅｘｃｉｔｅ）、ミルウォーキー、米国）における８チャネル表面送信コイルおよび並列処理（センス）受信ボディコイルを使用して検査した。記録中の乳房の動きは、１ボクセルより小さかった。位置合わせは、不要であった。心臓および腎臓画像の位置ずれを最小限にするために、中程度に息を吸ってから息を完全にとめた状態で前造影画像シーケンスを記録した。数回強制的に息を吸って吐いてをしてから中程度に息を吸った状態で後造影画像シーケンスを記録した。必要に応じて、検査志願者は、記録中に数回息をし、かつ中程度に息を吸った状態で息をとめ続けた。検査中はずっと、鼻カテーテルによって酸素（３リットル・分^−１）を与えて、なるべく呼吸を低減した。

画像データは、オフラインコンピュータに転送され、後処理された。その手順は、Ｃ言語で実現された（いくつかの数値レシピ（Ｎｕｍｅｒｉｃａｌ Ｒｅｃｉｐｅｓ）ルーチンを含む）。フリーのＸＩＴＥパッケージ（オスロ大学インフォマティクス部、ノルウェー）およびＭＡＴＬＡＢ（マスワークス（Ｍａｔｈｗｏｒｋｓ）、ネイティック、ＭＡ、米国）に基づく画像処理ツールを画像ハンドリング、視覚化および評価のために使用した。動きのアーチファクトをなるべく低減するために、画像化ボリュームのすべての前造影および後造影画像は、１つのボリューム画像時間シーケンスにマージし、そして最大尤度位置合わせアルゴリズムを使用して自動的に平行移動して位置合わせした。平面における両方向の最大平行移動は、１０〜２０ボクセルであった。各画像ボリュームに対するすべての可能な平行移動は、現代のデスクトップコンピュータ上で１２〜２４時間の計算時間をかけて力任せに（ｂｙ ｂｒｕｔｅ ｆｏｒｃｅ）試験した。

呼吸の生じた心臓または腎臓の画像は棄却した。呼吸のない対応の新しい画像は、前後のボリューム画像を使用して線形補間によって生成した。位置合わせのためのボリューム対象領域は、手動で描いた。

位置合わせに続いて、同じバターワース低域通過フィルタリングを各時間シーケンスに対して行って、帯域外ノイズを低減した。カットオフ周波数は、−π／２およびπ／２ラジアンであり、かつ相対カットオフ周波数は、０：３ラジアンであった。

実施例Ｉ：記録
乳房 磁化率コントラスト増強灌流イメージングのために、マルチセクション２次元マルチショットエコープラナー二重エコーＴ２^＊−重み付け（ＳＴ１／Ｔ２^＊．ＰＥＲＦ）平行画像取り込みシーケンスを乳房の軸方向平面において使用した。すべてのセクションに対して、フリップ角度α＝４０度、ＴＲ＝１２１０ｍｓ、ＴＥ_１＝６．０ｍｓ、ＴＥ_２＝３０．０ｍｓである。２４個の平行セクションは、スライス厚が５ｍｍであり、かつ面内ボクセルサイズが１．１×１．１ｍｍ^２であった。これは、ボクセル体積が６．０５ｍｍ^３であることに相当する。取り込みマトリクスサイズは、８５個の位相エンコードステップおよび２５６個の周波数エンコードステップを有し、かつ最終の画像マトリクスは、２５６×２５６ピクセルであった。位相エンコードは、行方向に行った。

ＭＲ信号をトレーサ濃度に変換し、そしてＴ_１効果を除去するために、二重エコーＴ_２ ^＊重み付けシーケンスには、前コントラスト画像シーケンスが必要でない（後述を参照）。Ｇｄ−ＤＴＰＡ（Ｍａｇｎｅｖｉｓｔ）（０．１ｍｍｏｌ／ｋｇ）をパワーインジェクタ（２ｍｌ／ｓ、７．５ｓ）を使用してボーラス注入し、その直後に３０ｍｌの生理食塩水を流した後、ボーラスの通過を各セクションのその後の８０枚の画像（サンプリング間隔３．６３ｓ）によってモニターした。

心臓 緩和度コントラスト増強灌流イメージングのために、Ｔ_１重み付け反転回復グラジエントエコー（ＳＰＧＲ）平行画像取り込みシーケンスを心臓の軸方向平面内で使用した。すべてのセクションに対して、ＴＩ＝１０ｍｓ、ＴＲ＝６：３２ｍｓ、ＴＥ＝１：１９２ｍｓであった。６つの平行セクションは、スライス厚が８ｍｍ、スライス間隔が１７ｍｍ、かつ面内ボクセルサイズが１．２５×１．２５ｍｍ^２であった。これは、ボクセル体積が１２．５ｍｍ^３であることに相当する。取り込みマトリクスサイズは、１２８×１２８であり、かつ画像マトリクスサイズは、２５６×２５６であった。視野は、３２０×３２０ｍｍであり、かつ位相エンコードは、行方向に行った。

ＭＲ信号をトレーサ濃度に変換できるように、３つの短い前コントラスト画像シーケンス（各１０フレームであり、それぞれのフリップα＝１０、２５、３０度）を記録した。ガドジアミド（Ｇａｄｏｄｉａｍｉｄ）（Ｏｍｎｉｓｃａｎ）（０．１ｍｍｏｌ／ｋｇ）をパワーインジェクタ（４ｍｌ／ｓ、５ｓ）を使用してボーラス注入し、その直後に２０ｍｌの生理食塩水を流した後、ボーラスの通過を各セクションのその後の１１５枚の画像（サンプリング間隔１．７１ｓかつフリップ角度α＝２５度）によってモニターした。

腎臓 緩和度コントラスト増強灌流イメージングのために、Ｔ_１重み付けＬＡＶＡ画像取り込みシーケンスを冠状平面内で使用した。すべてのセクションに対して、ＴＲ＝２：６６８ｍｓ、ＴＥ＝１．１３２ｍｓであった。ＬＡＶＡは、両面内方向における高速低角度飽和回復グラジエントエコーシーケンス（センスファクター（Ｓｅｎｓｅ ｆａｃｔｏｒ）が２）である。１４個の平行セクションは、スライス厚が５ｍｍであり、かつ面内ボクセルサイズが１：３６７２×１：３６７３ｍｍ^２であった。これは、ボクセル体積が９：３５ｍｍ^３であることに相当する。取り込みマトリクスサイズは、１２８×１２８であり、かつ画像マトリクスサイズは、２５６×２５６であった。視野は、３５０×３５０ｍｍであった。同じセクションの２つの連続画像フレーム間のサンプリング間隔は、１．８５６ｓであった。マンガン硫酸塩を有する大きなパッドを表面コイルと体との間に使用して、誘電効果の劣化をなるべく低減した。

ＭＲ信号をトレーサ濃度に変換できるように、４つの短い前コントラスト画像シーケンス（各１０フレームであり、それぞれのフリップα＝４、８、１２、１５度）を記録した。低用量のガドブトロール（ｇａｄｏｂｕｔｒｏｌ）（Ｇａｄｏｖｉｓｔ）（０．０１ｍｍｏｌ／ｋｇ）をパワーインジェクタ（３ｍｌ／ｓ、５ｓ）を用いてボーラス注入し、その直後に２０ｍｌの生理食塩水を流した後、ボーラスの通過を各セクションのその後の１１２枚の画像（サンプリング間隔１．８５６ｓかつフリップ角度α＝１２度）によってモニターした。

実施例２：トレーサ濃度の推定
時間ドメインにおける組織残留関数のエイリアジング効果を回避するために、各トレーサ時間シーケンスの最後の２５個の時間サンプルについて、組織残留関数をゼロに設定した。各画像シーケンスのうちの最初の３つの画像を棄却して、安定な状態の信号強度を確保した。調査される器官に対して、血液および組織における緩和度を一定と仮定し、かつ細胞膜を横断する水交換が高速であると仮定した。ガドリニウムトレーサの緩和度値は、血管パラメータの定量化には必要でない。しかし、同じ検査中では、すべての組織に対して同じ緩和度が仮定されなければならない。

組織緩和度から、トレーサ濃度を計算した。ガドリニウムトレーサのガドジアミド（Ｏｍｎｉｓｃａｎ）の緩和度の値は、Ｂ_０＝３．０Ｔｅｌｓａで３．２ｓ^−１ｍｍｏｌ^−１である。これは、以下のトレーサ濃度を与える。

ここで、κは、定常磁場の所定の強度Ｂ_０におけるコントラスト媒体の緩和度の値である。

乳房 二重エコーT_２ ^＊重み付けイメージングによって、Ｔ_１に関する仮定なしに、緩和速度（Ｒ_２ ^＊［ｔ］＝Ｔ_２ ^＊−１［ｔ］）を直接定量することが可能になる。エコープラナーパルスシーケンスの記録信号ｓ［ｔ］は、以下によって与えられる。

式（１４）を操作すると、以下の横緩和度変化が得られる。

ここで、下付き文字のｐｒｅは、対応する前コントラスト信号を示す。面内ボクセル添え字（ｉ、ｊ）は、理解しやすくするために省略した。

心臓 緩和速度Ｒ_１［ｔ］＝Ｔ_１ ^−１［ｔ］は、以下の式にしたがって、記録信号（ｓ）から、各画像ｔ（ｔ＝１、２、．．．、Ｔ）におけるボクセルｉ，ｊ（ｉ＝１、２、．．．、２５６、ｊ＝１、２、．．．、２５６）ごとに計算した。

ここで、ａ＝Ｅｃｏｓα、Ｅ＝ｅ^{−Ｒ１・ＴＲ}、αは、フリップ角度、およびｎは、ｋ空間の中心に対応するｎ番目の位相エンコードラインである。ここでは、ｎ＝１だけを使用した。面内添え字（ｉ，ｊ）は、理解しやすくするために省略した。ＴＩは、反転時間であり、およびＴＲは、繰り返し時間である。

ｓ［０］＝ｃｏｎ・ρ・ｓｉｎ（α）は、フリップ角度依存性キャリブレーション定数であり、受信側ゲインｃｏｎおよびボクセルのプロトン密度ρを含む。受信側ゲインｃｏｎもまた、Ｂ_１場不均一性および空間的に変化する受信側感度により、ボクセル別である。ボクセル別ｃｏｎおよびρ値は、トレーサ実験の間、一定であると仮定した。推定ボクセル別ｃｏｎ値および式（１６）を使用して、各ボクセルの緩和速度Ｒ_１［ｔ］を式（１６）の線形化式を使用して推定した（Ｐｒｅｓｓ ｅｔ ａｌ．、Ｎｕｍｅｒｉｃａｌ Ｒｅｃｉｐｅｓ ｉｎ Ｃ）。

腎臓 緩和速度Ｒ_１［ｔ］＝Ｔ_１ ^−１［ｔ］は、以下の式にしたがって、記録信号ｓ［ｔ］から、各画像ｔ（ｔ＝１、２、．．．、Ｔ）におけるボクセルｉ，ｊ（ｉ＝１、２、．．．、２５６、ｊ＝１、２、．．．、２５６）ごとに計算した。

ここで、αは、フリップ角度であり、および面内添え字（ｉ、ｊ）は、理解しやすくなるように省略した。ｓ［０］＝ｃｏｎ・ρ・ｓｉｎ（α）は、フリップ角度依存性キャリブレーション定数であり、受信側ゲインｃｏｎおよびボクセルのプロトン密度ρを含む。受信側ゲインｃｏｎもまた、Ｂ_１場不均一性および空間的に変化する受信側感度により、ボクセル別である。ボクセル別ｃｏｎおよびρ値は、トレーサ実験の間、一定であると仮定した。値は、ネルサーおよびミード（Ｎｅｌｓｅｒ ａｎｄ Ｍｅａｄ）の多次元の滑降シンプレックス法を使用して推定した。

動脈入力関数 血液ｃ_ｐ［ｔ］だけをボクセルからの有する動脈入力関数の積分値は、画像化ボリュームにおける大きな動脈内に対象領域を描くことによって推定した。同じ患者を連続して検査した場合に測定される流速、移動、飽和およびＴ_２ ^＊効果が大きいと、大きな血管の動脈入力関数が大きくばらつく。したがって、同じ検査からのすべてのセクションに対して、１セクションからの１つの動脈入力関数を使用した。乳房および心臓記録に対して、マスクをそれぞれ肺動脈および左心室内に描いた。腎臓記録に対して、マスクを腰下行大動脈内に描いた。

実施例３：インビボ実験
調査した器官の概観を与えるために、造影剤あり（図１ａ、ｃ）または造影剤なし（図１ｂ）の灌流シリーズの平均時間値を取ることによって解剖学的画像を得た。

腫瘍を有する乳房 ａａＪＷモデルに対して、腫瘍の対象領域内の造影剤の観測時間経過、および対応する推定動脈入力関数および組織残留関数の畳み込みは、重なり合った（図２ａ、ｂ）。例外は、２つの曲線の最近の部分であり、そこでは推定造影剤濃度が急峻に低減した。この低減は、離散畳み込みの境界効果によって生じた。この重なり合いは、かなり実際に近い動脈入力関数および組織残留関数が見つかったことを示す（図２ｃ、ｄ）。

ａａＪＷモデルを使用すると、血漿流量および血漿体積（図３ａ、ｂ）は、正常乳房および上皮内癌と異ならなかった（表１；１ａ、１ｄ）。反対に、通過時間（図３ｅ）、摂取率（図３ｄ）、透過性表面積の積および血管外細胞外空間（図３ｃ）は、上皮内癌に対するよりも正常乳房組織に対するほうが明らかにより低かった（表１；１ａ、１ｄ）。これらの結果は、非常に漏れの大きい毛細管が乳房癌腫に存在するという観察結果と一致する。

最後に、ピークまでの時間は、正常乳房組織に対するよりも上皮内癌に対するほうがより短かった（表１；１ａ、１ｄ、図３ｆ）。これは、腫瘍細胞が動脈供給に対して逆行効果を有することを意味するに違いない。

ブラインド逆畳み込みを用いたコンセンサスモデルは、ａａＪＷモデルよりも血漿流量および血漿体積に対してずっと低い推定値を与えた。しかし、最も実質的な変化は、平均通過時間および毛細管移動定数において見られた（表１；１ａ、１ｂ、１ｄ、１ｅ）。これらの大きな違いは、おそらくコンセンサスモデルにおいて存在する簡略化によって引き起こされ得る近似誤差によって生じたものである。

ブラインド逆畳み込みを用いたコンセンサスモデルおよびそれを用いないコンセンサスモデルの血漿流量および通過時間は、反対方向に約２倍異なり、他方血漿体積は、同様であった（表１；１ｂ、１ｃ）。これらの違いは、ブラインドコンセンサスモデルの組織残留関数の実質的な付加的ばらつきを除去することによって最もよく説明され得る。遅延のＳ．Ｄ．／平均比においては、著しい違いはなかった。

コンセンサスモデルを使用して正常乳房および上皮内癌の血管パラメータ値を比較すると、それらの値における変化は漸次大きくなった（表１；１ｃ、１ｆ）。これは、癌腫が大きくなり、そしてその下の間質性組織を侵襲する場合に予期されるものである。

梗塞を有する心臓 ａａＪＷモデルを使用した心臓に対しても、腫瘍の対象領域内の造影剤の観測された時間経過、およびそれに対応する推定動脈入力関数および組織残留関数の畳み込みは、よく重なり合った（図４ａ、ｂ）。この観察結果は、かなり実際に近い動脈入力関数および組織残留関数が見つかったことを示す（図４ｃ、ｄ）。

ａａＪＷモデルを使用すると、血漿流量、血漿体積、摂取率（図５ａ、ｂ、ｄ）、透過性表面積の積および血管外細胞外体積（図５ｃ）は、正常心筋層におけるよりも陳旧性梗塞におけるほうが明らかにより小さかった。反対に、ピークまでの時間は、正常心筋層におけるよりも陳旧性梗塞領域におけるほうが実質的により長かった（表２；１ａ、１ｄ、図５ｆ）。通過時間においては、違いがなかった（図５ｅ）。血管パラメータにおけるこれらすべての変化は、正常心筋層および陳旧性瘢痕組織の微小血管特性における差異と一致した。

平均通過時間は、コンセンサスモデルに対するよりもａａＪＷモデルに対するほうが明らかにより短く、かつ対応する標準偏差もより小さかった（表２；１ａ、１ｂ、１ｄ、１ｅ）。これらの大きな違いは、おそらくコンセンサスモデルにおいて存在する簡略化によって引き起こされ得る近似誤差によって生じたものである。左心室のマスク内（４５０ボクセル）の前方駆出率の平均値は、０：２１±０：００であった。

腎臓皮質 ａａＪＷモデルに対して、皮質の対象領域内の造影剤の観測された時間経過、およびそれに対応する推定動脈入力関数および組織残留関数の畳み込みは、完全に重なり合った（図６ａ、ｂ）。やはり、この観察結果は、かなり実際に近い動脈入力関数および組織残留関数が見つかったことを示す（図６ｃ、ｄ）。

ａａＪＷモデルを使用すると、平均糸球体ろ過速度は、５６ｍｌ／組織１００ｍｌ／分であった（図６ｃ）。大人の健常な男性の腎臓のそれぞれの重さは、約１５０ｇであり、かつ皮質は、腎臓の約７０％を構成する。これらの値および組織密度１．０３ｇ／ｍｌを使用すると、総糸球体ろ過速度は、１１４ｍｌ／分となった。

短い平均通過時間は、腎臓動脈から腎臓皮質へ血液が通過することによって生じる動脈入力関数の付加的なばらつきを除去することによって説明される（図７ｅ）（表３；１ａ）。小さな血漿体積は、内部に中心体積を使用することによって生じる（表３；１ａ、図７ｂ）。

摂取率（図７ｄ）は、使用される逆畳み込みまたはパラメータ推定方法にあまり依存しないほうがよい。なぜなら、漏れ項および血管内血漿率項は、血漿流量および動脈入力関数によって同程度に影響されるからである（式（１、２））。腎臓皮質の透過性表面積の積は、大きい（表３；１ａ、図７ｂ）。

実施例４：左心室の研究
ＳｏｎｏＶｕｅ（登録商標）をボーラス注入により投与された２人の志願者の左心室を標準のＧＥの超音波スキャナおよび標準のフィリップス（Ｐｈｉｌｉｐｓ）の超音波スキャナを使用して検査した。ＧＥのスキャナを使用して４回スキャンおよびフィリップスのスキャナを使用して１回のスキャンを行い、そしてデータセットを記録した。

左心室に対する前方駆出率は、記録されたデータセットから、ブラインド逆畳み込みおよび薬物動態コンセンサスモデルを使用し、およびまた、心収縮期および心拡張期における心内膜境界を描く公知の方法を使用することによって計算した。前方駆出率を計算するための本発明の方法を使用すると、すべての記録に対してよい一致が得られた。

マルチパス動的コントラスト増強灌流イメージングにおけるボクセル別動脈入力関数および組織残留関数のブラインド逆畳み込みのための方法が初めて提示された。本発明によって、血管パラメータの推定を線形、簡単かつノイズに対して頑強にすることが可能にされる。本発明の方法によって、公知の方法よりも臨床的により実際に近い血管パラメータを導出できる。

表１ａ。乳房の灌流パラメータ。平均＋Ｓ．Ｄ．。血漿流量（単位：ｍｌ／組織１００ｍｌ／分）。ａａＪＷモデルを用いた血漿体積（単位：ｍｌ／組織１００ｍｌ）。研究１ｂ：コンセンサスモデルを用いた正常乳房組織の結果。研究１ｃ：コンセンサスモデルおよび公知の動脈入力関数を用いた正常乳房組織の結果。研究１ｃ：コンセンサスモデルおよび公知の動脈入力関数を用いた正常乳房組織の結果。研究１ｄ：ａａＪＷモデルを用いた乳管内癌腫腫瘍上皮の結果。研究１ｅ：コンセンサスモデルを用いた乳管内癌腫腫瘍上皮の結果。研究１ｆ：コンセンサスモデルおよび公知の動脈入力関数を用いた乳管内上皮内癌腫瘍の結果。グループ１における標準偏差は、正常乳房の２２０個のボクセル内のばらつき、および乳管内上皮内癌腫瘍の４８個のボクセル内のばらつきを示す。１．０３ｇ／ｍｌに等しいと仮定される組織密度は、グループ２、３、４および６の組織１００ｍｌ当たりに与えられる単位に換算する際に使用される。毛細管においてヘマトクリットを使用して推定したグループ５における血漿流量およびグループ６における血漿体積は、０．７５である。

表１ｂ 乳房の灌流パラメータの続き。透過性表面積（ＰＳ）の積（単位：ｍｌ／組織１００ｍｌ／分）。血管外細胞外液体積（ＥＥＶ）（単位：ｍｌ／組織１００ｍｌ）。ＦＥ（単位：ｍｌ／１００ｍｌ／分 組織）。ピークまでの時間（単位：秒）。

表２ａ 心臓の灌流パラメータ。表１と同じ単位。研究１ａ：ａａＪＷモデルを使用した、左心室側壁に梗塞を有する１人の心臓の正常心筋層における結果。研究１ｂ：コンセンサスモデルを使用した、左心室側壁に梗塞を有する１人の心臓の正常心筋層における結果。研究１ｃ：公知の動脈入力関数を用いたコンセンサスモデルを使用した、左心室側壁に梗塞を有する１人の心臓の正常心筋層における結果。研究１ｄ：ａａＪＷモデルを使用した、１人の心臓における左心室側壁の梗塞心筋層における結果。研究１ｅ：コンセンサスモデルを使用した、１人の心臓における左心室側壁の梗塞心筋層における結果。研究１ｆ：公知の動脈入力関数を用いたコンセンサスモデルを使用した、１人の心臓における左心室側壁の梗塞心筋層における結果。グループ１における標準偏差は、正常心筋層の６２個のボクセル内のばらつき、および大きな梗塞を有する心臓の８５個のボクセル内のばらつきを示す。その他の標準偏差は、患者間のばらつきを表す。１．０３ｇ／ｍｌに等しいと仮定される組織密度は、グループ２、３、４、６および７の組織１００ｍｌ当たりに与えられる単位に換算する際に使用される。毛細管においてヘマトクリットを使用して推定したグループ５における血漿流量およびグループ６における血漿体積は、０．７５である。

表２ｂ 心臓の灌流パラメータの続き

表３ａ ヒトの腎臓皮質の灌流パラメータ。平均＋Ｓ．Ｄ．。表１と同じ単位。皮質マスクにおける７３個のボクセル。研究１ａ：ａａＪＷモデルを用いた正常ヒト腎臓皮質の結果。研究１ｂ：コンセンサスモデルを用いた正常ヒト腎臓皮質の結果。研究１ｃ：公知の動脈入力関数を使用したコンセンサスモデルを用いた常ヒト腎臓皮質の結果。

表３ｂ ヒトの腎臓皮質の灌流パラメータの続き。糸球体ろ過速度（ＧＦＲ）（単位：ｍｌ／組織１００ｍｌ／分）。

Claims (28)

1. 灌流イメージングデータを処理する方法であって、ブラインド逆畳み込みアルゴリズムならびにジョンソン（Ｊｏｈｎｓｏｎ）およびウィルソン（Ｗｉｌｓｏｎ）の断熱近似モデル（ａａＪＷ）の使用を含む、方法。
2. 灌流イメージングの方法であって、ブラインド逆畳み込みアルゴリズムおよびａａＪＷモデルの使用を介して、一続きの時間値［ｔ］において、被験者の対象領域の画像を前記対象領域のボクセルから生成することを含む、方法。
3. 灌流画像の画像処理の方法であって、ブラインド逆畳み込みアルゴリズムおよびａａＪＷモデルの使用および画像の生成を含む、方法。
4. 前記ブラインド逆畳み込みアルゴリズムおよび前記ａａＪＷモデルは、パラメータであるボクセル別動脈入力関数ｃ_ｐ［ｔ］およびボクセル別組織残留関数ｒ［ｔ］の生成において使用される、請求項１〜３のいずれか１つに記載の方法。
5. 前記方法は、ボクセル別パラメータである血漿流量、血漿体積、平均通過時間、摂取率（ｅｘｔｒａｃｔｉｏｎ ｆｒａｃｔｉｏｎ）、透過性表面積の積、血管外細胞外液体積、およびピークまでの時間のうちの少なくとも１つの生成をさらに含む、請求項１〜４のいずれか１つに記載の方法。
6. 血漿流量、血漿体積、平均通過時間、摂取率、透過性表面積の積、血管外細胞外液体積およびピークまでの時間が生成される、請求項５に記載の方法。
7. 前記パラメータのうちのいずれか１つの個別の画像が生成される、請求項４または５に記載の方法。
8. 請求項４または６に記載のすべてのパラメータの個別の画像が生成される、請求項４または６に記載の方法。
9. 前記ブラインド逆畳み込み方法が繰り返し使用される、請求項１〜８のいずれか１つに記載の方法。
10. 前記ブラインド逆畳み込み方法を繰り返し使用して、前記ボクセル別動脈入力関数ｃ_ｐ［ｔ］および拡大・縮小される（ｓｃａｌｅｄ）ボクセル別組織残留関数ｕ［ｔ］の値を得る、請求項９に記載の方法。
11. 前記繰り返しループは、４回以上の繰り返し回数を含む、請求項１０に記載の方法。
12. 前記繰り返しループは、６〜１０回の繰り返し回数を含む、請求項１１に記載の方法。
13. ａａＪＷを使用して、前記拡大・縮小されるボクセル別組織残留関数ｕ［ｔ］からボクセル別組織残留関数ｒ［ｔ］を生成する、請求項１０〜１２のいずれか１つに記載の方法。
14. 前記ボクセル別動脈入力関数ｃ_ｐ［ｔ］の初期推定値は、観測されるボクセル別トレーサ濃度ｃ［ｔ］である、請求項１０〜１３のいずれか１つに記載の方法。
15. 前記ボクセル別動脈入力関数ｃ_ｐ［ｔ］の初期推定値は、第１のパス動脈入力関数の積分値である、請求項１０〜１３のいずれか１つに記載の方法。
16. 前記拡大・縮小されるボクセル別組織残留関数ｕ［ｔ］の積分値の初期推定値は、１である、請求項１０〜１５のいずれか１つに記載の方法。
17. 前記灌流イメージング方法は、磁気共鳴イメージング（ＭＲＩ）、医用超音波検査法（音波造影法）、陽電子断層撮影法（ＰＥＴ）、およびコンピュータ断層撮影法（ＣＴ）である、請求項１〜１６のいずれか１つに記載の方法。
18. 前記灌流イメージング方法は、磁気共鳴イメージング（ＭＲＩ）である、請求項１７に記載の方法。
19. 前記ブラインド逆畳み込み方法は、ルーシー−リチャードソン（Ｌｕｃｙ−Ｒｉｃｈａｒｄｓｏｎ）逆畳み込み方法、またはランドウェバー（Ｌａｎｄｗｅｂｅｒ）方法とウィーナーフィルタリングとの組み合わせである、請求項１〜１８のいずれか１つに記載の方法。
20. 前記ブラインド逆畳み込み方法は、ルーシー−リチャードソン逆畳み込み方法である、請求項１９に記載の方法。
21. 低速水交換方式モデルがａａＪＷモデルに組み込まれる、請求項１〜２０のいずれか１つに記載の方法。
22. 平均細胞内滞留時間のボクセル別の値が生成される、請求項２１に記載の方法。
23. 前記ボクセル別平均細胞内滞留時間の個別の画像が生成される、請求項２２に記載の方法。
24. 前記灌流イメージングデータは、マルチパスデータである、請求項１〜２３のいずれか１つに記載の方法。
25. 請求項１〜２４のいずれか１つに記載の方法において使用されるデータ処理装置であって、前記装置は、一続きの時間値（ｔ）に対して対象領域のボクセルについての１セットの灌流イメージング信号値を入力として受け取る手段と、ブラインド逆畳み込みアルゴリズムならびにジョンソンおよびウィルソンの断熱近似モデル（ａａＪＷ）の使用を含む、前記１セットの灌流イメージング信号値からパラメータを生成する手段とを有する、装置。
26. 命令を搭載するコンピュータソフトウェア製品であって、前記命令は、請求項１〜２４のいずれか１つに記載の方法において使用されるデータ処理装置上で実行される場合に、前記装置が、一続きの時間値（ｔ）に対して対象領域のボクセルについての１セットの灌流イメージング信号値を入力として受け取り、かつブラインド逆畳み込みアルゴリズムならびにジョンソンおよびウィルソンの断熱近似モデル（ａａＪＷ）の使用を含む、前記１セットの灌流イメージング信号値からパラメータの値を生成する手段を有するように構成されるようにする、コンピュータソフトウェア製品。
27. 灌流イメージングデータを処理して、前方駆出率を生成する方法であって、薬物動態コンセンサスモデルおよびブラインド逆畳み込みアルゴリズムの使用を含む、方法。
28. 左心室の前方駆出率が生成される、請求項２６に記載の方法。

Images