JP3819032B2 - 組織およびその他のランダム媒体における蛍光寿命に基づく撮像および分光分析 - Google Patents

Description

発明の背景
本発明は、異質光散乱媒体(heterogeneous light scattering media)の分光撮像に関し、更に特定すれば、生物学的組織の生体撮像に関するものであり、時変(time-varying)光源からの励起光に応答して放出される光の検出による、当該組織の蛍光特性をマッピングすることによって得るものである。なお、これのみに限定される訳ではない。
病気の早期発見は、治療介入(therapeutic intervention)の有効性向上を約束する。近年、非破壊的な技術が開発され、患者の組織における生化学的変化を検出することによって、種々の病気(affliction)を信頼性高くしかも早期に診断する機能が向上した。例えば、磁気共鳴撮像（ＭＲＩ：Magnetic Resonance Imaging）は、常磁性核のスピン状態の緩和を監視することによって、組織の生物医学的撮像および生化学的分光分析を提供することに成功している。しかしながら、ＭＲＩ診断は複雑であり費用がかかることから、特に、病気に対する監視手段として、その適用用途が限定されている。
生物科学への応用が増大しつつある他の強力な分析技術に、蛍光分光分析(fluorescence spectroscopy)がある。その用途には、生物医学的診断、遺伝シーケンス(genetic sequencing)、および流動血球計算法(flow cytometry)が含まれる。今日までに、Ｃａ++，ｐＨ，グルコース(glucose)，ｐＯ₂，およびｐＣＯ₂のような、患者の新陳代謝および環境状態を観察するための蛍光性および燐光性化合物を開発した工業的および学術的研究所がいくつかある。近赤外線（ＮＩＲ：near-infrared red）波長レジームにおいて励起および再放出する染料および光動作的蛍光剤(photodynamic fluorescent agent)の開発により、組織内の深い位置にある病気に冒された組織の非破壊的検出も可能となった。何故なら、赤色光の励起および再放出光は、組織−空気界面間のかなりの距離の伝搬を可能とするからである（Wilson et al.の"Time-Dependent Optical Spectroscopy and Imaging for Biomedial Applications"(80 Proceedings IEEEE,第９１８−３０頁（１９９２年）)を参照のこと）。
Richards-Kortum et al.の米国特許第５，４２１，３３７号、およびWu et al.の米国特許第５，４５２，７２３号によって例示されるように、幾人かの研究者が、非破壊的外部測定または破壊を最小限に抑えた内視鏡測定技術(endoscopic measuring techniques)を用い、病気等に冒された組織および正常組織の区別化を蛍光の放出に基づいて行う種々の手法を提案している。しかしながら、これらの手法は、有効な空間撮像手順を提供するには至っていない。蛍光に基づく撮像が十分に実用化されない理由の１つとして、組織のようなランダムに多重発散する媒体から有意な相応の蛍光特性測定値を得ることが困難であることが上げられる。例えば、蛍光性化合物（即ち、蛍光体）の濃度即ち「吸収量」の関数である蛍光強度(fluorescent intensity)は、撮像に可能な１つの候補である。しかしながら、微粒子（セル）懸濁液、粉体、または組織のような、光学的密度が高い媒体においてこの特性を用いた場合、局部的な散乱および吸収特性が、測定された蛍光強度に影響を及ぼす。
強度以外に、蛍光量子効率(fluorescent quantum efficiency)や寿命のような、選択された蛍光体の他の特性にも、局部的な生化学的環境に感応するものがある。ここで用いる場合、「蛍光量子効率」とは、吸収された核励起光子に対して再放出された蛍光光子の割合（分数）、または蛍光光子の放出が結果的に得られた崩壊事象(decay event)の率を意味するものとする。また、ここで用いる場合、「蛍光寿命(fluoresent lifetime)」は、活性化された蛍光体の平均残存時間、または励起光子の吸収および蛍光光子の再放出間の平均時間として定義する。強度と同様、他の蛍光特性の測定も、多くの場合、研究所における十分に定義された生体への応用、または散乱、吸収および変化する蛍光濃度のような発生事象(issue)の制御または測定が可能な、流動血球計算法(flow cytometry)に限定される。更に、これらの限定のために、腫瘍や目視検査では検出不可能なその他の冒された組織領域のような、隠れた異質組織の意味のある蛍光に基づく撮像も、一般に妨げられている。
したがって、組織の内在的な光学特性に関する情報をそれほど必要とせずに、蛍光降伏(fluorescence yield)や寿命特性によって与えられる濃淡表現機能(contrast capability)を利用し、１つ以上の蛍光特性に基づいて、多重散乱組織を非破壊的に撮像し、異質組織の識別に役立てることが、引き続き求められている。本発明はこの要求を満たすものである。
発明の概要
本発明は、異質な光散乱物質の分光撮像に関するものである。本発明の態様は、新規であり、進歩性があり、種々の利点を提供する。本発明によって包含される実際の特徴は、添付の請求の範囲を参照することによってのみ定義され得るものであるが、本発明の特徴となるいくつかの特徴を端的に述べると次の通りである。
本発明の特徴の１つは、異質光散乱物質を撮像する技術である。このプロセスは、物質の表面を光源からの光に露出させ、それに応じて放出される光を検出することを含む。この物質の蛍光特性の空間変動が、プロセッサによって、放出の関数として定義される。空間変動は、位置の関数として蛍光特性を表す１組の値によって特徴付けることができる。画像の生成は、物質の異質組成に対応する空間変動にしたがって行われる。この技術は、外部または内視鏡的な機器を用いて、生体内の生物学的組織に適用し、病気を示す異質性を検出することができる。この技術は、物質内に蛍光剤を導入することを含んでもよい。検出される蛍光特性は、蛍光寿命、蛍光量子効率、蛍光体吸収係数、蛍光降伏（蛍光量子効率および蛍光体吸収度の関数）、または、当業者には公知のその他の蛍光特性とすることができる。
本発明の他の特徴では、異質組成を有する光散乱物質の空間変動の判定において、光学特性または蛍光特性の変動推定値を確定し、物質からの計算上の放出を、前記推定値の関数として判定し、計算上の放出を放出検出値と比較し、対応する誤差を判定する。変動推定値を修正し、この修正した推定値を用いて計算上の放出を判定し直し、誤差が所望の最小値に到達するまで、比較を繰り返す。物体の画像は、異質組成に対応する修正推定値から生成する。
したがって、本発明の目的の１つは、物質の異質組成と共に変化する光散乱物質の蛍光特性をマッピングすることにより、対応する画像を生成することである。
本発明のその他の目的は、生体器官内に隠れている組織部分(tissue volume)の蛍光特性を非破壊的に監視し、生体中の器官の選択された代謝産物(metabolite)を監視する分光分析技術を提供することである。
更に他の目的は、内在性または外来性蛍光体を造影剤(contrast agent)として用い、冒された組織を識別する蛍光撮像システムおよび方法を提供することである。この撮像は、濃度、寿命、または正常組織および冒された組織間の量子効率差によって得ることができる。
本発明の更に他の目的は、局部的な蛍光体濃度には依存しない光学特性のコントラストに基づく、撮像技術およびアルゴリズムを提供することである。
本発明の更に別の目的、特徴、対応、および利点は、図面および発明の詳細な説明から明らかとなろう。
【図面の簡単な説明】
図１は、本発明の一実施例のシステムの構成図である。
図２は、図１のシステムにおいて実行される方法のフロー・チャートである。
図３は、本発明の種々の態様を例証するために用いる組織模型構成を表す図である。
図４〜図７は、本発明に用いられる数式の選択された特性を表すグラフである。
図８および図９は、本発明の一実施例を利用し、蛍光降伏および残存時間の空間変動判定シミュレーションにおける収束を示すグラフである。
図１０〜図１４は、本発明の実験例１〜３から得られた画像である。
図１５は、本発明の別の実施例のシステムの構成図である。
好適実施例の説明
本発明の原理の理解を促進する目的のために、図面に示す実施例をこれより参照し、具体的な用語を用いてこれを説明する。しかしながら、これによって本発明の範囲を限定する意図はないことは理解されよう。ここに記載する装置におけるあらゆる改造やその他の変更、および個々に記載する本発明の原理のあらゆる別の応用も、本発明に関連する当業者には通常に想起され得るものであろう。
図１は、組織１００の蛍光撮像のための、本発明のシステム１１０を示す。組織１００は、表面１０１を有し、かつ表面１０１の下に位置する領域１０２，１０３で表される異質組成を有する。異質部分１０２，１０３は、通常、表面１０１の目視検査では検出不可能である。
システム１１０は、所定の周波数および波長の強度変調励起光を、光ファイバ１２３を通じて組織１００に供給する変調光源１２０を含む。好ましくは、該光源１２０は、変調出力が１〜５００ＭＨｚ周波数範囲であり、白黒出力が１００〜１０００ｎｍ波長範囲の、従来からの構造のレーザ・ダイオードである。組織１００内において指定された特定の蛍光体を励起するために、特定の波長が選択される。ビーム・スプリッタ１２６を用いて、励起信号の一部分を、処理の目的のために基準センサ１２８に差し向けてもよい。
また、システム１１０は検出サブシステム１４０も含み、この検出サブシステム１４０は、組織１００の多数の検出部位から放出された光子を検出するための光ファイバ１４３を有する。サブシステム１４０は、１つ以上の放出（放射光）センサ１４８を含む。また、検出サブシステム１４０はインターフェース・フィルタも含み、これにより、組織１００内の指定された蛍光体の放出光に対応する所定の放出光波長を得るようにしている。一実施例では、サブシステム１４０は、単一のセンサ１４８を含み、ファイバ群１４３からの信号を多重化して検出する。好ましくは、センサ１２８，１４８は、光増倍管（ＰＭＴ：Photo-multiplier Tube）、またはフォトダイオードであるが、画像増強器(image itensifier)および電荷結合素子のような、その他の様々なセンサも考えられる。
センサ１２８，１４８および光源１２０は、ヘテロダイン・サブシステム１３０に動作的に結合されている。サブシステム１３０は、従来からのレーザ・ヘテロダイン技術を用いて、センサ１２８によって検出された光の、センサ１４８によって検出された光に対する位相、ＡＣ（交流）、ＤＣ（直流）強度に関する情報を得るように構成されている。一実施例では、ヘテロダイン・サブシステム１３０は、光源１２０に用いるレーザの繰り返し周波数（反復率：repetition rate）に位相ロックされた信号シンセサイザを含む。この実施例では、サブシステム１３０は、増幅器を含み、該増幅器は、レーザ反復率（パルス状レーザを用いる場合）の高調波、または変調周波数（変調レーザ・ダイオードを用いる場合）にオフセットを加えたものにおいて、センサ１２８，１４８の利得を変調し、所望のヘテロダインを与える。本実施例の変形においては、８０ＭＨｚのパルス状レーザ反復率を１０ＭＨｚに分周してシンセサイザに入力し、１００ｋＨｚのヘテロダイン・オフセットをセンサ１２８，１４８用の増幅器に入力する。
センサ１２８，１４８は、動作的にプロセッサ１６０に結合されている。プロセッサ１６０は、入力／制御装置１６２、出力装置１６４、およびメモリ１６６を含む。プロセッサ１６０は、１つ以上の構成要素から成る電子回路とすることができる。同様に、プロセッサ１６０は、デジタル回路、アナログ回路、またはそれらの双方から成るものとすることも可能である。また、プロセッサ１６０は、プログラム可能とすることも、一体化された状態マシンとすることも、あるいはそれらの混成体とすることも可能である。好ましくは、入力装置１６２はキーボード等の従来からの様々な入力制御装置であり、出力装置１６６は、ビデオ・ディスプレイである陰極線管（ＣＲＴ：Cathode Ray Tube）、プリンタ、または当業者には公知のその他の画像表示システムである。メモリ１６６は、好ましくは、種々の電子的（例えば、固体）、磁気的、または光学的なもので、電子コントローラまたはプロセッサとの使用のために容易に入手可能な形式のものである。更に、メモリ１６６は、光ディスク・メモリ（ＣＤ）、電磁ハード・ディスク媒体またはフロッピ・ディスク媒体、あるいはそれらの組み合わせを含むことも可能である。
図２は、システム１１０の一つの動作態様を、プロセス２１０として示している。プロセス２１０は、プロセッサ１６０を用いて蛍光降伏および残存時間の空間変動をマッピングし、このマッピングにしたがって画像信号を発生することを含む。出力装置１６４は、この画像信号に応答して画像を表示するように構成されている。プロセス２１０は、ステップ２１２において、蛍光造影剤を組織１００に導入することにより開始される。この造影剤は、サブシステム２４０による検出のための、蛍光放出源を提供する。変調光源１２０、ヘテロダイン・サブシステム１３０、および検出サブシステム１４０の構成は、選択された蛍光剤の励起および放出特性に適応するように設計される。他の実施例では、内在性蛍光体(endogenous fluorophore)を代わりにまたは追加的に用いて、それに応じてシステム１１０を調整するようにしてもよい。
ステップ２１４において、組織１００は、選択された蛍光体に応じて構成された光源１２０によって励起される。ステップ２１６において、光源１２０からの励起光に対する、各検出部位「ｉ」における放出光の位相θ_obsおよびＡＣ強度の対数(log)Ｍ_obsを、ヘテロダイン（またはオフセット）周波数において判定する。検出部位「Ｄｉ」箇所に対して、検出即ち観察した位相およびＡＣ強度に、それぞれ（θ_obs）_iおよび（Ｍ_obs）_iを用い、「ｉ」をインデックスとする。プロセッサ１６０は、相対位相およびＡＣ強度の情報をメモリ１６６に格納する。
ステップ２１８において、撮像に選択した組織１００の領域に対して二次元格子を確定し、格子点のマトリクスを確定し、「ｊ」を格子点のインデックスとする。蛍光降伏に対する均一シード値(uniform seed value)、すなわちｙ_j＝（ημ_ax→m）_j、および蛍光寿命（τ）_jを、各格子点ｊに指定する。これらの値は、後のステップにおいて変更される、降伏および寿命値の初期均一推測値である。「η」項は、蛍光体の量子効率であり、蛍光体の周囲の環境によって変動する。「μ_ax→m」項は、蛍光体に対する吸収係数であり、自然対数に基づく蛍光体の吸収光度(extinction)と蛍光体の濃度の積である。その結果、降伏ｙ＝（ημ_ax→m）は、周囲の新陳代謝および蛍光体の吸収量による影響を受ける。ある公知の蛍光体の吸収量は、宿主（ホスト）組織(host tissue)の種類や状態によって異なり、病気の検出に有用な別の蛍光特性を与える。これらの特性によって提供されるコントラストは、蛍光体濃度には概ね無関係である。蛍光降伏および寿命の初期推定値は、プロセッサ１６０によって、後の使用のためにメモリ１６６に格納される。
降伏ημ_ax→m、および寿命τの蛍光特性の初期推定値をこのように確定した後、処理ループ２２０はステップ２３０に入る。好ましくは、処理ループ２２０のステップは、プログラムされたソフトウエア、専用ハードウエア、または双方の組み合わせの適切なものを用い、プロセッサ１６０が実行する。プロセス２１０およびループ２２０の多様な数学的態様の理解に役立てるために、選択した変数の表を以下に纏めておく。
ｃ 光速
Ｄ（ｒ） 光学的拡散係数
Ｄｉ 検出部位数
ｆ 変調周波数
Ｉ 識別マトリクス
ｉ 検出部位インデックス
Ｊ 各格子点ｊにおける感度を各検出部位における応答に関係付けるヤコビの行列
ｊ 格子点インデックス
Ｊ_j,i ヤコビの行列Ｊの個々の要素
ｋ ソース・インデックス
Ｍ 変調蛍光位置のＡＣ強度の対数
ｍ 多数の変調周波数に対するインデックス
ｎ 平均屈折率
ｒ 位置（二次元または三次元）
Ｓｋ 変調光源数
Ｓ（ｒ，ω） 位置ｒおよび周波数ωにおける変調光に対する原始項
ギリシア文字
χ² 最小二乗誤差を表す効果関数
Φ_x（ｒ，ω） 周波数領域内の位置ｒおよび周波数ωにおける格子束を表す複素数
η 蛍光プローブまたは染料（ダイ）の量子効率
μ_a 平均吸収係数
μ_am 非蛍光発色団(non-fluorescing chromophores)および発光体双方による蛍光の吸収係数
μ_ax 非蛍光発色団および発光体双方による励起光の吸収係数
μ_ax→c 非蛍光発色団による励起光の吸着係数
μ_ax→m 発光体による励起光の吸着係数
μ′_s 有効散乱係数
θ ある変調光波から他への位相シフト
τ 位置ｒにおける活性化されたプローブまたは染料の寿命
ω ２πｆで与えられる角度変調周波数
添字
ｏｂｓ 観察データまたは経験データ
ｘ 励起光
ｍ 蛍光または放出光
ステップ２３０において、各検出部位「ｉ」における位相および相対的ＡＣ強度を、各格子点ｊに対する、降伏および残存時間の初期推定値の関数として計算する。計算した位相および強度は、各検出部位ｉにおいて、それぞれ、（θ_m）_iおよび（Ｍ_m）_iとして表される。（θ_m）_iおよび（Ｍ_m）_iの値は、放射伝達（通過）方程式(radiative transport equation)の近似拡散方程式(diffusion equation approximation)を用いて決定する。近似拡散方程式は、組織または多重散乱媒体の空間的および時間的伝達を記述する。結合された周波数領域の拡散方程式を用いて、式（１）および（２）によって、組織１００の選択された格子内のいずれの位置ｒにおいても、励起および放出流束率(emission fluence rates)Φ_x（ｒ，ω）およびΦ_m（ｒ，ω）をそれぞれ予測することができる。
∇・［Ｄ_x（ｒ）∇Φ_x（ｒ，ω）］
−［μ_ax（ｒ）＋ｉω／ｃ_n］Φ_x（ｒ，ω）＋Ｓ_x（ｒ，ω）＝０ （１）
∇・［Ｄ_m（ｒ）∇Φ_m（ｒ，ω）］
−［μ_am（ｒ）＋ｉω／ｃ_n］Φ_m（ｒ，ω）＋Ｓ_m（ｒ，ω）＝０ （２）
励起光Ｓ_x（ｒ，ω）に対する原始項(source term)は、各周波数ω＝２ωｆにおける正弦変調光によるものであり、ここでｆは一般にＭＨｚ周波数範囲内である。拡散方程式（１）および（２）双方における第１項は、光の拡散的即ち「ランダム・ウオーク」伝達を表す。なお、Ｄ_x,mは、以下の式（３）に示す光拡散係数である。
Ｄ_x,m＝［３（μ_ax,m＋μ'_sx,m）］^-1 （３）
そして、μ_aおよびμ'_sは、組織１００、即ち、対象媒体のそれぞれ吸収および等方的散乱係数である。光学的特性は光の波長に依存するので、そのため、光源１２０（添字ｘ）からの励起光およびサブシステム１４０（添字ｍ）によって検出される蛍光放出では異なる。励起波長μ_axにおける全吸収係数は、励起波長に応答する非蛍光発色団および蛍光体に関連するものである。全吸収係数は、非蛍光発色団μ_ax→cおよび発光体μ_ax→mによる吸収係数の和によって与えられる。通常、蛍光波長において発生する吸収は、主に非蛍光発色団によるものであると仮定することができる。組織内の光速はｃ_n＝ｃ／ｎであり、ここでｎは平均屈折率である。蛍光放出に対する原始項は、励起光流束Φ_x（ｒ，ω）に依存し、以下の式（４）で与えられる。
Ｓ_m（ｒ，ω）
＝ημ_ax→m（ｒ）Φ_x（ｒ，ω）［（１−ｉωτ（ｒ））
／（１＋ω²τ（ｒ）²）］ （４）
この項は、励起光の入射パルスに続く、時間領域における蛍光崩壊項のフーリエ変換から得られ、τは蛍光体の寿命、ηは量子効率、および吸収係数μ_ax→mは自然対数に基づく励起係数と基底状態にある発光体の濃度との積である。先に示したように、組み合わせた積μ_ax→mは蛍光降伏ｙであり、生成される蛍光流束に比例する。式（４）を式（２）に代入することによって、各格子点「ｊ」に対するΦ_mの決定が容易になる。格子点「ｊ」によって定義される二次元領域に対する拡散方程式（１）および（２）の解は、容易に三次元に拡張可能であり、三次元における位置に対応する「ｒ」を有する選択された三次元部分(volume)における１つ以上の蛍光特性の空間変動を推定することができる。
拡散方程式（１）および（２）は双方共、線形複素楕円方程式であり、複素数Φ_x（ｒ，ω）およびΦ_m（ｒ，ω）についての境界値問題として解くことができる。この解は、有限差(finit difference)の方法を用い、対応する有限差方程式を作成する。これら有限差方程式を利用して、各格子点ｊにおける近似解を得る。この解法は、Fulton et al.のMultigrid Method for Elliptic Problems.A Review,114 American Meteorological Society pp.943-59（１９８６年５月）、およびB.W.Pogue et al.のInitial Assessment of a Simple System for Frequency Domain Diffuse Optical Tomography,40 Physics in Medicine and Biology pp.1709-1729（１９９５年）において、記載されている。この解決法を実行する好適な方法の１つは、Adams,J.C.のMUDPACK:Multigrid Portable Fortran Software for the Efficient Solution of Linear Elliptic Partial Differential Equations,34 App.Math Comp.p.133（１９８９年）に記載されているMUDPACKルーチンを使用することである。有限差方程式の解については、組織１００の表面１０１において、Φ_m,x（ｒ，ω）＝０と仮定する。これは、ゼロ流束境界条件として知られている。他の境界条件も選択可能であり、解法もそれに応じて変化することは認められよう。
拡散方程式（１）および（２）は、各格子点ｊにおけるΦ_mの複素数について解くことができる。表面で検出される信号は光子流束の傾斜の垂直成分(normal component)に比例する。組織１００の表面１０１に位置する検出器部位「ｉ」における信号を近似するために、当該部位に最も近い内部格子点におけるΦ_m値を選択する。これは、光子流束傾斜の垂直成分は、表面１０１の直ぐ内側のΦ_mに比例するという関係によるものである。検出部位における計算した位相ラグθ_m、およびＡＣ強度の対数Ｍ_mを、光源１２０の位相およびＡＣ強度に関して、複素数Φ_mの虚部および実部から計算する。
拡散方程式（１）および（２）は、検出器部位ｉにおいて測定したθ_mおよびＭ_mについて、組織１００の蛍光光学特性を変更する感度に対する洞察を与える。この洞察は、拡散方程式（１）および（２）の種々のパラメータを固定した一連の計算から得られる。これらの計算は、図３に示すように、模型の背景３０３に隠されている包埋異質部分３０２を有する円形の組織模型３００を想定している。模型３００に対して二次元格子を確定するが、三次元に容易に拡張することができる。これらのシミュレーション条件の下で、シミュレートされた組織模型の外側の全ての格子点における励起および蛍光双方に対する吸収係数に、大きな値を指定する。図３の４つの光源（Ｓｋ＝４）は、各光源に最も近い表面付近の格子点において任意の複素数を指定することによってシミュレートする。検出部位に最も近い格子点「ｊ」におけるΦ_mから決定された計算値を用いて、図３の２０カ所の検出部位Ｄ１〜Ｄ２０（Ｄｉ＝２０）をシミュレートする。拡散方程式（１）および（２）に対するシミュレーションの解は、６５ｘ６５格子に対して、二次元で得られた。この格子は、組織模型３００の中央に位置する直径２０ｍｍの円形包埋異質部分（この位置は、図３の異質部分３０２の構成とは多少異なる）を有する、直径１００ｍｍの円形組織模型３００を覆う。蛍光位相シフトおよびＡＣ強度のシミュレーション測定値は、等間隔で円周上に位置する２０カ所の検出部位Ｄ１〜Ｄ２０に対して報告される。変調周波数ｆは１５０ＭＨｚに等しく設定した。異質部分および背景の光学特性を、以下の表１に示す。

ημ_ax→mの影響を評価するために、異質部分内のημ_ax→mの値を１０^-4ｍm^-1から１０^-1ｍm^-1まで増加させ、背景３０３におけるημ_ax→mの値を一定に維持して、各検出部位Ｄ１〜Ｄ２０においてθ_mおよびＭ_mを計算した。寿命τは、物体および背景ともに１ｎｓに等しく設定し、ημ_ax→mにおける差によるコントラストを発生させた。１つのアクティブな光源Ｓ１について、θ_mおよびＭ_mのプロットを図４および図５にそれぞれ示す。異質部分１０２のημ_ax→mが大きな値に増加するに連れて、ＡＣ強度は、希釈非散乱溶液(dilute non-scattering solution)において予測されるものと同様、上限に近づいていく。図５は、蛍光体μ_ax→mによる吸収係数が背景の１０〜１００倍に減少する際に、蛍光位相シフトθ_mがどのように減少するのかを示す。これらのシミュレーションから、Ｍ_mはシミュレートした組織異質部分１０２のημ_ax→mにおける変化に直接依存するように見え、これに対してθ_mは、光子の移動の変化によって、ημ_ax→mに間接的に依存する。
τの影響を評価するために、異質部分におけるτの値を１０^-1ｎｓから１０³ｎｓまで変化させ、背景におけるτの値を１ｎｓに保持して、検出部位Ｄ１〜Ｄ２０の各々におけるθ_mおよびＭ_mを計算した。背景のημ_ax→mは１０^-5ｍm^-1に設定し、異質部分に対するημ_ax→mは１０^-3ｍm^-1に設定した。図６に示すように、検出したＡＣ強度は、τが減少するに連れて増加している。図７は、異質部分の寿命を０．１ｎｓ〜１０００ｎｓに変化させた場合の、各検出部位における蛍光位相シフトの値を示す。所定の変調周波数（この計算では１５０ＭＨｚ）において、τを０．１ｎｓから１０００ｎｓに増大させると、θ_mは最初に減少し、極大値に到達し、その後減少する。したがって、各検出部位Ｄ１〜Ｄ２０におけるθ_mおよびＭ_mは、異質部分における寿命の値に直接影響を受けるように見える。
図２に戻り、ステップ２４０において、各検出部位「ｉ」について計算した放出光位相および強度（θ_m）_iおよび（Ｍ_m）_iを、測定した放出光位相および強度（θ_obs）_iおよび（Ｍ_obs）_iと比較し、測定値および計算値間の差即ち「誤差」を確認する。（ημ_ax→m）_jは（Ｍ_m）_iに影響を与えるので、この比較は、以下の式（５）の効果関数χμ²の形態で示される。

なお、σＭは、Ｍ_mにおけるノイズの典型的な標準偏差であり、０．０１とする。Ｓｋはｋをインデックスとする励起源部位の数であり、Ｄｉはｉをインデックスとする検出部位の数である。このアルゴリズムの目的は、（ημ_ax→m）_jの適切な更新によって、χμ²を最小化することである。（ημ_ax→m）_jの最初の更新の後、（τ）_jに関する他の効果関数がステップ２４０の比較に関与する。この効果関数χτ²は、次の式（６）のように表される。

ここで、σ_θは（θ_m）_iにおけるノイズの典型的な標準偏差であり、１度とする。Ｓｋはｋをインデックスとする励起源部位の数であり、Ｄｉはｉをインデックスとする検出部位の数である。寿命は（θ_m）_iおよび（Ｍ_m）_i双方に影響を与えるので、位相およびＡＣ強度が式（６）に用いられる。
効果関数χμ²およびχτ²を計算することによってステップ２４０の比較を行った後、制御は条件ステップ２５０に進み、効果関数による観察値（θ_obs）_iおよび（Ｍ_obs）_iの計算値（θ_m）_iおよび（Ｍ_m）_iに対する比較が、選択された収束基準を満たすか否かについて検査を行う。この基準は、降伏値および寿命値を判定する際の許容誤差の度合いに対応する。一実施例では、以下の３つの量（ｉ）χ²、（ｉｉ）ループ２２０の連続する繰り返しにおけるχ²の変化、（ｉｉｉ）ループ２２０の連続する繰り返しにおけるχ²の相対的変化のいずれかが、所定のスレシホルド値１．０ｘ１０^-2未満となった場合に、収束が達成される。他の実施例では、当業者に想定されるような、別の比較計算および関連する条件を用いることも可能である。ステップ２５０の条件が満たされた場合、制御はステップ２７０に進み、ループ２２０から抜け出す。しかしながら、基準が満たされない場合、ステップ２６０においてループ２２０の実行を継続する。
ステップ２６０において、降伏（ｙ）_j＝（ημ_ax→m）_jおよび寿命（τ）_jを、各格子点ｊについて更新し、これらの値が、比較ステップ２４０および条件２５０の検査に対応する最小誤差に到達するようにする。これらの値を更新するために、ヤコビの行列を用い、各検出位置ｉにおける応答の（ｙ）_j＝（ημ_ax→m）_jの変化に対する感度、および各格子点ｊにおける寿命（τ）_jを記述する。

という３つのヤコビの行列を用いる。これらヤコビの行列の要素Ｊ_i,jは、それぞれ、

で与えられる。これらの要素は、各格子点ｊについて拡散（１）および（２）を４回解くことによって計算することができ、（τ）_jおよび（τ＋δτ）_j、ならびに（ημ_ax→m）_jおよび（ημ_ax→m＋δημ_ax→m）_jによって計算されるＭ_m，ｉおよびθ_m，ｉを得る。最小二乗最小化から、降伏および寿命に対する更新を計算する。一好適実施例では、この更新アルゴリズムは、Yorkey,et al.,Comparing reconstruction Algorithms for Electrical Impdedance Tomography,34 Transactions in Biomedical Engineering pp.843-52（１９８７年）によって提案されたアルゴリズムと同様の、電気インピーダンス断層撮影法によって得られる画像を再構成するために用いられるアルゴリズムから改造したものである。ヤコビの行列を用いて更新ベクトル

を解き、降伏および寿命ベクトル

をそれぞれ推定する。これらのベクトルは、格子点の数に対応する次元のものである。ループ２２０を通過する各繰り返し毎に、以下のヤコビの方程式（７）および（８）を解き、推定した降伏および寿命ベクトルに対する更新を決定する。

は、それぞれ、ｉ箇所の検出部位各々における、ＡＣ強度の対数の観察ベクトルおよび計算ベクトルである。

は、それぞれ、ｉ箇所の検出部位各々における、位相ラグの観察ベクトルおよび計算ベクトルである。ヤコビの行列の不良条件的性質のために、項λ₁Ｉまたはλ₂Ｉがマルカルト最小化方式の一部として加えられており、Ｉは識別マトリクスである。パラメータλ₁またはλ₂は、Press et al.,Numerical Recipes:The Art of scientific Computing,(Cambridge University Press,1992)に開示されている形式のマルカルト−レーベンバーグ型アルゴリズム(Marquardt-Levenberg type algorithm)によって調整する。従来の数値方法を用いて、ヤコビの行列方程式（７）および（８）から得られる連立線形代数方程式を解く。ループ２２０を通過する各繰り返しにおいて、ヤコビの行列を再度計算する。式（７）および（８）は、降伏および寿命の推定値に対する適切な変更を選択する方法を提供することがわかっている。しかしながら、当業者に想起される他の周知手法によって繰り返し計算し、容認可能な推定値を得ることも考えられる。一旦更新が完了したなら、制御はステップ２３０に戻る。
収束基準が条件ステップ２５０において満たされた場合、格子点に対する降伏および寿命の推定は、容認可能な最小値に到達したことになり、制御はステップ２７０に進む。ステップ２７０において、降伏および／または寿命蛍光特性の空間変動から、プロセッサ１６０が画像信号を発生する。この画像信号は出力装置１６４に送られ、それに応答して画像を表示する。降伏および寿命の蛍光特性は、通常、蛍光体の生物的環境と共に変化するので、この画像は、概ね組織変化を示し、異質部分１０２，１０３を検出する機能を提供する。例えば、組織に数センチメートル貫入可能な近赤外線（ＮＩＲ）光を供給可能なレーザ・ダイオード、およびＮＩＲ光に応答する蛍光造影剤を用いれば、実行可能な撮像システムを提供することができる。一実施例では、このシステムは、内視鏡と共に用いるように改造する。
降伏および寿命以外にも、拡散方程式（１）および（２）を用いて、冒された組織を区別するのに有用な他の蛍光特性の空間変動をマッピングすることが可能である。このような代わりの蛍光特性には、降伏の効果(yield product)には依存しない別個の特性として決定される、量子効率ηおよび蛍光吸収係数μ_ax→mの少なくとも一方が含まれるが、これらに限定される訳ではない。
本発明の他の実施例では、光子流束方程式およびヤコビの推定プロセスを適応化し、指定された蛍光体の吸収濃度のマッピングを判定する。この実施例では、降伏および寿命の推定値の代わりに、各光子点における発行団吸着係数μ_ax→cおよび散乱係数μ'_sを推定することによって、指定された発光体がない場合の発光団吸着係数μ_ax→cの第１のマップおよび散乱係数μ'_sを判定する。Φ_x（ｒ，ω）に対する拡散方程式（１）を、変形したヤコビの方程式（７）および（８）と共に用いて、この第１のマップを作成することができる。この変形で、降伏の代わりに、発光団吸着および散乱係数が代入され、これらの新しい特性に適応させる処理(adaptation)後では、次のようになる。

採用した４つのヤコビの行列、

の要素は、それぞれ、

によって与えられる。吸収および散乱マップに対する更新は、効果関数χ２を最小化するように行われる。

ここで、ｎ_s＝Ｓｋおよびｎ_d＝Ｄｉである。
第１のマップを生成した後、指定された蛍光造影剤を導入し、式（９）〜（１１）においてμ_ax→cの代わりにμａｘを代入することによって全吸着係数μ_axを判定し、全吸着係数の第２のマップを得る。尚、μ_ax＝μ_ax→m＋μ_ax→cであり、蛍光造影剤の吸収量は、μ_ax→mに直接比例し、吸収量濃度は、第１および第２のマップに対する吸着係数の変化の差を判定することによってマッピングすることができる。この「差分マップ」を用いて、次に吸収量濃度に対応する画像を生成することができる。
他の代替実施例では、多数の光源変調周波数ｆの各々に応答して、放出を測定する。用いる異なる周波数の全数をＭｆで表す。この追加データを得るために、ループ２２０の繰り返しを、各周波数ｆについてインデックスがｍになるまで実行する。光源の数Ｓｋおよび検出部位Ｄｉは、それぞれ、ｋおよびｉをインデックスとする。この追加データを用いて、システム１１０によって得られる撮像結果の改良や、評価における検出部位数または励起源部位数の減少を図ることができる。この追加データに対応する代表的な効果関数は、以下の式（１２）のように与えられる。

蛍光降伏および寿命の他に、多重周波数方法(multi-frequency method)を用いて、関心のある他の光学特性をマッピングすることも可能である。正弦変調光源の他に、本発明は、パルス状または他の時間によって変化する（時変）励起光源を用いて動作ように改造することも、代替実施例では可能である。
図１５は、本発明の他の実施例の光学システム４１０を示す。このシステムは、レーザ・ドライバ４２２、動作的に結合されているレーザ・ダイオード４２４、および基準周波数発生器４２６を有する変調光源４２０とを含む。光源４２０は、組織模型４００に変調光を送出するように構成され、該模型からの再放出光は、５０ｍｍレンズ４３２を通過して、可変利得の画像増感器(intensifier)４３０上で合焦する。画像増感器４３０は、光子を電子に変換するフォトカソード面(photocathode face)、電子信号をアバランシェ倍増(avalanche multiplication)によって倍増させる多チャネル・プレート（ＭＣＰ：Multi_Channel Plate）、および電子を光学画像に変換する燐光スクリーンを含む。好ましくは、増感器４３０は、Litton Electronics,Inc.,によって製造される製品群の内の高速増感器であり、フォトカソードおよびＭＣＰ間に、増幅器４２８からのＤＣバイアスおよびＲＦ信号を印加することによって、変調を可能にする。例えば、増感器４３０からの画像の変調は、シンセサイザ４２６からの１０ＭＨｚ出力信号によって、レーザ・ダイオード４２４に位相ロックされる。レーザ・ダイオード４２４および画像増感器４３０を同じ周波数で変調することにより、燐光スクリーン上に定常状態画像が得られる。Gratton et al.の米国特許第５，２１３，１０５号は、この技術のある態様に関する別の背景を与える。燐光スクリーンからの画像は、干渉フィルタ４３３を通過し、１５０ｍｍマクロ・レンズ４３６を通過して、電荷結合素子（ＣＣＤ）カメラ４３４上で合焦する。カメラ４３４は、５１２ｘ５１２アレイのＣＣＤ検出器を有し、対応する画素化された画像(pixelated image)を与えるように構成されている。カメラ４３４は、先に述べたプロセッサ１６０と同様の構成のプロセッサ４６０に、動作的に結合されている。
獲得した各画像にしたがって、プロセッサ４６０の制御の下で、周波数シンセサイザ４５２を用いて、画像増感器４３０の位相を０〜２６０度の間の値にステップ的に変化させることによって、画像増感器４３０およびレーザ・ダイオード４２４間に位相遅延を導入する。画像増感器４３０およびレーザ・ダイオード４２４の利得変調は同じ周波数で行われるので、ホモダイン(homodyning)によって、位相に依存する定常的な燐光画像が増感器４３０上に得られる。好ましくは、シンセサイザ４５２およびプロセッサ４６０間の制御は、従来からのＧＰＩＢインターフェースによって得るようにする。画像増感器４３０の燐光スクリーンからの画像は、各位相遅延毎に集められる。次に、位相遅延が徐々に大きくなる画像を用いて、励起光および模型４００からの再放出光間の位相シフトおよび強度変調率のマップを生成する。干渉フィルタまたは適切な光学フィルタを適用することによって、励起光から選択的に放出光を分離し、測定することができる。カメラ４３４の出力は、プロセス２１０を用いて、プロセッサ４６０によって処理することができる。
次に、以下の具体例１〜３を参照しながら、更に本発明を説明する。これらの例は性質上例示的であり、限定的なものではないことは理解されよう。例１〜３は、プロセス２１０のコンピュータ・シミュレーションを含む。この種のシミュレーションは、組織のシミュレーションを含み、当業者に蛍光分光撮像の性能を実証する手段として受け入れることができるものである。これらの例は、以下の表２の条件の下で、θ_mおよびＭ_mについて拡散方程式（１）および（２）を解くことによって得られたシミュレーション値を用いる。

本例は、直径１００ｍｍの図３の組織模型３００のシミュレーションを行う。値θ_mおよびＭ_mは、周囲に位置する４つの変調光源Ｓ１〜Ｓ４に応答する、図３のＤ１〜Ｄ２０検出部位の各々において計算した。励起光変調周波数ｆは１５０ＭＨｚでシミュレーションを行った。拡散方程式（１）および（２）を解いて、検出部位および光源部位の種々の組み合わせに対応する、８０個のθ_mおよびＭ_mのシミュレーション値を得た（Ｓｋ＊Ｄｉ＝４ｘ２０＝８０）。θ_mにおける０．１度（即ち、自由１度）およびＭ_mにおける１％の標準偏差を有するガウス・ノイズを、拡散方程式の解に重ね合わせた。適応化したMUDPACKルーチンを用いて、SunSparc10コンピュータ上で、拡散方程式（１）および（２）を解いた。こうして得られたデータ集合を、例１〜例３に対するプロセス２１０のシミュレーション入力データとして用いた。その結果を、以下の表３および表４に示す。

例１
例１は、非蛍光発色団による、吸収のない場合の蛍光降伏および寿命を再構成する。本例のための実験データをシミュレートするために、背景および異質部分３０２に対する蛍光降伏（ημ_ax→m）_jを、それぞれ、１ｘ１０^-5ｍm^-1および１ｘ１０^-3ｍm^-1、そして背景および異質部分３０２に対する蛍光寿命（τ）_jを、それぞれ、１０ｎｓおよび１ｎｓに選択した。ループ２２０の実行の間、異質部分３０２の位置や背景の蛍光特性等の先見的知識を想定せず、１ｘ１０^-5および１０ｎｓという均一の推測値を、それぞれ、蛍光降伏（ημ_ax→m）_jおよび寿命（τ）_jに与えた。二次元の１７ｘ１７格子では、５０回未満のループ２２０の繰り返し（SunSparc10上での計算時間は２時間）で収束に達した。シミュレートされた物体が占める格子点におけるημ_ax→mおよびτの平均値は、５０回以内の繰り返しで、ημ_ax→m＝０．９３ｘ１０^-3ｍm^-1およびτ＝１．０３ｎｓの場合が、図８および図９にそれぞれ示されている。図１０および図１１は、それぞれημ_ax→mおよびτのマッピングされた値から再構成された画像を示し、予測画像を表す。画像は、例１〜例３における補間によって平滑化が施されており、非物理的に高い値を有していたが、物理的に達成可能な範囲の値に包囲されていた擬似点を除去した。これらの擬似値は、ループ２２０のシミュレーションから得られた、平均背景蛍光降伏および寿命で置き換えた。
シミュレートされた背景が占める格子点におけるημ_ax→mの平均値は、５０回以内の繰り返しで、９ｘ１０^-5ｍm^-1に収束した。背景の値は５．４ｎｓに収束する。最終画像の初期推測値に対する依存性を、（ημ_ax→m）_jおよび寿命（τ）_jに対して１ｘ１０^-4ｍm^-1および１０ｎｓという初期均一推測値を与えることによって検査した。その結果、図１０および図１１において得られたものと同様の画像が得られた。
異質部分３０２の位置は、ημ_ax→mのピーク値の３５％（任意に選択）よりもημ_ax→mが高い格子点全てから成るものとして識別された（図１０）。識別された物体の格子点全ての座標の平均は、位置（６０．８，５８．５）であり、実験データをシミュレートするために用いた位置（６０，６０）に近い値である。表３に纏めたように、識別のための我々の任意の定義に基づく異質部分の面積は、７２ｍｍ２であり、我々のシミュレートした実験データを生成するために用いたものに近い。
例２
例２は、組織を真似るように構成された発色団吸収のシミュレーションによって、蛍光降伏および寿命を再構成する。均一背景発色団吸収係数ημ_ax→に１ｘ１０^-3ｍm^-1を用いて、シミュレートされた実験データを生成したことを除いて、例１において記載したものと同じ隠された異質部分ならびに光学特性およびシミュレーション機器を用いた。励起光伝搬は、画像再構成のためには用いなかったが、この光学特性を考慮した。これは、生理学的条件の下で逆画像再構成のための可能な限り最良の性能を評価することが知られている。蛍光降伏（ημ_ax→m）_jおよび寿命（τ）_jの二次元再構成空間マップを、図１２および図１３にそれぞれ示す。表３に示したように、μ_ax→mに基づく我々の基準による物体の位置の平均値は、実験データをシミュレートするために用いた条件と一致する、位置（５９．４，５８．３）として得られた。識別のための任意の定義に基づく異質部分の寸法（ημ_ax→mがその最大値の３５％よりも高い格子点全て）は７０３ｍｍ²であり、これは、我々がシミュレートした実験データを生成するために用いたものと近い値である。シミュレートされた物体が占める格子点におけるμ_ax→mおよびτの平均値は、５０回以内の繰り返しで、μ_ax→m＝０．８ｘ１０^-3ｍm^-1およびτ＝０．７ｎｓに収束した。これらは、シミュレートした実験データ（表３参照）を生成するために用いた値と一致する。シミュレートされた背景が占める格子点におけるμ_ax→mおよびτの平均値は、５０回以内の繰り返しで、例１について報告したものと同様の値に収束した。
例３
例３は、組織模型（図３には示されていない）に、２つの異質部分が隠された場合をシミュレートしたものである。この場合、物体１および２について、蛍光降伏μ_ax→mを、それぞれ、１ｘ１０^-3ｍm^-1および２ｘ１０^-3ｍm^-1に選択し、異質部分に対する寿命τをそれぞれ１ｎｓおよび２ｎｓに選択したことを除いて、例１に記載したものと同じ光学パラメータを用いた。また、１７ｘ１７格子の代わりに、３３ｘ３３格子を用いた。降伏のマッピングに対応する画像を図１４に示す。
この明細書において引用した全ての刊行物および特許出願書類は、個々の刊行物または特許出願書類各々が、具体的かつ個別的に、その言及によって本願に含まれることが示されたかのように、この言及によって本願にも含まれるものとする。本発明について図面および上述の説明において詳細に示しかつ記載したが、これらは性質上限定的ではなく例示的に解釈すべきものであり、好適実施例のみについて示しかつ記載したこと、および本発明の精神に該当する変更や改良は全て保護対象となることを望むことは理解されよう。

Claims (11)

1. 異質組成を有し蛍光体を含んでいる、光を散乱させる組織を画像化するシステムにおいて、
   （ａ）前記組織の前記蛍光体を励起する励起手段と、
   （ｂ）前記励起手段による励起により前記組織から生じる蛍光放出に対応する光信号を検出する検出手段と、
   （ｃ）前記検出手段に接続され、前記光信号に応答して動作するプロセッサであって、
   前記組織の蛍光特性の空間変動の推定値を確定し、
   前記推定値の関数として計算上の放出を決定し、
   前記計算上の放出を前記光信号と比較して、それらの間の誤差を検出し、
   前記誤差が所定の基準値である最小誤差より大きい場合に、前記組織の前記蛍光特性の空間変動の推定値を更新し、
   前記誤差が前記最小誤差以下となるまで、前記更新された推定値の関数として計算上の放出を決定し、該計算上の放出を前記光信号と比較して、それらの間の誤差を検出し、かつ、前記蛍光特性の空間変動の推定値を更新することを反復的に実行し、前記誤差が前記最小誤差以下となったときに、そのときの前記空間変動の推定値を出力信号として生成する
   手段を含んでいるプロセッサと、
   （ｄ）前記出力信号に応答して、前記組織の前記異質組成に対応する画像を提供する手段と
   からなることを特徴とするシステム。
2. 請求項１記載のシステムにおいて、前記励起手段は、時間的に強度が変動する光を提供する光源を含んでいることを特徴とするシステム。
3. 請求項２記載のシステムにおいて、前記検出手段は、前記組織内の前記蛍光体の励起により該組織から放出される多重散乱光に対応する光信号を提供するセンサを含んでいることを特徴とするシステム。
4. 請求項１記載のシステムにおいて、前記プロセッサは、前記組織の前記異質組成によって変化する蛍光特性のレベルを表す複数の値を、位置の関数として提供し、該複数の値の関数として前記出力信号を生成するよう構成されていることを特徴とするシステム。
5. 請求項１記載のシステムにおいて、前記励起手段は、複数の変調光源を備えていることを特徴とするシステム。
6. 請求項１記載のシステムにおいて、前記蛍光特性は、蛍光寿命、蛍光量子効率、蛍光収量、及び蛍光吸収度の内の少なくとも１つに対応していることを特徴とするシステム。
7. 請求項１記載のシステムにおいて、前記検出手段は、前記組織の表面に沿った複数の位置における前記蛍光放出を検出するよう構成されていることを特徴とするシステム。
8. 請求項４記載のシステムにおいて、前記プロセッサは、前記蛍光特性の推定空間変量の関数として決定される計算上の放出と、前記光信号から得られる測定による放出との比較によって前記複数の値を決定し、前記計算上の放出と前記測定による放出との間の前記誤差が前記所定の最小値以下となるまで、前記空間変動の推定値を更新して比較を繰り返すよう構成されていることを特徴とするシステム。
9. 請求項４記載のシステムにおいて、前記プロセッサは、前記組織を通過した多重散乱光の伝播をモデル化した式から、前記複数の値を決定するよう構成されていることを特徴とするシステム。
10. 請求項２記載のシステムにおいて、前記光源は、所定の周波数で光の強度を変調することにより、時間によって変動する強度の光を生じるよう構成されていることを特徴とするシステム。
11. 請求項１記載のシステムにおいて、前記プロセッサは、前記光放出の測定された振幅又は測定された位相の少なくとも一方に基づいて、前記蛍光特性を決定するよう構成されていることを特徴とするシステム。

Images