JP6373388B2 - 灌流指数を推定する方法 - Google Patents

Description

本発明は、例えば血流といった、哺乳動物の灌流指数を推定する方法、対応する医療画像システムおよび本発明を実行可能な対応するコンピュータプログラム製品に関する。

血液が組織を流れる過程は、灌流と呼ばれ、最も基本的な生理学で定量化可能なものの一つである。灌流の不調は、哺乳動物を病気や死に至らせる過程である。

正常な脳機能は、活動の代謝要求を満たすために酸素の連続供給を必要とする。脳組織内における酸素の局部的有用性は、慣例的に脳血流（ＣＢＦ）の量および動脈血中の酸素濃度から推定される。ＣＢＦは、神経血管結合として知られる神経活動の局部レベルに影響されやすく、したがってＣＢＦの変化を検出する方法は、脳機能のマッピングの有力な手段を提供するものである。

神経細胞の壊死や永続的な梗塞の発症に関して明確に定義された閾値により、疾患において、ＣＢＦの測定は、急性虚血性脳卒中の疑いのある患者および頸動脈の流れを制限する狭閉塞疾患の患者の評価に広く利用されている。

現在、ＣＢＦを推定する標準的な方法は、いくつかの変形において利用可能な特異値分解（ＳＶＤ）である。要約すると、ＣＢＦを推定するＳＶＤ技術は本質的に、（ＣＢＦにより乗じられる残留関数と定義される）インパルス応答を推定するため、特異値分解（ＳＶＤ）を用いる濃度時間曲線と動脈入力関数（ＡＩＦ）の解析に基づくものである。その最大値は、それゆえＣＢＦである。

しかし、ＳＶＤは、高い流れ成分を過小評価し、動脈入力関数（ＡＩＦ）における遅延に影響されやすいことが知られており、これはＳＶＤにより取得された灌流指数の生理学的解釈を損なう可能性がある。これは、虚血性領域を描写するためＣＢＦが用いられる、例えば急性脳卒中といった、深刻な状態にある患者の診断にとっては、言うまでもなく、非常に危険なものとなる可能性がある。

最近、Kim Mouridsen他著, Neurolmage, 33（2006）570-579ページが、毛細管内の通過の生理学的モデルを用いる代替モデルのアプローチを提案した。ベイズモデルを用いて、血行力学パラメータの推定値を取得することができる。特に、ガンマ関数を有する残留関数のパラメトリックモデルが用いられる。取得された結果は、低流量のＳＶＤ方法によくたとえられるが、高流量を過小評価せず、またＳＶＤのようなノイズの影響が比較的ない。しかし、結果はまだノイズにやや影響されやすく、また推定値の精度を損ない、臨床的に実行不可能な計算時間という結果に終わる可能性がある一般的な数値近似スキームに依存しているため、実用的な実装は制限されている。

したがって、改善された灌流指数を推定する方法は有益であり、灌流指数を推定する、より効率的および／または信頼性のある方法は、特に有益である。

本発明の目的は、従来技術の代替物を提供することである。

特に、ノイズに対して不十分なロバストネスを有する灌流指数および／または長時間の画像／測定を必要とする灌流指数を有する上記の従来技術の問題を解決する、灌流指数を推定する方法を提供することを、本発明の目的として見ることができる。

したがって、上記の目的および他のいくつかの目的は、プロセッサにより哺乳動物における灌流指数（ＰＩ）を推定する方法を提供することによって、本発明の第一の態様において取得されることが意図されており、前記方法は、注入された造影剤の時間の関数としての造影剤濃度を表すデータ（ＤＡＴ，ＤＡＴ’）を取得することと、取得したデータに対して、時間の関数としての毛細管通過時間分布のパラメトリックモデル（ＰＭ）を適用することと、最適化によって灌流指数を推定するための正則化を施して、好ましくは期待値最大化（ＥＭ：Expectation-Maximization）型法であるマイナライズ最大化（ＭＭ：Minorize-Maximization）型法を適用することと、
を備え、好ましくは前記期待値最大化型法である前記マイナライズ最大化型法は、最大化ステップにおいて用いられる前記造影剤濃度の非線形観測モデルにおける造影剤の観測誤差の分散（Ｃ_ε）の適切な解析表現を有する方法。

本発明は、特に、しかし非排他的に、より効果的および／または確実な方法で灌流指数を取得するのに有益である。本発明の有益性は、特に、
・最新鋭のＳＶＤ推定と比較して、実際の患者の灌流指数の改善された推定と、
・少ない画像取得物で虚血性の病変が検出可能なことによる、例えば急性脳卒中におけるスキャン時間の低減と、
・神経学的評価と急性脳卒中における画像との間のより良い対応関係と、
・毛細管通過時間分布の、より生理学上実行可能でロバストネスのある推定値と、
・統計の利用性をより高める、毛管流レベルの検出のこの技術の高い識別能力による、臨床試験に必要とされる患者数の低減と、
・灌流指数における、造影剤到達の遅延の影響の受けにくさと、
を含む。

これらの、また他の利点については、以下に、より詳細に説明および図示される。

本発明の文脈において、注入された造影剤の時間の関数としての造影剤濃度を表すデータを取得することは、造影剤自体を注入する工程とは別個独立していることに留意すべきである。本発明は、造影剤を注入する工程を含むことを意図しておらず、実用的な目的のため、そのようなデータの取得は、一般的に哺乳動物への造影剤の注入が試験されてからかなりの時間の後、行われる。造影剤の注入は、一般的に、標準的な医療手順であり、重大な健康リスクと関連付けられるものではないことにも留意すべきである。

灌流マーカまたは血行力学パラメータとも呼ばれる灌流指数は、例えば脳血流（ＣＢＦ）、脳血液量（ＣＢＶ）または平均通過時間（ＭＴＴ）といった従来の指数に限定されないだけでなく、例えば同質ではない毛細管分布の血流を表す毛細管通過の不均一性（ＣＴＨまたはＣＴＴＨ）または毛細管から抽出することができる酸素の最大量に関する酸素抽出量（ＯＥＣ）といった、適切に定義されてもよい他の灌流指数でもあることが、本発明の文脈にある。

本発明は、例えば前臨床試験および臨床目的の両方の試験のために、ヒト、サル、ウマ、ウシ、ブタ、げっ歯類等を含むが、これらに限定されない、いずれかの哺乳動物の灌流指数を推定するために適用されてもよい。

本発明は、体内において灌流が起こる脳、腎臓、肝臓、心臓、筋肉の一部等を含むが、これらに限定されない、身体のいずれかの部分の灌流指数を推定するために適用されてもよい。

本発明の文脈において、「推定する」という用語は、広く解釈されるものとして理解されるものであり、１以上の灌流指数の直接的および／または間接的な測定値を得ることを含んでもよい。これは、特に、１以上の指数のおおよその測定値を含む。

本発明の文脈においては、「時間の関数としての毛細管通過時間分布のパラメトリックモデル（ＰＭ）」という用語は、１以上のパラメータによる毛細管通過分布の生理学的モデル、すなわち時間の経過とともに毛細管領域を流れる血液に関するモデルが選択されると、広く解釈されるものとして理解される。実際のモデル選択は、データを取得するための医療用画像モダリティ、疾患の種類および／または画像される組織のタイプ、計算時間／精度／正確さ／ロバストネス／収束性等を含むが、これらに限定されない多くの要因に依存する。いくつかの用途のために、いくつかの性質の異なるモデルは、適用および／または比較されてもよいよう意図されている。

本文脈において、「毛細管領域」、「毛細管構造」および「毛細血管床」という用語は、器官または組織の指定された部位に供給をする毛細管が絡み合う網状組織を指す。毛細血管床は、本発明の文脈において、毛細血管床を通る造影剤濃度を表すデータを得るために適用される手段の性質および／または毛細血管床における造影剤濃度を表すデータを測定するための手段の性質および／または測定される組織に依存する様々な空間の拡張を有していてもよい。

毛細血管床または構造は、マイクロメートル範囲の基本寸法を有する毛細管の網状組織からなり、一般的な脳の毛細管は、例えば１２０マイクロメートルの長さ、および８マイクロメートルの直径を有している。毛細血管床の延長は、したがって、造影剤濃度を表すデータの有意義な測定値を引き出すための複数の毛細管の測定の必要性により、下から制限される。同様に、医療用画像技術を扱う当業者によって容易に理解されるように、毛細血管床の拡張は、一般的に、造影剤濃度を表すデータの測定に適用される測定手段の利用可能な空間解像度により制限される。例えば二光子顕微鏡画像といった、いくつかの光検出技術は、厳密には、回折限界（〜λ／２，λはプロービング放射線の波長）のみによって制限され、この制限は、一般的に単一の毛細管の空間拡張より小さいため、より多くの毛細管の測定または平均を取ることが必要となるのに対し、磁気共鳴画像（ＭＲＩまたはＭＲ）で標本を加工することでサブミリメートルオーダーの空間解像度が生まれる。

数値解析の中で、マイナライズ最大化（ＭＭ）型問題は、通常、対応するメジャライズ最大化（ＭＭ）型問題として再定義できることが、よく知られている。本文脈において、本発明は単に、例えば期待値最大化（ＥＭ）といったマイナライズ最大化型モデルとして明確にするために考案されるものだが、数値解析の当業者に容易に理解されるように、本発明の教義および原理は、メジャライズ最小化（ＭＭ）型のモデルもまた同様に含むことができる。この二重性の詳細については、当業者は、その全てを参照により本明細書に援用されるKenneth Lange著, Numerical Analysis for Statisticians, 2010年, Springer（189-221の特定ページ）を参照する。

期待値最大化（ＥＭ）法は、閉形式の尤度方程式は利用できないが、最尤の反復推定が利用可能な場合に、特に有用である、特別な種類のマイナライズ最大化（ＭＭ）法として数値解析から知られている。ＥＭ法の中核には、欠損データ（記録失敗または理論的欠落）の概念がある。ＥＭアルゴリズムの期待ステップを欠損データの「充填」とする考え方がある。代理関数は、一般的に尤度自体よりも見つけ出すことが容易であるため、多くの場合、Ｍステップは、解析的に行うことができる。しかしながら、代償は、反復が実行されなければならないことである。本発明は、以下にさらに詳細に説明するように計算を著しく改善し、観測誤差の適切な解析表現が最大化ステップで使用されている点で独特である。一般的なＥＭ法のさらなる詳細については、当業者は、その全てを参照により本明細書に援用されるKenneth Lange著, Numerical Analysis for Statisticians, 2010年, Springer（特定の189-221ページ）を参照する。

数値解析において、不良設定問題の解決または過剰適合を防ぐために、追加の制約情報を提供することにより、いわゆる正則化を頻繁に行う必要がある。いくつかの実施形態にかかる本発明については、正則化なしで機能する可能性があることが意図されているが、本発明は、正則化を用いても実施される。一般に、ベイズモデルは、モデルパラメータに課せられた制約によって正則化を適用してもよい。特に、正則化は、灌流パラメータの事前の期待値と分散を介して促進されてもよい。これらの数量は、灌流測定値から経験的に取得することができるため、これは、残留関数の振幅が制限されるＳＶＤに標準的に実行されるより、生理学的に、正則化へのより適切なアプローチであってもよい。

本発明の文脈において、「適切な解析表現」という用語は、広く解釈されるべきであると理解されるものである。この用語は、特に直接的かつ明確に分散値を生成することができる造影剤の観測誤差の分散の閉形式を意味するが、これに限定されるものではない。これは、近似値のみを提供する四則演算、または反復法によってのみ取得することができる値の陰関数記述とは通常反対のものである。したがって、本発明は、分散値を生成する数値近似法、反復法またはさらなる数学的手順を必要としない。

本発明の原理および教義を考慮すると、医療従事者や医療用画像業者により容易に理解されるように、本発明の文脈において、「造影剤」という用語は、広く解釈されるべきであり、特に「ボーラス」、「ボーラス剤」、「トレーサ」および／または「トレーサ剤」という用語またはこれらの技術的等価物のいずれかを含むが、これに限定されるものではないことをさらに理解されるものである。

本発明の文脈では、「造影剤の漏出」という用語は、例えば、脈管が無傷ではない時に脈管構造から脈管外および細胞外空間へ造影剤が漏れることであると広く解釈されるべきであるとさらに理解されるものである。これは、例えば、血液脳関門の破損、穿刺または穿孔、または内皮の堅い接合状態の損失である可能性がある。また、血液脳関門が存在しない、中枢神経系の外部組織の動態にも関連する。血管領域への脈管外および細胞外空間からの逆流については、臨床測定の基準に比べて、はるかに長い時間範囲で起こるものであるため、ほとんどの実施形態に含まれない。例えば、逆流はおよそ５〜１０分以上であるが、典型的な短いＭＲＩスキャンは、およそ１分である。

一の実施形態において、時間の関数としての毛細管通過時間分散のパラメトリックモデル（ＰＭ）は、少なくとも２つのパラメータ、任意で少なくとも３つのパラメータを有し、前記モデルは、さらに正またはゼロの通過時間に対して定義されてもよい。あるいは、本発明に、２つのパラメータのみ適用されてもよく、さらにまた、３つのパラメータのみ適用されてもよい。少数のパラメータは、同時に生理学の現実的なモデルでありながら、シンプルで、それゆえ高速な計算の実装を容易にするものである。特に、時間の関数としての毛細管通過時間分散のパラメトリックモデル（ＰＭ）は、ガンマ分布、Ｆ分布、逆ガウス分布、ガウス分布、半正規分布、レイリー分布、ワイブル分布、負の二項分布、ポアソン分布、ボルツマン分布、またはそれらの任意の組み合わせ、またはそれらから導かれる任意の等価物からなる、非網羅的なグループから選択されてもよい。

有利な実施形態では、毛細管通過時間分布のパラメトリックモデル（ＰＭ）は、好ましくは、Ｒが残留関数で、αおよびβが２つのパラメータである同次多項式：
で、残留関数Ｒ（ｔ）の推定を可能にするように、時間の関数（Г（ｔ））としてのガンマ関数であってもよい。

いくつかの実施形態においては、経時的な造影剤濃度の非線形観測モデル（ｙ_ｔ）は、第一要素（ｆ）および第二要素（ε）を有してもよく、初項は、好ましくは脳血流（ＣＢＦ）または他の身体部分の血流と、残留関数と畳み込まれる組織濃度（Ｃ_ａ）と、の積を表し、第二項は、その観測誤差（ε_ｔ）を表しており、下記の式（２．３）参照、この分離は改善されたモデリングを容易にする。１つの変形例において、生理学上、研究中の哺乳動物および／または潜在的な疾患に関連する場合、経時的な造影剤濃度の非線形観測モデル（ｙ_ｔ）は、脈管構造から脈管外および細胞外空間への経時的な造影剤の漏出要素を表す第三項をさらに備えてもよい。

特に、トレーサ溢出（漏出）の影響を受けない灌流モダリティのロバストな推定値を取得する能力は、血液脳関門を破壊する可能性のある腫瘍、脳水腫および多発性硬化症などの病的状態に必要である。周辺組織へのトレーサの損失をモデル化し損ねてしまった場合、灌流の定量化にバイアスが生じる。

特に、非線形観測モデルは、θが、好ましくは脳血流（ＣＢＦ）である血流、パラメトリックモデルのパラメータ（α，β）、およびＣ_ａの測定部位とｙ_ｔの測定部位との間の遅延を表す時間遅延（δ）の関数である未知パラメータである、同次多項式：
を有する。

遅延は、一例として、例えば患者といった特定の哺乳動物の従前のデータに基づいて、および／または同一、類似または関連する疾患を持つ多数の被験者の集団から取得されたデータに基づいて、推定、算出および／または予測することができる。本発明が、この遅延の影響を組み込むのに、特に有利であることが実証されている。図４Ａおよび対応する以下の説明を参照。

有利なことに、
である場合に、非線形観測モデルの線形化は、
の通りに行われる。

また、ｙの分布が、以下に説明される、より簡単な後続のモデリングを促進する可能性のある統計的線形モデルと定義されてもよいように、ｙが造影剤濃度曲線を表し、
である場合、残差ｒは、
と定義されてもよい。

有利には、統計的線形モデルの１以上のパラメータは、尤度関数または事後確率等の灌流に関する関数を最適化するため、好ましくは期待値最大化（ＥＭ）型法であるマイナライズ最大化（ＭＭ）型法を用いて推定されてもよい。有利なことに、上述の関数は、事後確率であってもよく、好ましくは期待値最大化（ＥＭ）型法であるマイナライズ最大化（ＭＭ）型法の収束が、非最適定常状態に対するアルゴリズム収束が得られたときに、特にパラメータの平均値が更新される、統計的線形モデルのパラメータの事前分布の選択を評価するために監視されてもよい。この方法は、事前のパラメータの再初期化も示す。

有利には、マイノリゼーションは、好ましくは共同対数尤度関数である関数の期待値を用いて、より好ましくは期待値最大化（ＥＭ）型法：
を用いて取得されてもよい。特定の好ましい実施形態において、造影剤の観測誤差の分散（Ｃ_ε）は、灌流指数の改善された推定に適用されてもよい、少なくともＪ、Ｃ_θ｜ｙおよびに残差ｒに関数的に依存してもよい。

より有利には、造影剤の観測誤差の分散（Ｃ_ε）は、Ｃ_ε＝σ^２Ｉという仮定のもと、Ｎが、本発明により取得される灌流指数の速度および品質に大きな影響を与える、観測の適切な解析表現を生む造影剤濃度測定値の数である場合に、
と記述することができる。複数の分散要素が存在する場合、Fisherのスコアリング法を適用することができる。

いくつかの実施形態において、本方法は、最近、生理学的にかなり重要であることが実証された、ＣＴＨまたはＣＴＴＨ等の哺乳動物の毛細血管床における血流の不均一性を表す灌流指数を推定してもよい。他の実施形態において、本方法は、哺乳動物の毛細血管床における血液からの、好ましくは酸素（ＯＥＣ，ＯＥＦ）および／またはグルコースまたは他の生化学物質である物質の抽出能力（ＥＣ）を表す灌流指数を推定してもよい。

本発明は、平均通過時間（ＭＴＴ）、好ましく脳血流（ＣＢＦ）である血流、好ましくは脳血液量（ＣＢＶ）である血液量および漏出（Ｋ_{ｔｒａｎｓ}またはＫ２）からなる、非網羅的なグループから１以上の灌流指数を推定してもよい。

本発明の第一の態様の変形は、哺乳動物における灌流指数を推定する方法の提供に関するものであって、前記方法は、注入された造影剤の時間の関数としての造影剤濃度を表すデータを取得することと、時間の関数としての毛細管通過時間分布のパラメトリックモデル（ＰＭ）と、灌流指数を推定する正則化とともに、好ましくは期待値最大化（ＥＭ）型法であるマイナライズ最大化（ＭＭ）型法と、を適用することによって、灌流指数を推定することを備え、好ましくは前記期待値最大化型法である前記マイナライズ最大化型法は、前記最大化ステップにおいて用いられる前記造影剤濃度の非線形観測モデルにおける造影剤の観測誤差の分散（Ｃ_ε）の適切な解析表現を有する。

第二の態様において、本発明は、医療用画像システムに関するものであり、前記システムは、注入された造影剤の造影剤動態を表すデータと、任意で前記データを記憶する記憶手段とを取得する画像モダリティと、第一の態様にかかる灌流指数を推定するよう準備されたプロセッサと、を備える。

哺乳動物の毛細管構造における造影剤濃度を表すデータの取得に適した医療用画像システムのいずれも、本発明の全般的な教義および原理内に適用されることが意図されている。

画像モダリティは、ＤＳＣ−ＭＲＩやＤＣＥ−ＭＲＩ（心臓ＭＲＩ）などの磁気共鳴映像法（ＭＲＩ）、陽電子放射断層撮影法（ＰＥＴ）、単光子放射断層撮影法（ＳＰＥＣＴ）、コンピュータ断層撮影法（ＣＴ）、超音波法、二光子画像法を含む多光子分光法、直交偏波スペクトル（ＯＰＳ）画像法、サイドストリーム暗視野（ＳＤＦ）画像法、共焦点型画像法、レーザードップラー測定、網膜のレーザ走査検眼、毛細管顕微鏡検査、およびこれらの任意の組合せからなる、非網羅的なグループから、特に選択されてもよい。

第三の態様において、本発明は、接続されたデータ記憶手段を有する少なくとも１つのコンピュータを備えるコンピュータシステムに、本発明の第二の態様にかかる関連する哺乳動物における灌流指数を推定する画像モダリティおよび／または医療用画像システムを制御させる、コンピュータプログラム製品に関する。

本発明の態様は、特に、しかし非排他的に、コンピュータシステムにダウンロードまたはアップロードされた場合に、本発明の第二の態様の装置／システムの操作をするためのコンピュータシステムを可能にするコンピュータプログラム製品によって、本発明が完成される場合に有益である。そのようなコンピュータプログラム製品は、いずれの種類のコンピュータ可読媒体またはネットワークを通して備えられてもよい。

本発明の個々の態様は、他のいずれの態様と組み合わせられてもよい。本発明のこれらの、また他の態様については、説明される実施形態に照らし、以下の説明にて明らかにされる。

本発明にかかる方法は、添付の図面に関連して、より詳細に説明される。図面は、本発明の実施の１つの方法を示すものであり、添付の特許請求の範囲に記載されない、他の実施形態の可能性を限定するものとして解釈されるものではない。

デジタルファントムの正確なＣＴＴＨおよびＯＥＦ^ｍａｘマップの図である。さらに示されているのは、ＳＮＲの２つのレベルに対応する算出されたマップのカラー版である。ＯＥＦ^ｍａｘマップは無次元である一方、ＣＴＴＨマップのカラーバーの単位は秒である。 図１Ａに対応するグレートーン版である。 ＳＮＲ＝２０およびＳＮＲ＝１００の、提案されたパラメトリック、ｓＳＶＤおよびｏＳＶＤ法により算出されたＭＴＴマップのカラー版を示す。正しいＭＴＴ値の正方形からなるマップは、左側に示される。 図２Ａに対応するグレートーン版である。 ＳＮＲの関数として、ｓＳＶＤ、ｏＳＶＤおよびパラメトリックモデル（Ｐａｒａｍ）といった、異なるモデルそれぞれのＭＴＴの平均二乗誤差、ＭＳＥを示すグラフである。３０秒の最大ＭＴＴ値が、実施されている。挿入図は、対数スケールでのパラメトリックモデルの収束を示す。 本発明にかかるパラメトリック（実線）を用いた遅延の関数としてのＲＭＳＥとＳＶＤ（２種類の破線）とを示すグラフである。ＣＢＦが相対的な尺度であるのに対し、ＯＥＦ^ｍａｘ値は、１０００倍に拡大されている（テキスト参照）。遅延は、秒で示されている。 パラメトリック（実線）およびＳＶＤ（２種類の破線）を用いたＡＩＦのＦＷＨＭ（「幅」）の関数としてのＭＳＥを示すグラフである。ＣＢＦが相対的な尺度であるのに対し、ＯＥＦ^ｍａｘ値は、１０００倍に拡大されている（テキスト参照）。ＦＷＨＭは、秒で示されている。 事前の共分散マトリックスでの縮尺変更の影響を示すグラフである。灌流観測値のＭＳＥにおける急激な増加が、反復アプローチ（破線）で観測されるのに対し、本発明による正確な最大化（実線）は、縮尺変更へのロバストネスを提供している。ＯＥＦ^ｍａｘ値は、このグラフで視覚的に比較するために、１０００倍に拡大されたという点に留意されたい。 臨床的に取得された一連の脳卒中データの灌流マーカを示す。右に示されるカラースケールは、行のすべての画像に関するものである。 カラー版の図６Ａに対応するグレースケール版である。 漏出期間を含む、本発明にかかるアルゴリズムを造影剤漏出のレベルが異なるファントムに適用することで得られるＭＴＴマップを示す。 カラーの図７Ａに対応するグレースケール版である。 本発明の第二の態様における、医療用画像システム１００の概略図である。 本発明にかかる方法またはコンピュータプログラム製品の詳細の概要を表す概略フローチャートである。

図８は、本発明の一の態様における、医療用画像システム１００の概略図である。システム１００は、例えばその頭部を模式断面図で見られるヒトといった、哺乳動物２００ＭＡＭに注入された造影剤の時間の関数としての造影剤濃度を表すデータを取得可能な、例えば本明細書に示されるように、ＤＳＣ ＭＲＩモダリティといった医療用画像モダリティ１０ＩＭを有する。例えば脳といった、１以上の哺乳動物の部位を画像することができる。ＤＳＣ ＭＲＩは、一般的に、強磁場ＭＲスキャナ１２と、対応する計測・制御システム１１とを有する。スキャナは、計測・制御システム１１に伝達される一次データＤＡＴを取得し、その後の工程で計測・制御システム１１から二次データＤＡＴ´がプロセッサ２０ＰＲＯＣに伝達される。

プロセッサ２０ＰＲＯＣにおいて、注入された造影剤の時間の関数としての造影剤濃度を示す取得されたデータＤＡＴおよびＤＡＴ´を用いて、哺乳動物２００ＭＡＭの灌流指数を推定する方法が実施される。灌流指数は、例えば脳内といった毛細管構造における時間の関数としての毛細管通過時間分布に、特にパラメトリックモデルＰＭ２１を適用することで推定される。さらに、好ましくは期待値最大化ＥＭ法２２であるマイナライズ最大化ＭＭ法は、以下で、より詳細に説明されるように用いられる。灌流指数を推定する際には、正則化も行われる。

特に、好ましくは期待値最大化ＥＭ法であるマイナライズ最大化ＭＭ法は、矢印２３により模式的に示されるとおり、最大化ステップにおいて用いられる造影剤濃度の非線形観測モデルにおける、造影剤の観測誤差の分散Ｃ_εの適切な解析表現を有する。

プロセッサは、模式的に示されるように、結果として得られた灌流指数ＰＩの画像を表示可能なように、コンピュータスクリーンまたはモニタ等のディスプレイ装置３０ＤＩＳＰに使用可能に接続されてもよい。あるいは、または加えて、灌流指数ＰＩは、後の解析および診断目的のために、適した種類のいずれかの記憶装置４０ＳＴＯに伝達されてもよい。

２．理論
２．１．トレーサ動態モデル
定められた組織ボクセルに関して、ｔが時刻を示す場合、脈管内トレーサは、Ｃ_ａ（ｔ）と示される特徴的な動脈血濃度を持つ毛細血管床に運ばれることを前提としている。インジケータ希釈理論を用いて、時刻tにおける毛細管内のトレーサの濃度は、時刻tにおいて毛細管中に依然として存在するトレーサの割合を表す、残留関数Ｒ（ｔ）と動脈供給の畳み込みに比例する。つまり、
である。定数κは、細動脈および毛細管の微小血管ヘマトクリットレベル、および脳組織の密度に依存するが、本実施形態では一定と仮定する。

残留係数のパラメトリックの推定を可能にするため、毛細管通過時間に関するガンマ分布を採用する。
ＡＩＰ測定の部位に関する毛細血管に対するトレーサ送達における、いずれの遅延にも対応するため、ＡＩＦタイミングのシフトを可能にするための追加のパラメータがモデルに含まれている。すなわちＣ_ａ（ｔ）→Ｃ_ａ（ｔ−δ）である。

２．２．パラメータ推定
経時的な造影剤濃度の測定に関し、（２．１）の右手側を表すｆ、および未知パラメータθ＝（ＣＢＦ，α，β，δ）を表すθを用いて非線形観測モデルを仮定する。
観測誤差εが共分散Ｃ_εのゼロ−平均ガウス分布であると仮定している。パラメータベクトルについても、平均η_θおよび共分散マトリックスＣ_θを用いて、確率的であると仮定される。パラメータ推定は、ｙ＝（ｙ１，．．．，ｙＮ）の周辺分布ｐ（ｙ）を用いて頻度論的視点から行なわれる。（２．４）における、観察されていない確率的要素による、期待値最大化（ＥＭ）アルゴリズムを用いて自由パラメータを推定することは、自然である。Dempster他著（1977）参照。以下において、θへの事前分布を許可する、ランダム効果の現在の推定値の周りに一般的に適用される線膨張に基づいて、ＥＭ型アルゴリズムを採用する。Mouridsen他著（2006）参照。従前の方法とは対照的に、本発明は、Ｍステップにおける正確な最大化に関し、陽関数表示の数式を適用する。詳細については、アペンディックスＡにて説明される。

モデルを線形にするため、
である場合、ランダム効果（２．３）、（２．４）の現在の推定値
についてのｆの一次膨張が、
と考える。

ｙが組織濃度曲線を表す場合、残差は、
であり、以下の式が成り立つ。

εおよびε_θが互いに独立していると仮定すると、データの周辺分布は、
である。効果的には、これは、（２．６）の共分散構造における第二項がランダム効果に起因する、統計的線形混合モデルの特別な場合である。しかし、パラメータ推定値の正則化を可能にするため、θの超パラメータＮ（η_θ、Ｃ_θ）が与えられ、Ｎが正規分布である、ベイズアプローチが用いられる。実施において、それらはモデル自由アプローチによって決定される。セクション３．２を参照。ＥＭ型アルゴリズムは、誤差共分散Ｃおよびθの事後平均値η_θ｜ｙと共分散Ｃ_θ｜ｙを推定するため、当然適用される。

Ｅステップにおける平均値は、観察結果ｙが与えられた場合のθのガウス事後密度であるｑ（θ｜ｙ）を基準にして取られる。θの事後モーメントは、
である。Ｍステップにおいて、最大化される式は、
である。

一度に単一のボクセルのみが分析されるため、誤差共分散Ｃ_εの簡単な式がパラメータの識別可能性を確保すると仮定されなければならならず、結果的にＣ_ε＝σ２Ｉと考えられる。ここでＮが、観察の回数である場合、分散要素の推定が以下の式の通りであるように、Ｍステップを分析的に行えることがアペンディックスＡに示される。

Mouridsen他著（2006）に示される通り、θの事後密度を最大化するニュートンアルゴリズムで事後並数を推定するため、ＥＭアプローチの等価性を利用することで、膨張点は、レーベンバーグ・マーカート減衰率を使用する正則化を用いてリセットされる。

３．物質および方法
３．１．酸素抽出能力の計算
毛細管通過時間と最大酸素抽出割合との関係ＯＥＦ^ｍａｘは、JespersenおよびOstergaard著（2012）において、より詳細に説明される。ｈ（τ）が毛細通過時間分散（２．２）であり、Ｑが、通過時間τの単一の、毛細管の酸素抽出量である場合、ＯＥＦ^ｍａｘは、
と表される。Ｑは、通常の微分方程式の数値解法として取得される（詳細には、JespersenおよびOstergaard著（2012）を参照）。ここでｋが毛細管壁の酸素に対する浸透性を表し、Ｃは、酸素の濃度であり、τが毛細管通過時間である場合、
である。残りの定数は、JespersenおよびOstergaard著（2012）から採用されたものであり、α_Ｈ＝３．１×１０^−５ｍｍＨｇ^−１，Ｂ＝０．１９４３ｍＬ／ｍＬ，Ｐ_５０＝２６ｍｍＨｇおよびＣ_ｔ＝２５ｍｍＨｇである。速度定数ｋは、本明細書の臨床の脳卒中データにおける、正常に見える白質におけるＯＥＦ^ｍａｘ＝３０％を要求することで定められる。セクション３．５を参照。結果の数値は、ｋ＝３８ｓ^−１となる。

３．２．事後平均値および共分散の特定
ＣＢＦおよびＭＴＴの計算上効率的な推定値は、特異値分解（ＳＶＤ）を用いて取得することができる。これらの推定値にはバイアスがかかっていることが知られているが、Mouridsen他著（2006）、それらは、実行可能な最初の近似値を提供することができ、したがって、本実施形態において事前確率を集中させるのに用いられる。遅延は、残留関数のノンパラメトリック推定の最大値に対する時間Ｔｍａｘとして推定される。事前確率に関し、残留関数が指数であるとして、β＝ＭＴＴを意味する（２．２）においてα＝１を設定する。

ノンパラメトリック推定値の不正確さを認識し、ＥＭアルゴリズムの推定値を用いるフィッティングの質は、二乗平均平方根誤差値によって監視される、すなわち、ｙ_ｉがデータポイントｉを表し、
が対応する推定値である場合、以下の通りとなる。
二乗誤差基準が、最終ＥＭ推定に関して定められた値（例えば、ε_ｔｈｒ＞０．０３）を越えた場合、αがα＝１および９＝ＭＴＴと設定される一方で、ＣＢＦおよびδの事前平均値を最近の推定値にリセットする。そしてＥＭ法が繰り返される。ε_ｔｈｒの直接の最適化がなかったとしても、このアプローチはうまく動作する。

パラメータの事前共分散すなわちＣ_θに関して、対応するパラメータの、一般的な、相対的な大きさを反映する対角行列が用いられる。具体的には、Ｃ_θ＝ｄｉａｇ（０．１，１，１，１０）を設定する。Mouridsen他著（2006）にある通り、対数スケールおよびＣ_θにおける分散にパラメータのモデルを作り、対数変換パラメータに適用する。セクション３．４において、Ｃ_θのミススペシフィケーションに対するロバストネスを試験する。

３．３．シミュレーション
信号の強さ曲線は、Mouridsen他著（2006）にある通りに（２．１）を用いて生成され、造影剤濃度と横緩和速度Ｒ^＊ _２との間の直線関係を仮定した。特に、ＡＩＦは、標準的な臨床注入スキームの一般的な形状および偏角を持つ動脈濃度時間曲線を生成するため、ａ＝１，ｔ_０＝１０，ｂ＝３およびｃ＝１．５である以下のガンマ変量関数を用いてシミュレートされた。Mouridsen他著（2006）。
ＣＢＶ＝４％に関して、シミュレーションが行われた。信号曲線は、Ｓ（０）＝１００およびＴ_Ｅ＝５０ｍｓであるとして、ｒ_２が横緩和速度である場合の以下の式を用いて生成された。
式（３．５）では、薬剤の脈管内濃度と横緩和速度との間の直線関係を仮定した。６０ｍｌ／１００ｇ／ｍｉｎのＣＢＦおよび４％のＣＢＶが、正常な灰白質には典型的な基準線に対して４０％の信号低下を生成するよう、ＤＳＣデータにて一般的に観察される信号低下を生成するため、定数κｒ_２は、Mouridsen他著（2006）にあるように較正された。ＡＩＦに関して、κｒ_２は、基準線に対して６０％信号低下を生成するよう固定された。

ゼロ平均を有するガウスノイズは、特定の基準線ＳＮＲｓを生成するため、信号曲線（３．５）に追加された。セクション３．４を参照。サンプリングレートは、ＴＲ＝１．５ｓに固定され、１０．５ｓ基準線を含む９９ｓの取得時間が、シミュレートされた。

３．４．デジタルファントム画像
式（３．１）における毛細管通過時間分布ｈ（τ）は、その平均値、ＭＴＴ＝αβおよび標準偏差ＣＴＨ＝β√αにより完全に特定される。ＭＴＴおよびＣＴＨを個別に推定することができるかを評価するため、グリッド上でこれらの数量が変化するデジタルファントムを生成した。これは、現在の研究では、ＭＴＴおよびＣＴＨ値のさらに広い範囲が考察されるものであるが、ＯＥＦ^ｍａｘの通過時間分布への依存関係を示すJespersenおよびOstergaard著（2012）における図４Ａを模倣するものである。各ファントムは７×７の正方形からなり、すべての正方形は、ＭＴＴおよびＣＴＨの、ある組み合わせに対応する。各正方形は、そのすべてが同一のＭＴＴおよびＣＴＨを持つが、異なるノイズ認識を持つ、１４×１４のボクセルに細分される。ＭＴＴは、ｘ方向に沿って２ｓから２０ｓに変化し、ＣＴＨは、ｙ方向に沿って２ｓから２０ｓに変化した。正方形は、２つのボクセルにより分離された。ＣＴＨおよびＯＥＦ^ｍａｘの真の値は、図１「正確」の左の行に示され、真のＭＴＴの値は、図２「正確」の左の行に表示される。

ファントムは、１０から５００の間のＳＮＲで生成された。加えて、ＡＩＦの開始と組織曲線の開始Ｃ（ｔ）との間の時差を導入することで、動脈トレーサ到達遅延がシミュレートされた。遅延は、０ｓから１０ｓの間で変化した。１ｓから６ｓの間で変動する半値全幅（ＦＷＨＭ）の値を取得するため、速度パラメータ（ｃ）が変化する一方、ＡＩＦ（３．４）における形状パラメータをｂ＝３．０に固定することで、ＡＩＦの測定部位から観察濃度への通過の間の分散がシミュレートされた。

ノンパラメトリックＳＶＤが、事前分布の中心として寄与する血流力学パラメータの適切な初期推定を提供する一方、事前の共分散の決定は困難である。全長の流れ（ＣＢＦ）と通過時間分布（α）とを共に表す、ＣＢＦの尺度および通過関数パラメータαを変化させることで、ミススペシフィケーションに対するロバストネスを試験する。ＣＢＦの分散が１０の因数で拡大され、αの分散が同じ因数で縮小される、代替の共分散マトリックス
を検討する。これは、ＣＢＦおよびαの尺度を効果的に交換するものである。以下の事前の共分散マトリックスを用いて、中間的な尺度変化のある共分散マトリックスに関するデジタルファントムにおけるアルゴリズム性能を求め、
また、１／２０単位で増分されるζ＝０からζ＝１の間で変動するζの値に対するファントムにおける二乗平均誤差を計算する。正確な最大化（２．１２）を持つ現在のアルゴリズムの結果を、Mouridsen他著（2006）における反復するFisherのスコアリング法と比較する。

３．５．患者データ
定められたアルゴリズムの性能と、標準的なＳＶＤ方法の性能とを比較するために、前部循環脳卒中である患者７名のＭＲＩが含まれた。データは、Ｉ−Ｋｎｏｗ ｓｔｒｏｋｅ ｉｍａｇｉｎｇ ｐｒｏｊｅｃｔの一部として取得された。プロトコルは、関連する国内倫理委員会によって承認され、インフォームドコンセントは患者の入会の前に取得された。標準ＤＳＣ−ＭＲＩは、ＧＥ社製Ｓｉｇｎａ Ｅｘｃｉｔｅ １．５、ＧＥ社製Ｓｉｇｎａ Ｅｘｃｉｔｅ ３、ＧＥ社製Ｓｉｇｎａ ＨＤｘ １．５、Ｓｉｅｍｅｎｓ社製Ｔｒｉ−ｏＴｉｍ ３、Ｓｉｅｍｅｎｓ社製Ａｖａｎｔｏ １．５、Ｓｉｅｍｅｎｓ社製Ｓｏｎａｔａ １．５、Ｐｈｉｌｉｐｓ社製Ｇｙｒｏｓｃａｎ ＮＴ １．５またはＰｈｉｌｉｐｓ社製Ｉｎｔｅｒａ １．５で行われた。ダイナミックエコープラナー画像法（３／１．５テスラの電界強度でエコー時間（ＴＥ）３０／５０ｍｓ、繰り返し時間（ＴＲ）１５００ｍｓ、視野（ＦＯＶ）２４ｃｍ、マトリックス１２８×１２８、スライス厚５ｍｍ）が、後に同じ速度での生理食塩水３０ｍｌの注入が続く、ガドリニウムベースの造影剤（０．１ｍｍｏｌ／ｋｇ）の５ｍｌ／ｓの速度での静脈注射中に行なわれた。経過観察の液体減衰反転回復（ＦＬＡＩＲ）画像は、最終的な梗塞の範囲の識別を可能にするために、１週間または１ヶ月後に取得された。

３．６．画像処理
ファントム画像ならびに患者データの分析は、「Ｐｅｎｇｕｉｎ」と呼ばれる、自ら開発した灌流分析ソフトウェア（ＡＩＦ選択、パラメトリックおよびノンパラメトリック・デコンボリューションおよび画像表示ツールを備える）で行なわれた。ソフトウェアは、オーフス大学の神経科学機能統合センター（ＣＦＩＮ）機関で開発された。臨床データにおいて、動脈入力関数は、熟練の放射線科医師によって半自動的に識別された。血行力学モダリティのパラメトリック推定は、セクション３．２で説明された、定められた値で行われた。標準的なＳＶＤは０．２の閾値に適用され、ｏＳＶＤとともに、振幅指数が０．０９５に設定された。ＳＶＤ法の遅延は、残留関数の最大化の時間（Ｔｍａｘ）として推定された。

４．結果
４．１．特定のＳＮＲｓでのＣＴＨおよびＯＥＦ^ｍａｘの推定
図１は、ＳＮＲ＝１００およびＳＮＲ＝２０のデジタルファントムにおける、ＯＥＦ^ｍａｘおよびＣＴＨの推定値を示す。真の値は、左の列のファントムに表示される。ＳＮＲ＝１００（中央の列）において、ＣＴＴＨまたはＣＴＨ、およびＯＥＦ^ｍａｘの推定値は、測定可能なバイアスおよび小さな視野内強度変化を実証する真の値と十分に一致している。

ＳＮＲ＝２０において、ＣＴＨ推定値は、高ノイズの状況にあっても水平方向における正方形に沿ったバイアスが低くかつ変化が小さく、ＭＴＴのＣＴＨを個別に推定することができることを示す。予想通り、正方形内の変化は、より低いＳＮＲで増加するが、ＣＴＨの差異は、ＭＴＴレベル全体にわたり識別できる状態で残る。より低いパネルに示される通り、非常に高い流れ（低いＭＴＴ）および高い通過時間の不均一性を有する場合のＯＥＦ^ｍａｘのいくつかのバイアス（過大評価）を除き、ＯＥＦ^ｍａｘ推定値の分散およびバイアスは、同様に機能した。

図３において、１０から５００の範囲におけるＳＮＲに関して、ＯＥＦ^ｍａｘおよびＣＴＨパラメータの相対的な二乗平均誤差（ＭＳＥ）が描かれている。両方のパラメータは、ＳＮＲの増進に関する高速収束を表示する。

４．２．特定のＳＮＲにおけるＭＴＴの推定
ＣＴＨおよびＯＥＦ^ｍａｘの推定値が、我々の技術に固有のものである一方、ＭＴＴ推定値は、ｓＳＶＤおよびｏＳＶＤと比較することができる。図２は、パラメトリックモデルによる、またＳＮＲ＝１００（最上部パネル）およびＳＮＲ＝２０（最下部パネル）の２つのＳＶＤ技術による、ＭＴＴ推定値を示している。左端に示される真のＭＴＴの値と比較すると、ＳＮＲ＝１００において、パラメトリックモデルが識別できるバイアスがないことと領域内の分散が小さいことが見てとれる。図１におけるＣＴＨの同様の結果と相まって、これは、ＭＴＴおよびＣＴＨの独立した推定値を提供するアルゴリズムの性能を解明するものである。比較において、ｓＳＶＤとｏＳＶＤの両方には、ＭＴＴ推定値における実質的なバイアスが現れている。ファントムの低い列のように、ＣＴＨが低い場合、このバイアスは顕著ではないが、適度なＣＴＨ（中間の行）および高いＣＴＨに対して増加し、ＭＴＴ推定値は、行全体にわたって均一に近い形で現れる。

ＳＮＲ＝２０において、パラメトリック法では、中程度のＣＴＨの値および高いＣＴＨ値に関して、ＭＴＴにおけるわずかなバイアスが現れ、これは、低いＭＴＴ値に関してより顕著なものである（左上角）。この領域における低いＭＴＴは、同時に高いＣＴＨが残留関数に長期にわたってゼロ以外を残すよう要求する間に、残留関数が、その最大値から急速な低下をしなければならないことを意味している。したがって、この領域におけるＭＴＴ過大評価は、低いサンプリングレートＴＲ＝１．５ｓによるものであってもよい。ＳＶＤ法は、非常に低いＣＴＨにあっても、実質的なバイアスを示す。

下部左側から上部右側の角までの対角線に沿う、基本的な残留関数は、指数関数的である。両方のＳＶＤ技術は、以前に公開された結果と一致する、この方向に沿ってMTTの合理的な推定値を与える。

４．３．全体的なＳＮＲのセンシティビティ
設定のより広い範囲にわたる灌流アルゴリズムの性能を比較するため、ＲＭＳＥによって、推定された値および真の値の間の一致を要約する。結果は、図３に示される。ＳＮＲが１０から５０に増大するのに従い、パラメトリックアルゴリズムのＲＭＳＥは０に向けて急激に減少する。ＳＮＲ＝５０において、パラメトリックアルゴリズムに関するＲＭＳＥは、ＳＶＤ技術のものよりもおよそ一桁ほど低い。ｏＳＶＤに関するＲＭＳＥがＳＮＲ＝５０まで減少するのに対し、ｓＳＶＤに関するＲＭＳＥは、増加するＳＮＲとともに増大しないが、ＳＮＲがさらに増加することによってわずかに増大する。

４．４．トレーサ遅延に対するセンシティビティおよびＡＩＭ測定部位からの分散
図４Ａは、固定されたＳＮＲ＝５０の、０から１０秒の間のトレーサ遅延に関する、全てのモダリティのＲＭＳＥまたは灌流指数ＰＩを示す。

一般に、ＲＭＳＥに関するパラメトリックマップは、遅延へのわずかな体系的依存を表す（図４Ａ参照）。同図は、またｏＳＶＤ、ＣＢＦおよびＭＴＴ値の遅延鈍感性を確認するものであるが、技術間におけるＣＢＦ推定誤差が類似しているのに対し、ＭＴＴのＲＭＳＥはパラメトリック推定値に比例して値が２倍される。パラメトリックアプローチにおいて非常にわずかな推定誤差が観測される一方、ｏＳＶＤ遅延測定値における誤差は、ここで抽出される遅延のすべての値に関して完全に一定している。ｏＳＶＤ方法により作成された遅延マップ（本明細書にマップは図示されない）を精査することで、ＭＴＴの増加により遅延が増加するという傾向が明らかになる。これは、ｏＳＶＤがＭＴＴと遅延とを的確に分離することができないことを示唆するものである。

図４Ｂは、パラメトリックモデルによって生成されるマップは、ＡＩＦにおける分散の度合いを和らげるのに全般的に鈍感であることを示すものであり、ボーラスが、脳に到達する前に、脈管の疾患によってさらなる分散を起こしている場合の、対象物２００における灌流推定値または灌流指数ＰＩを確実に推定することができることを示唆するものである。同様の結果が、精度の低いＭＴＴ推定値へ向かう傾向を除き、ｏＳＶＤに関して観測される。１秒から２秒へ向かうパラメトリックモデルに関して観測される大きな低下は、ＴＲ＝１．５ｓに固定されるサンプリングレートよりも低いＡＩＦの幅（ＦＷＨＭ）に関連する。ｏＳＶＤ方法については、この特徴を示すことが観測されない。

４．５．パラメトリックモデル：事前の共分散に関するロバストネス
図５において、ζ＝０は、予め用いられた事前の共分散マトリックスに対応し、ζ＝１は、ＣＢＦおよびαの尺度が１０の因数により変化した共分散マトリックスを表す。

正確な最適化アプローチは、事前の共分散マトリックスのミススペシフィケーションに対して、ζ＝０．５周辺のすべてのモダリティにわたってＲＭＳＥの急激な増加が観測される反復方法よりもロバストであることが見てとれる。

４．６．患者データ
ＡからＧに分類された７名の脳卒中患者の臨床データが表１に示される。図６は、パラメトリック法およびＳＶＤ技術、ならびに永続的な病変を示す経過観察のＦＬＡＩＲ画像からの灌流マップを示す。

患者Ａにおいては、すべてのメトリクスは神経病的欠損（ＮＩＨＳＳ＝１１）と一致する、灌流が不完全である大きな領域を例示しているが、パラメトリックＭＴＴおよびＣＴＨ画像が、灌流が不完全な組織を、周辺の影響のない組織から、より明確に区別している様子が見てとれる。患者Ｂにおいては、病変が、パラメトリックＭＴＴ、ＯＥＦ^ｍａｘおよびＣＴＨに比べて、ＳＶＤ ＭＴＴでは不明瞭に見える。特に、病変の前部の境界については、画像ノイズのため、ＳＶＤ画像では識別が困難となっている。

患者Ｃについては、ＳＶＤマップの、特に脳半球の反対側で、より大きな度合いで分散している高い鮮明度の領域が観察される。これらの領域は、ＳＶＤベース画像において灌流が不完全であると見られるが、デコンボリューションにより生成されたアーチファクトを反映していると思われる。これは、負の濃度曲線の部分をレンダリングし、ＳＶＤと同様に直接適合する場合に高振幅の残留関数をもたらす、低コントラスト対雑音比の結果（この患者で一般的に３から５）であってもよい。パラメトリックモデルは、ボーラスの導管からの意味のある値に適合可能であってもよく、構造上、負になることができない。この発見は、低いＳＮＲがＳＶＤ法に関するＭＴＴの著しい過大評価を引き起こすという、我々のファントムシミュレーションと一致する。図２を参照。

患者Ｄにおいては、ＳＶＤマップの鮮明度のレベルが低いことを示唆する可能性があるのに対し、とりわけ最終結果と一致する高い鮮明度量を表示するパラメトリックＭＴＴにおいては、モダリティ全体にわたる病変の推定値における一致度が低い。病変全体にわたってＣＴＨの上昇が見てとれるため、そのシミュレーション結果によれば、通過時間の不均一性によってＳＶＤ ＭＴＴマップが混乱させられたものと推測することができる。

ＳＶＤ法による灌流病変の明らかな過大評価は、患者Ｅにおいて、わずかにのみ上昇したパラメトリックＭＴＴをもって、より明白なものとなっている。小さな灌流病変は、低いＮＩＨＳＳの値と一致する。パラメトリック法およびＳＶＤ法はすべて、ＳＶＤ ＭＴＴ病変に対応する拡張（データは図示されない）をもって、ボーラス遅延の大きな領域を示す。これは、ＭＣＡ領域における効率的な二次的な供給によって、神経機能が維持されたことを示唆するものである。

患者Ｅとは対照的に、患者ＦおよびＧについては、パラメトリックＭＴＴ、ＯＥＦ^ｍａｘおよびＣＴＨでは病変が明確に見てとれるのに対し、ｏＳＶＤ ＭＴＴにて不完全な灌流を示していない。患者Ｇについては、ｓＳＶＤ ＭＴＴで病変を見てとることができるが、患者Ｆについては不明瞭である。両方の場合において、パラメトリック画像において見られる病変は、１ヶ月の経過観察ＦＬＡＩＲに見られる病変と一致する。基本的な信号曲線の精査は、灌流の不完全な組織における切り捨てられた信号曲線という結果をもたらす短いデータ取得時間が、この効果に寄与した可能性を明らかにした。パラメトリックモデルおよび標準ＳＶＤは、これらのケースでは、この現象に影響されないものと思われる。

５．考察
我々は、独創的な血行力学に関する数量、ＣＴＨおよび導き出される酸素抽出量、ＤＳＣ−ＭＲＩに基づくＯＥＦ^ｍａｘを推定する統計学技術を発明した。この方法は、巨視的な組織灌流指数ＣＢＦ、ＣＢＶおよびＭＴＴを推定する現存の技術を拡張するものである。我々のシミュレーション結果は、ＳＮＲの乏しいデータにおいてさえ、残留関数の広い範囲にわたるＣＴＨおよびＯＥＦ^ｍａｘ推定値におけるバイアスおよび分散が低いことを示すものである。重要なことには、これらの推定値は、急性脳卒中の画像での共通の懸念である、遅延や分散アーチファクトに対して比較的鈍感である。

我々は、ｓＳＶＤやｏＳＶＤと比較して、ＭＴＴ推定における実質的な改善を実証した。とりわけ、我々は、広く用いられるｓＳＶＤおよびｏＳＶＤデコンボリューションアルゴリズムによるＭＴＴ推定値は、低いノイズレベル（ＳＮＲ＝１００）であってもＣＴＨにより混乱させられる（図２）ことを実証した。

本発明の一つの主な動機となるのは、パラメトリックモデルの手段により毛管流分布の独創的なマーカの有用な推定が可能になることである。原理的に、ｓＳＶＤやｏＳＶＤ等のノンパラメトリックアプローチは、第２ステップにてそのようなパラメータが抽出できる残留関数を推定する。しかし、モデル独立技術が未知の残留関数のすべての時点を推定するのに対し、モデルベースのアプローチは関連する灌流量だけを推定し、結果、より有用な確実な推定をもたらすものである。デコンボリューション問題の不良設定具合（２．１）は、ＳＶＤ技術が残留関数における発振に関する様々な制約を通して実施される正則化を必要とする。本発明にかかる、定められたパラメトリックアプローチをもって、灌流指数の事前の推定および分散を通して正則化が促進される。これらの数量は、灌流測定から経験的に取得することができるため、これは、生理学的に正則化に対してより適切なアプローチであってもよい。

アペンディックスＡ．期待値最大化
このアペンディックスにおいて、上記のセクション２の推定アルゴリズムにおける、ＥおよびＭステップに関する正確な式を取得する。
以下のモデルにおいて、
η_θおよびＣ_θを固定された事前確率と見なしてＣ_εを推定し、式（２．５）における膨張点をリセットするため、事後平均値η_θ｜ｙを使用する。

ＥＭアルゴリズムは、θ|yの事後密度の最近の推定に関するθおよびｙの結合密度の推定値を最初に計算し、続いて、未知のパラメータＣ_εに関するこの推定値を最大化する。

Ａ．１．Ｅステップ
式（Ａ．１）における結合密度ｐ（θ，ｒ｜Ｃ_ε）は、
と記載することができる。両方の密度がガウスであるため、以下の通りとなる。
式（２．９）および（２．１０）に記載のパラメータη_θ｜ｙおよびＣ_θ｜ｙを有する、θ｜ｙの事後密度もガウスである。θが事後密度ｑを有する場合にｌｎｐ（θ，ｒ｜Ｃ_ε）の平均値を計算するため、Ａがマトリックス、Ｅｚ＝μおよびＣｏｖ（ｚ）＝Σである二次形式ｚ’Ａｚの期待値が以下の通りであることを利用する。
それゆえ、事後分布のもと
であり、
に続くものである。類似するように、
は、以下の式を与えるものである。
それゆえ、ＭステップにおいてＣに対して最大化される式は、以下の通りである。

Ａ．２．Ｍステップ
Ｍステップにおける正確な最適化を導くため、最初に、Ｑが共分散の構造（例えば、重みづけ）を説明するＣ＝σ^２Ｑとして、（Ａ．５）の偏導関数を計算する。
である。
である場合、Ｍステップにおける式の一次導関数は、
と記載できる。単一の分散成分、すなわちＣ＝σ^２Ｉのみが存在する場合を考慮する。この場合、式（Ａ．６）は、Ｎが観察の回数である、
にまで減少する。以下の式を生む微分をゼロと等しく設定することで、最適値を直接得ることができる。
アペンディックスＡの終了

図７Ａは、造影剤漏出の異なるレベルを持つファントムに対し、漏出期間を含んだ、本発明にかかるアルゴリズムを適用した結果生まれたＭＴＴマップを示す。同様の推定マップが右側の行に提示される一方、左側の行において、ＭＴＴパラメータ（上側パネル）および漏出パラメータＫ２（下側パネル）に関する正確なマップが提示される。漏出パラメータがヘルツ（秒の逆数）であるのに対し、ＭＴＴマップの単位は秒である。図７Ｂは、カラーの図７Ａに対応するグレースケール版である。

漏出は、血液脳関門の破損、穿刺または穿孔、または内皮の堅い接合状態の損失に関連する。また、血液脳関門が存在しない、中枢神経系の外部組織の動態にも関する。

図９は、本発明にかかる方法またはコンピュータプログラム製品の詳細の概要を表す概略フローチャートである。哺乳動物２００における灌流指数ＰＩを推定する方法は、以下のステップを有している。
Ｓ１動脈入力関数ＡＩＦを推定する。
Ｓ２注入された造影剤の時間の関数として造影剤濃度を表す、取得されたデータＤＡＴおよびＤＡＴ’に基づく時間の関数としての毛細管通過時間分布にパラメトリックモデルＰＭを適用する。
Ｓ３期待値最大化（ＥＭ）型法の期待ステップを適用する。
Ｓ４造影剤濃度の非線形観測モデルにおける、造影剤の観測誤差の分散（Ｃ_ε）の適切な解析表現とともに期待値最大化（ＥＭ）型法の最大化ステップを適用する。
Ｓ５例えば、ε_ｔｈｒ等の収束条件と比較して、ＥＭ型法が収束であるかどうか判定し、もしＥＭ型法が収束ではない場合、ステップＳ３およびＳ４を再び行い、もしステップＳ３およびＳ４のＥＭ型法が、収束条件に基づいて、収束である場合、
Ｓ６灌流指数ＰＩを推定、また任意で表示または記憶する。

要約すると、本発明は、例えば、ヒトといった哺乳動物２００における灌流指数ＰＩを推定する方法に関する。注入された造影剤の時間の関数として造影剤濃度を表すデータ、ＤＡＴおよびＤＡＴ’は、医療用画像システム１００から取得される（図８を参照）。灌流指数ＰＩは、時間の関数としての毛細管通過時間分布のパラメトリックモデルＰＭ、および期待値最大化（ＥＭ）型法といったマイナライズ最大化（ＭＭ）型法を正則化とともに適用することで見出される。マイナライズ最大化型方法は、最大化ステップにおいて用いられる、造影剤濃度の非線形観測モデルにおける造影剤の観測誤差の分散Ｃ_εの適切な解析表現を有する。７名の患者に行われた臨床試験（図６参照）は、特異値分解（ＳＶＤ）および遅延への減少したセンシティビティと比較して改善されたＭＴＴマッピングを示す。

＜参考文献＞
Dempster, A. P. , Laird , N. M. , Rubin, D. B. 著,1977年.Maximum likelihood from in-complete data via the EM algorithm. J. Royal Stat. Soc. B39(1), 1-38ページ
Jespersen, S. N., Ostergaard, L.著, 2012年.The roles of cerebral blood flow, capillary transit time heterogeneity, and oxygen tension in brain oxygenation and metabolism. J Cereb Blood Flow Metab 32, 264-277ページ
Kenneth Lange著, Numerical Analysis for Statisticians, 2010年, Springer．
Mouridsen, K., Friston, K., Hjort, N., Gyldensted, L., Ostergaard, L., Kiebel, S.著, 2006年. Bayesian estimation of cerebral perfusion using a physiological model of microvasculature. NeuroImage 33(2), 570-579ページ．
Ostergaard, L., Aamand, R., Gutierrez-Jimenez, E., Ho, Y.-L., Blicher, J. U., Madsen, S. M., Nagenthiraja, K., Dalby, R. B., Drasbek, K. R., Moller, A., Braendgaard, H., Mouridsen, K., Jespersen, S. N., Jensen, M. S., West, M. J.著, 2012年. The capillary dysfunction hypothesis of alzheimer’s disease.Neurobiol Aging Submitted.
Ostergaard, L., Jespersen, S. N., Mouridsen, K., Mikkelsen, I. K., Jonsdottir, K., Tietze, A., Blicher, J. U., Aamand, R., Hjort, N., Iversen, N. K., Cai, C., Hougaard, K. D., Simonsen, C. Z., Von Weitzel-Mudersbach, P., Modrau, B., Nagenthiraja, K., Riisgaard Ribe, L., Hansen, M. B., Bekke, S. L., Dahlman, M. G., Puig, J. ,Pedraza, S., Serena, J., Cho, T. H., Siemonsen, S., Thomalla, G., Fiehler, J., Nighoghossian, N., Andersen, G.著, May 2013a. The role of the cerebral capillaries in acute ischemic stroke: the extended penumbra model. J. Cereb. Blood Flow Metab.33(5), 635-648ページ．

上記の参考文献は、その全てを参照により本明細書に援用されるものである。

本発明は、ハードウェア、ソフトウェア、ファームウェアまたはこれらの組み合わせの手段によって実施することができる。本発明またはその特徴のいつくかについても、１以上のデータプロセッサおよび／またはデジタル信号プロセッサで走るソフトウェアとして実施することができる。

本発明の実施形態の別個の構成要素は、単一の単位、複数の単位または分離した機能単位等のいずれかの適切な方法で、物理的、機能的また論理的に実施されてもよい。本発明は、単一の単位で実施されてもよく、または物理面および機能面の両方が、異なる単位およびプロセッサに割り当てられてもよい。

本発明を特定の実施形態に関連して説明してきたが、決して、示される実施例に限定されるものとして解釈されるべきものではない。
本発明の範囲は、添付の請求項一式に照らし合わせて解釈されるものである。特許請求の範囲の文脈において、「備える」または「有する」という用語は、他の構成要素または手段の可能性を排除するものではない。また、「ａ」または「ａｎ」等の記述は、複数を除外すると解釈されるべきものではない。図面に示された要素に対する、請求項における符号の使用も本発明の範囲を限定するものとして解釈されるべきものではない。さらに、異なる請求項に記載される個々の特徴は、有利に結合されてもよく、また異なる請求項において、これらの特徴を記述することは、特徴の組み合わせの可能性や有益性を除外しない。

Claims (18)

1. プロセッサ（２０）によって哺乳動物（２００）における灌流指数（ＰＩ）を推定する方法であって、
   前記方法は、
   注入された造影剤の時間の関数としての造影剤濃度を表すデータ（ＤＡＴ，ＤＡＴ’）を取得することと、
   取得したデータに対して、時間の関数としての毛細管通過時間分布のパラメトリックモデル（ＰＭ）を適用することと、
   最適化によって灌流指数を推定するための正則化を施して、マイナライズ最大化（ＭＭ）型法を適用することと、
   によって灌流指数を推定することと、
   を備え、
   前記マイナライズ最大化型法は、最大化ステップにおいて用いられる前記造影剤濃度の非線形観測モデルにおける造影剤の観測誤差の分散（Ｃ_ε）の適切な解析表現を用いることを含む方法。
2. 前記時間の関数としての毛細管通過時間分布のパラメトリックモデル（ＰＭ）は、少なくとも２つのパラメータを有し、前記パラメトリックモデルは、さらに正またはゼロの通過時間に対して定義される、
   請求項１に記載の方法。
3. 前記時間の関数としての毛細管通過時間分布のパラメトリックモデル（ＰＭ）は、ガンマ分布、Ｆ分布、逆ガウス分布、ガウス分布、半正規分布、レイリー分布、ワイブル分布、負の二項分布、ポアソン分布、ボルツマン分布、またはそれらの任意の組み合わせ、またはそれらから導かれる任意の等価物からなるグループから選択される、
   請求項１または２に記載の方法。
4. 前記毛細管通過時間分布のパラメトリックモデル（ＰＭ）は、Ｒが残留関数で、αおよびβが２つのパラメータである同次多項式：
   で、残留関数Ｒ（ｔ）の推定を可能にするように、時間の関数（Г（ｔ））としてのガンマ関数である、
   請求項１または２に記載の方法。
5. 経時的な前記造影剤濃度の非線形観測モデル（ｙ_ｔ）は、第一要素（ｆ）および第二要素（ε）を有し、初項は、血流と、残留関数と畳み込まれる組織濃度（Ｃ_ａ）と、の積を表し、第二項は、その観測誤差（ε_ｔ）を表している、
   請求項１から４のいずれか一項に記載の方法。
6. 経時的な前記造影剤濃度の非線形観測モデル（ｙ_ｔ）は、脈管構造から脈管外および細胞外空間への経時的な造影剤の漏出要素を表す第三項をさらに備える、
   請求項５に記載の方法。
7. 前記非線形観測モデルは、θが、前記血流、前記パラメトリックモデルのパラメータ（α，β）、および前記Ｃ_ａの測定部位と前記ｙ_ｔの測定部位との間の遅延を表す時間遅延（δ）の関数である未知パラメータである、同次多項式：
   を有する、
   請求項５に記載の方法。
8. 前記非線形観測モデルの線形化は、
   である場合に、
   の通りに行われる、請求項７に記載の方法。
9. ｙが造影剤濃度曲線を表し、前記ｙの分布が統計的線形モデルであるように、
   である場合、残差ｒは、
   と定義される、
   請求項８に記載の方法。
10. 前記統計的線形モデルの１以上のパラメータは、尤度関数または事後確率等の灌流に関する関数を最適化するため、前記マイナライズ最大化（ＭＭ）型法を用いて推定される、
    請求項９に記載の方法。
11. 前記関数は前記事後確率であって、前記マイナライズ最大化（ＭＭ）型法の収束が、前記統計的線形モデルのパラメータの事前の分布の選択を評価するために監視される、
    請求項１０に記載の方法。
12. マイノリゼーションは、以下の関数の期待値
    を用いて取得される、
    請求項１０に記載の方法。
13. 前記造影剤の観測誤差の分散（Ｃ_ε）は、少なくともＪ、Ｃ_θ｜ｙおよびに前記残差ｒに関数的に依存する、
    請求項１０または１２に記載の方法。
14. 前記造影剤の観測誤差の分散（Ｃ_ε）は、Ｃ_ε＝σ^２Ｉという仮定のもと、Ｎが前記造影剤濃度測定値の数である場合に、
    と記述することができる、
    請求項１３に記載の方法。
15. 前記方法は、前記哺乳動物の毛細血管床における血流の不均一性を表す灌流指数をさらに推定する、
    請求項１から１４のいずれか一項に記載の方法。
16. 前記方法は、前記哺乳動物の毛細血管床における血液からの、物質の抽出能力（ＥＣ）を表す灌流指数をさらに推定する、
    請求項１から１５のいずれか一項に記載の方法。
17. 関連する哺乳動物における灌流指数を推定する医療用画像システム（１００）であって、
    注入された造影剤の造影剤動態を表すデータ（ＤＡＴ，ＤＡＴ’）を取得する画像モダリティ（１０）であって、任意で前記データを記憶する記憶手段を有する画像モダリティ（１０）と、
    請求項１から１６のいずれか一項にかかる灌流指数を推定するよう準備されたプロセッサ（２０）と、を備えるシステム。
18. 請求項１から１６のいずれか一項にかかる方法を適用することで、接続されたデータ記憶手段を有する少なくとも１つのコンピュータを備えるコンピュータシステムに、請求項１７にかかる関連する哺乳動物における灌流指数を推定する画像モダリティおよび／または医療用画像システムを制御させる、コンピュータプログラム製品。