生物組織における光子輸送に関する放射伝達方程式と拡散理論

生物組織における光子輸送は、モンテカルロシミュレーションで数値的にモデル化することも、放射伝達方程式 (RTE) で解析的にモデル化することもできます。ただし、RTE は近似を導入せずに解くのは困難です。ここでまとめた一般的な近似は拡散近似です。全体的に、光子輸送の拡散方程式の解は、モンテカルロシミュレーションよりも計算効率は高いですが、精度は劣ります。^[1]

RTE は、光子が組織内を移動する際のエネルギー伝達を数学的にモデル化できます。放射場内の小面積要素を通る放射エネルギーの流れは、単位 の 放射輝度で特徴付けることができます。放射輝度は、単位法線面積、単位 立体角、単位時間あたりのエネルギーの流れとして定義されます。ここで、は位置、は単位方向ベクトル、 は時間を表します (図 1)。 他のいくつかの重要な物理量は、放射輝度の定義に基づいています。^[1]${\displaystyle L({\vec {r}},{\hat {s}},t)}$${\displaystyle {\frac {\mathrm {W} }{\mathrm {m} ^{2}\mathrm {sr} }}}$${\displaystyle {\vec {r}}}$${\displaystyle {\hat {s}}}$${\displaystyle t}$

• フルエンス率または強度 ${\displaystyle \Phi ({\vec {r}},t)=\int _{4\pi }L({\vec {r}},{\hat {s}},t)d\Omega \quad \left[{\frac {\mathrm {W} }{\mathrm {m} ^{2}\mathrm {sr} }}\right]}$
• フルエンス ${\displaystyle F({\vec {r}})=\int _{-\infty }^{+\infty }\Phi ({\vec {r}},t)dt\quad \left[{\frac {\mathrm {J} }{\mathrm {m} ^{2}}}\right]}$
• 電流密度（エネルギー流束）。これは、エネルギーの流れの一般的な方向を指すフルエンス率のベクトル対応です。${\displaystyle {\vec {J}}({\vec {r}},t)=\int _{4\pi}{\hat {s}}L({\vec {r}},{\hat {s}},t)d\Omega \quad \left[{\frac {\mathrm {W}}{\mathrm {m}^{2}}}\right]}$

RTE は、放射輝度 を記述する微分方程式です。これは、エネルギー保存則によって導出できます。簡単に言うと、RTE は、光線が発散と消滅(光線からの吸収と散乱の両方を含む) によってエネルギーを失い、媒体内の光源と光線に向かう散乱からエネルギーを得ることを述べています。コヒーレンス、偏光、非線形性は無視されます。屈折率、吸収係数μ _a、散乱係数 μ _s、散乱異方性などの光学特性は、時間的に不変ですが、空間的に変化する可能性があります。散乱は弾性であると仮定されます。したがって、RTE (ボルツマン方程式) は次のように記述されます。^[1]${\displaystyle L({\vec {r}},{\hat {s}},t)}$ ${\displaystyle n}$${\displaystyle g}$

${\displaystyle {\frac {\partial L({\vec {r}},{\hat {s}},t)/c}{\partial t}}=-{\hat {s}}\cdot \nabla L({\vec {r}},{\hat {s}},t)-\mu _{t}L({\vec {r}},{\hat {s}},t)+\mu _{s}\int _{4\pi }L({\vec {r}},{\hat {s}}',t)P({\hat {s}}',{\hat {s}})d\Omega '+S({\vec {r}},{\hat {s}},t)}$

どこ

• ${\displaystyle c}$組織内の光の速度であり、相対屈折率によって決定される。
• μ _t μ _a +μ _{s は}消衰係数です${\displaystyle =}$
• ${\displaystyle P({\hat {s}}',{\hat {s}})}$は位相関数であり、伝播方向 の光が立体角の周囲に散乱される確率を表す。ほとんどの場合、位相関数は散乱方向と入射方向の間の角度、すなわちのみに依存する。散乱異方性は次のように表される。${\displaystyle {\hat {s}}'}$${\displaystyle d\Omega }$${\displaystyle {\hat {s}}}$${\displaystyle {\hat {s}}'}$${\displaystyle {\hat {s}}}$${\displaystyle P({\hat {s}}',{\hat {s}})=P({\hat {s}}'\cdot {\hat {s}})}$${\displaystyle g=\int _{4\pi }({\hat {s}}'\cdot {\hat {s}})P({\hat {s}}'\cdot {\hat {s}})d\Omega }$
• ${\displaystyle S({\vec {r}},{\hat {s}},t)}$光源について説明します。

RTE では、6 つの異なる独立変数によって、任意の空間的および時間的ポイント ( からの 、 、および 、極角と方位角 、および ) での放射輝度が定義されます。散乱媒体内の光子の挙動について適切な仮定を行うことで、独立変数の数を減らすことができます。これらの仮定から、光子輸送の拡散理論(および拡散方程式) が導き出されます。2 つの仮定により、拡散理論を RTE に適用できます。 ${\displaystyle x}$${\displaystyle y}$${\displaystyle z}$${\displaystyle {\vec {r}}}$ ${\displaystyle \theta}$${\displaystyle \phi}$${\displaystyle {\hat {s}}}$${\displaystyle t}$

• 散乱イベントに比べると、吸収イベントは非常に少ないです。同様に、多数の散乱イベントの後には、吸収イベントはほとんど発生せず、放射輝度はほぼ等方性になります。この仮定は、方向の広がりと呼ばれることもあります。
• 主に散乱する媒体では、電流密度の大幅な変化にかかる時間は、1 つの輸送平均自由経路を通過する時間よりもはるかに長くなります。したがって、1 つの輸送平均自由経路では、電流密度の分数変化は 1 よりはるかに小さくなります。この特性は、時間的広がりと呼ばれることもあります。

これらの仮定は両方とも、高アルベド（主に散乱）媒体を必要とする。^[1]

放射輝度は、球面調和関数 _n、mの基底関数で展開できます。拡散理論では、放射輝度はほぼ等方性であると考えられるため、等方性項と 1 次の異方性項のみが使用されます。 ここで、_{n、m は}展開係数です。放射輝度は 4 つの項で表されます。n = 0 の場合は 1 つ (等方性項)、n = 1 の場合は 3 つ (異方性項) です。球面調和関数の特性と、フルエンス率および電流密度の定義を使用すると、等方性項と異方性項はそれぞれ次のように表すことができます。 ${\displaystyle Y}$${\displaystyle L({\vec {r}},{\hat {s}},t)\approx \\sum _{n=0}^{1}\sum _{m=-n}^{n}L_{n,m}({\vec {r}},t)Y_{n,m}({\hat {s}})}$${\displaystyle L}$${\displaystyle \Phi ({\vec {r}},t)}$${\displaystyle {\vec {J}}({\vec {r}},t)}$

• ${\displaystyle L_{0,0}({\vec {r}},t)Y_{0,0}({\hat {s}})={\frac {\Phi ({\vec {r}} ,t)}{4\pi }}}$
• ${\displaystyle \sum _{m=-1}^{1}L_{1,m}({\vec {r}},t)Y_{1,m}({\hat {s}})={\frac {3}{4\pi }}{\vec {J}}({\vec {r}},t)\cdot {\hat {s}}}$

したがって、放射輝度は次のように近似できる^[1]

${\displaystyle L({\vec {r}},{\hat {s}},t)={\frac {1}{4\pi }}\Phi ({\vec {r}},t)+ {\frac {3}{4\pi }}{\vec {J}}({\vec {r}},t)\cdot {\hat {s}}}$

上記の式を放射輝度に代入すると、RTEはそれぞれスカラー形式とベクトル形式で次のように書き直すことができます（RTEの散乱項は完全な立体角にわたって積分されます。ベクトル形式では、評価前にRTEに方向を掛けます）。^[1]${\displaystyle 4\pi }$${\displaystyle {\hat {s}}}$

${\displaystyle {\frac {\partial \Phi ({\vec {r}},t)}{c\partial t}}+\mu _{a}\Phi ({\vec {r}},t) +\nabla \cdot {\vec {J}}({\vec {r}},t)=S({\vec {r}},t)}$

${\displaystyle {\frac {\partial {\vec {J}}({\vec {r}},t)}{c\partial t}}+(\mu _{a}+\mu _{s} '){\vec {J}}({\vec {r}},t)+{\frac {1}{3}}\nabla \Phi ({\vec {r}},t)=0}$

拡散近似は、縮減散乱係数が吸収係数よりもはるかに大きく、最小層厚が輸送平均自由行程の数倍のオーダーであるシステムに限定されます。

拡散理論の2番目の仮定を用いると、 1輸送平均自由行程における電流密度の分数変化は無視できることがわかる。拡散理論RTEのベクトル表現はフィックの法則に帰着し、これは電流密度を流束率の勾配で定義する。フィックの法則をRTEのスカラー表現に代入すると、拡散方程式が得られる: ^[1]${\displaystyle {\vec {J}}({\vec {r}},t)}$ ${\displaystyle {\vec {J}}({\vec {r}},t)={\frac {-\nabla \Phi ({\vec {r}},t)}{3(\mu _{a}+\mu _{s}')}}}$

${\displaystyle {\frac {1}{c}}{\frac {\partial \Phi ({\vec {r}},t)}{\partial t}}+\mu _{a}\Phi ({\vec {r}},t)-\nabla \cdot [D\nabla \Phi ({\vec {r}},t)]=S({\vec {r}},t)}$

${\displaystyle D={\frac {1}{3(\mu _{a}+\mu _{s}')}}}$は拡散係数、μ' _s μ _{s は}縮減散乱係数です。 注目すべきことに、拡散方程式には散乱係数への明示的な依存性はありません。代わりに、縮減散乱係数のみが の式に現れます。これにより重要な関係が導かれます。縮減散乱係数が一定のままで散乱媒体の異方性が変化しても、拡散は影響を受けません。^[1]${\displaystyle =(1-g)}$
${\displaystyle D}$

境界（組織の層など）と光源のさまざまな構成に対して、適切な境界条件を適用し、状況に応じて 光源項を定義することで、拡散方程式を解くことができます。${\displaystyle S({\vec {r}},t)}$

このセクションでは、無限均質媒体内の短パルス点光源の単純なケースに対する拡散方程式の解を示します。拡散方程式の光源項は となり、ここで はフルエンス率が測定される位置、は光源の位置です。パルスは の時点でピークに達します。拡散方程式をフルエンス率について解くと、拡散方程式の グリーン関数が得られます。${\displaystyle S({\vec {r}},t,{\vec {r}}',t')=\delta ({\vec {r}}-{\vec {r}}')\delta (t-t')}$${\displaystyle {\vec {r}}}$${\displaystyle {\vec {r}}'}$${\displaystyle t'}$

${\displaystyle \Phi ({\vec {r}},t;{\vec {r}}',t)={\frac {c}{[4\pi Dc(t-t')]^{3/2}}}\exp \left[-{\frac {\mid {\vec {r}}-{\vec {r}}'\mid ^{2}}{4Dc(t-t')}}\right]\exp[-\mu _{a}c(t-t')]}$

項は、ビールの法則に従って吸収によるフルエンス率の指数関数的減少を表します。その他の項は散乱による広がりを表します。上記の解から、任意の光源は短パルス点光源の重ね合わせとして特徴付けることができます。拡散方程式から時間変化を取り除くと、時間に依存しない点光源について次の式が得られます。 ${\displaystyle \exp \left[-\mu _{a}c(t-t')\right]}$${\displaystyle S({\vec {r}})=\delta ({\vec {r}})}$

${\displaystyle \Phi ({\vec {r}})={\frac {1}{4\pi Dr}}\exp(-\mu _{\mathrm {eff} }r)}$

${\displaystyle \mu _{\mathrm {eff} }={\sqrt {\frac {\mu _{a}}{D}}}}$は有効減衰係数であり、影響度の空間減衰率を示す。^[1]

境界条件を考慮すると、拡散方程式を使用して、限られたサイズの媒体（媒体と周囲環境との境界を考慮する必要がある）における光の伝播を特徴付けることができます。境界に対処するには、媒体内の光子が境界（つまり表面）に到達したときに何が起こるかを検討します。境界で媒体に向けられた方向積分放射輝度は、境界で媒体から外に向けられた方向積分放射輝度に反射率 を乗じたものに等しくなります。 ${\displaystyle R_{F}}$

${\displaystyle \int _{{\hat {s}}\cdot {\hat {n}}<0}L({\vec {r}},{\hat {s}},t){\hat {s}}\cdot {\hat {n}}d\Omega =\int _{{\hat {s}}\cdot {\hat {n}}>0}R_{F}({\hat {s}}\cdot {\hat {n}})L({\vec {r}},{\hat {s}},t){\hat {s}}\cdot {\hat {n}}d\Omega }$

ここで は境界に垂直で、境界から遠ざかる方向を向いています。拡散近似により、放射輝度はフルエンス率と電流密度で表されます。上記の積分を代入して評価すると、次の式が得られます。^[3]${\displaystyle {\hat {n}}}$${\displaystyle L}$${\displaystyle \Phi }$${\displaystyle {\vec {J}}}$

${\displaystyle {\frac {\Phi ({\vec {r}},t)}{4}}+{\vec {J}}({\vec {r}},t)\cdot {\frac {\hat {n}}{2}}=R_{\Phi }{\frac {\Phi ({\vec {r}},t)}{4}}-R_{J}{\vec {J}}({\vec {r}},t)\cdot {\frac {\hat {n}}{2}}}$

• ${\displaystyle R_{\Phi }=\int _{0}^{\pi /2}2\sin \theta \cos \theta R_{F}(\cos \theta )d\theta }$
• ${\displaystyle R_{J}=\int _{0}^{\pi /2}3\sin \theta (\cos \theta )^{2}R_{F}(\cos \theta )d\theta }$

フィックの法則（）を代入すると、境界からの距離z=0において、^[3]${\displaystyle {\vec {J}}({\vec {r}},t)=-D\nabla \Phi ({\vec {r}},t)}$

${\displaystyle \Phi ({\vec {r}},t)=A_{z}{\frac {\partial \Phi ({\vec {r}},t)}{\partial z}}}$

• ${\displaystyle A_{z}=2D{\frac {1+R_{\mathrm {eff} }}{1-R_{\mathrm {eff} }}}}$
• ${\displaystyle R_{\mathrm {eff} }={\frac {R_{\Phi }+R_{J}}{2-R_{\Phi }+R_{J}}}}$

ゼロフルエンス境界を特定することが望ましい。しかし、物理的境界におけるフルエンス率は、一般にゼロではない。フルエンス率がゼロとなる_bの外挿境界を決定して、画像源を確立することができる。1次テイラー級数近似を使用すると、 ${\displaystyle \Phi (z=0,t)}$${\displaystyle z}$

${\displaystyle \left.\Phi (z=-A_{z},t)\approx \Phi (z=0,t)-A_{z}{\frac {\partial \Phi ({\vec {r}},t)}{\partial z}}\right|_{z=0}}$

これは なのでゼロと評価されます。したがって、定義により、_{b は上で定義されたz}でなければなりません。特に、境界の両側で屈折率が同じ場合、_{F は}ゼロであり、外挿された境界は_bになります。^[3]${\displaystyle \Phi ({\vec {r}},t)=A_{z}{\frac {\partial \Phi ({\vec {r}},t)}{\partial z}}}$${\displaystyle z}$${\displaystyle -A}$${\displaystyle R}$${\displaystyle z}$${\displaystyle =-2D}$

境界条件を使用すると、半無限媒体に垂直に入射するペンシルビームの拡散反射率を近似的に特徴付けることができます。ビームは、次のように無限媒体内の2つの点光源として表されます（図2）：^{[1] [4]}

1. 散乱媒体の散乱異方性を_{2に設定し、新しい散乱係数μ s2を}元のμ _s1に1を掛けた値に設定します。ここで、_1は元の散乱異方性です。${\displaystyle g}$${\displaystyle =0}$${\displaystyle (1-g}$${\displaystyle )}$${\displaystyle g}$
2. ペンシルビームを、表面から 1 輸送平均自由行程下の深さ、電力 = ' で等方性点光源に変換します。${\displaystyle l}$${\displaystyle a}$
3. ' _bの表面上に反対符号のイメージソースを追加して、外挿境界条件を実装します。${\displaystyle l}$${\displaystyle +2z}$

2つの点光源は、無限媒体内の点光源として特徴付けられる。

${\displaystyle \Phi _{\infty }(r,\theta ,z;r',\theta ',z')={\frac {1}{4\pi D\rho }}\exp(-\mu _{\mathrm {eff} }\rho )}$

${\displaystyle \rho }$は観測点から放射源位置までの距離であり、円筒座標で表されます。2つの画像放射源からの放射率寄与の線形結合は ${\displaystyle (r,\theta ,z)}$${\displaystyle (r',\theta ',z')}$

${\displaystyle \Phi (r,\theta ,z;r',\theta ',z')=a'\Phi _{\infty }(r,\theta ,z;r',\theta ',z')-a'\Phi _{\infty }(r,\theta ,z;r',\theta ',-z'-2z_{b})}$

これを使用して、フィックの法則を介して拡散反射率_dを取得できます。 ${\displaystyle R}$${\displaystyle (r)}$

${\displaystyle \left.R_{d}(r)=D{\frac {\partial \Phi }{\partial z}}\right|_{z=0}={\frac {a'z'(1+\mu _{\mathrm {eff} }\rho _{1})\exp(-\mu _{\mathrm {eff} }\rho _{1})}{4\pi \rho _{1}^{3}}}+{\frac {a'(z'+4D)(1+\mu _{\mathrm {eff} }\rho _{2})\exp(-\mu _{\mathrm {eff} }\rho _{2})}{4\pi \rho _{2}^{3}}}}$

${\displaystyle \rho _{1}}$は観測点から画像源までの距離であり、は観測点から画像源_bまでの距離である。^{[1] [4]}${\displaystyle (r,0,0)}$${\displaystyle (0,0,z')}$${\displaystyle \rho _{2}}$${\displaystyle (0,0,-z'-2z}$${\displaystyle )}$

光学的性質、の均質媒体に対する拡散方程式のグリーン関数解を とすると、となる光学的性質、のみが異なる均質媒体に対するグリーン関数解は、次の再スケーリングによって得られる: ^[5]${\displaystyle \Phi ({\vec {r}},t)}$${\displaystyle \mu _{a}}$${\displaystyle \mu _{s}'}$${\displaystyle {\bar {\mu }}_{a}}$${\displaystyle {\bar {\mu }}_{s}'}$${\displaystyle {\bar {\mu }}_{a}/\mu _{a}={\bar {\mu }}_{s}'/\mu _{s}'}$

${\displaystyle {\bar {\Phi }}({\bar {\vec {r}}},{\bar {t}})=\left({\frac {{\bar {\mu }}_{s}'}{\mu _{s}'}}\right)^{3}\Phi ({\vec {r}},t)}$

ここで、および。 ${\displaystyle {\bar {\vec {r}}}={\vec {r}}{\frac {\mu _{s}'}{{\bar {\mu }}_{s}'}}}$${\displaystyle {\bar {t}}=t{\frac {\mu _{s}'}{{\bar {\mu }}_{s}'}}}$

このような特性は、輸送係数、を消衰係数、で置き換えることによって、RTE のより一般的な枠組みにおける放射輝度にも拡張できます。 ${\displaystyle {\bar {\mu }}_{s}'}$${\displaystyle \mu _{s}'}$${\displaystyle {\bar {\mu }}_{t}}$${\displaystyle \mu _{t}}$

この特性の有用性は、実験室規模の設定に典型的な、与えられた形状と光学特性のセットに対して得られた結果を取り、それらを再スケールし、その広大な範囲やアクセス不能性のために測定を実行することが複雑になるような状況に拡張することにある。^[6]

を非吸収均質媒体の拡散方程式のグリーン関数解とします。すると、吸収係数が であるときの媒体のグリーン関数解は次のように得られます: ^[5]${\displaystyle \Phi ({\vec {r}},t,\mu _{a}=0)}$${\displaystyle \mu _{a}}$

${\displaystyle \Phi ({\vec {r}},t,\mu _{a})=\Phi ({\vec {r}},t,\mu _{a}=0)\exp(-\mu _{a}vt)}$

同様に、RTE 内の放射輝度についても同じ特性が当てはまります。

光子輸送のモンテカルロシミュレーションは時間がかかりますが、散乱媒体における光子の挙動を正確に予測します。拡散方程式で光子の挙動を特徴付ける際に必要な仮定は不正確さを生み出します。一般的に、吸収係数 μ _aが増加し、散乱係数 μ _sが減少するにつれて、拡散近似の精度は低下します。^{[7] [8]} 限られた深さの媒体に入射する光子ビームの場合、拡散近似による誤差は、光子の入射位置の 1 つの輸送平均自由行程内 (放射輝度がまだ等方性ではない) で最も顕著になります (図 3)。拡散
方程式を使用して半無限媒体に入射するペンシルビームを記述する手順のうち、媒体を異方性から等方性に変換する (手順 1) (図 4) およびビームを光源に変換する (手順 2) (図 5) は、単一の光源から一対のイメージ光源に変換する (手順 3) (図 6) よりも多くの誤差を生み出します。手順 2 は最も大きな誤差を生み出します。^{[1] [4]}

• 
  図 3: モンテカルロ シミュレーションによって決定された入射ペンシル ビームからの拡散反射率と半径 (赤)、および RTE の拡散理論解によって決定された 2 つの等方性点光源からの拡散反射率と半径 (青)。輸送平均自由行程は 0.1 cm です。
• 
  図 4: 異方性 (青) および等方性 (赤) 媒体における入射ペンシルビームからの拡散反射率と半径の関係。
• 
  図 5: ペンシルビーム (青) と等方性点光源 (赤) の光子源からの拡散反射率と半径の関係。
• 
  図 6: RTE の解 (青) とモンテカルロ シミュレーション (赤) によって特徴付けられた、等方性点光源の光子源からの拡散反射率と半径の関係。

• 光子輸送のためのモンテカルロ法
• 放射伝達

• LV Wang & HI Wu (2007).バイオメディカルオプティクス. Wiley. ISBN 978-0-471-74304-0。

• SG Proskurin (2011). 「量子電子. 41 402」.量子エレクトロニクス. 41 (5): 402–406. doi :10.1070/QE2011v041n05ABEH014597. S2CID  122946781.（2011）