光子輸送のためのモンテカルロ法

生体組織の光学的特性は、生物医学イメージングへのアプローチを提供します。血液やメラニンによる吸収、神経細胞や癌細胞核による散乱など、多くの内因性コントラストが存在します。さらに、蛍光プローブはさまざまな組織に標的化できます。顕微鏡技術（共焦点顕微鏡、二光子顕微鏡、光コヒーレンストモグラフィーなど）は、これらの特性を高空間分解能で画像化できますが、弾道光子に依存するため、深部への浸透は数ミリメートルに制限されます。光子が多重散乱される組織の深部を画像化するには、そのような環境における多数の光子の統計的挙動をより深く理解する必要があります。モンテカルロ法は、組織深部の光学的特性を再構築するためにさまざまな技術で使用されてきた柔軟なフレームワークを提供します。ここでは、これらの技術のいくつかについて簡単に紹介します。

• 光音響トモグラフィー（PAT）では、拡散レーザー光が吸収され、局所的な温度上昇が生じます。この局所的な温度変化は、熱弾性膨張によって超音波を発生させ、超音波トランスデューサによって検出されます。実際には、さまざまな設定パラメータ（光の波長、トランスデューサの開口数など）が変化するため、モンテカルロ法によるモデリングは、実験的手法を用いる前に組織の反応を予測するための貴重なツールとなります。
• 拡散光トモグラフィー（DOT）は、近赤外光源と検出器のアレイを用いて生体組織の光学的特性を測定するイメージング技術です。オキシヘモグロビンとデオキシヘモグロビンによる吸収（機能的神経イメージングや癌検出用）や蛍光プローブの濃度など、さまざまなコントラストを測定できます。画像を再構成するには、特定の光源から特定の検出器まで光がどのように伝わったか、そして測定値が光学的特性の分布と変化にどのように依存するか（順モデルと呼ばれる）を知る必要があります。生体組織は散乱性が高いため、このような経路は複雑で、感度関数は拡散的です。順モデルは、モンテカルロ法を用いて生成されることがよくあります。

放射線治療の目的は、周囲の正常組織を温存しながら、一般的には電離放射線の形でエネルギーを癌組織に照射することです。モンテカルロ法は、放射線治療において、患者組織からの散乱と、リニアアクセラレータの上流にあるコリメータからの散乱の両方によって患者が受ける周辺線量を決定するために一般的に用いられます。

光線力学療法（PDT）では、光を用いて化学療法剤を活性化します。PDTの特性上、化学療法剤を活性化するために適切な量の光が照射されるように、組織内での散乱と吸収をモデル化するためにモンテカルロ法を用いることが有効です。

本モデルでは、屈折率が一致しない媒質に入射する際の初期鏡面反射は無視します。そのため、光子パケットの初期位置と初期方向を設定するだけで済みます。グローバル座標系を用いると便利です。位置を決定するために3つの直交座標を用い、伝搬方向を決定するために3つの方向余弦を用います。初期開始条件は用途によって異なりますが、原点で初期化されたペンシルビームの場合、初期位置と方向余弦を次のように設定できます（等方性光源は、各パケットの初期方向をランダム化することで容易にモデル化できます）。

${\displaystyle {\begin{aligned}x&=0\\{\text{位置: }}y&=0\\z&=0\\\\\mu _{x}&=0\\{\text{方向余弦: }}\mu _{y}&=0\\\mu _{z}&=1\end{aligned}}}$

ステップサイズsは、光子パケットが相互作用サイト間を移動する距離です。ステップサイズの選択にはさまざまな方法があります。以下は、光子ステップサイズ選択の基本形式（逆分布法とベール・ランバートの法則を用いて導出）であり、これを均質モデルに使用します。

${\displaystyle s=-{\frac {\ln \xi }{\mu _{t}}}}$

どこ${\displaystyle \xi }$は乱数であり、${\displaystyle {\mu _{t}}}$は全相互作用係数（すなわち、吸収係数と散乱係数の合計）です。

ステップサイズが選択されると、光子パケットは方向余弦で定義される方向に距離sだけ伝播されます。これは、座標を次のように更新するだけで簡単に実現できます。

${\displaystyle {\begin{aligned}x&\leftarrow x+\mu _{x}s\\y&\leftarrow y+\mu _{y}s\\z&\leftarrow z+\mu _{z}s\end{aligned}}}$

光子の重みの一部は、各相互作用部位で吸収される。この重みの割合は、次のように決定される。

${\displaystyle \Delta W={\frac {\mu _{a}}{\mu _{t}}}W}$

どこ${\displaystyle {\mu _{a}}}$は吸収係数です。

吸収分布が特定の研究にとって重要な場合、重量分率は配列に記録することができます。その後、光子パケットの重量は次のように更新する必要があります。

${\displaystyle W\leftarrow W-\Delta W\,}$

吸収後、光子パケットは散乱されます。光子散乱角のコサインの加重平均は散乱異方性 ( g ) と呼ばれ、その値は-1から 1 の間です。光学異方性が 0 の場合、これは一般的に散乱が等方的であることを示します。g の値が 1 に近づくと、これは散乱が主に前方方向であることを示します。光子パケットの新しい方向 (したがって光子方向コサイン) を決定するには、散乱位相関数を知る必要があります。多くの場合、Henyey-Greenstein 位相関数が使用されます。次に、散乱角 θ は次の式を使用して決定されます。

${\displaystyle \cos \theta ={\begin{cases}{\frac {1}{2g}}\left[1+g^{2}-\left({\frac {1-g^{2}}{1-g+2g\xi }}\right)^{2}\right]&{\text{ if }}g\neq 0\\1-2\xi &{\text{ if }}g=0\end{cases}}}$

また、極角φは一般的に0から100の間で一様に分布していると仮定される。${\displaystyle 2\pi }$この仮定に基づき、以下のように設定できます。

${\displaystyle \varphi =2\pi \xi {\frac {}{}}}$

Based on these angles and the original direction cosines, we can find a new set of direction cosines. The new propagation direction can be represented in the global coordinate system as follows:

${\displaystyle {\begin{aligned}\mu '_{x}&={\frac {\sin \theta (\mu _{x}\mu _{z}\cos \varphi -\mu _{y}\sin \varphi )}{\sqrt {1-\mu _{z}^{2}}}+\mu _{x}\cos \theta \\\mu '_{y}&={\frac {\sin \theta (\mu _{y}\mu _{z}\cos \varphi +\mu _{x}\sin \varphi )}{\sqrt {1-\mu _{z}^{2}}}}+\mu _{y}\cos \theta \\\mu '_{z}&=-{\sqrt {1-\mu _{z}^{2}}}\sin \theta \cos \varphi +\mu _{z}\cos \theta \\\end{整列}}}$

For a special case

${\displaystyle {\begin{aligned}\mu _{z}=1\end{aligned}}}$

use

${\displaystyle {\begin{aligned}\mu '_{x}&=\sin \theta \cos \varphi \\\mu '_{y}&=\sin \theta \sin \varphi \\\mu '_{z}&=\cos \theta \\\end{aligned}}}$

or

${\displaystyle {\begin{aligned}\mu _{z}=-1\end{aligned}}}$

use

${\displaystyle {\begin{aligned}\mu '_{x}&=\sin \theta \cos \varphi \\\mu '_{y}&=-\sin \theta \sin \varphi \\\mu '_{z}&=-\cos \theta \\\end{aligned}}}$

C-code:

/*********************** Indicatrix ********************* *New direction cosines after scattering by angle theta, fi. * mux new=(sin(theta)*(mux*muz*cos(fi)-muy*sin(fi)))/sqrt(1-muz^2)+mux*cos(theta) * muy new=(sin(theta)*(muy*muz*cos(fi)+mux*sin(fi)))/sqrt(1-muz^2)+muy*cos(theta) * muz new= - sqrt(1-muz^2)*sin(theta)*cos(fi)+muz*cos(theta) *--------------------------------------------------------- *Input: * muxs,muys,muzs - direction cosine before collision * mutheta, fi - cosine of polar angle and the azimuthal angle *--------------------------------------------------------- *Output: * muxd,muyd,muzd - direction cosine after collision *--------------------------------------------------------- */ void Indicatrix (double muxs, double muys, double muzs, double mutheta, double fi, double *muxd, double *muyd, double *muzd) { double costheta = mutheta; double sintheta = sqrt(1.0-costheta*costheta); // sin(theta) double sinfi = sin(fi); double cosfi = cos(fi); if (muzs == 1.0) { *muxd = sintheta*cosfi; *muyd = sintheta*sinfi; *muzd = costheta; } elseif (muzs == -1.0) { *muxd = sintheta*cosfi; *muyd = -sintheta*sinfi; *muzd = -costheta; } else { double denom = sqrt(1.0-muzs*muzs); double muzcosfi = muzs*cosfi; *muxd = sintheta*(muxs*muzcosfi-muys*sinfi)/denom + muxs*costheta; *muyd = sintheta*(muys*muzcosfi+muxs*sinfi)/denom + muys*costheta; *muzd = -denom*sintheta*cosfi + muzs*costheta; } }

If a photon packet has experienced many interactions, for most applications the weight left in the packet is of little consequence. As a result, it is necessary to determine a means for terminating photon packets of sufficiently small weight. A simple method would use a threshold, and if the weight of the photon packet is below the threshold, the packet is considered dead. The aforementioned method is limited as it does not conserve energy. To keep total energy constant, a Russian roulette technique is often employed for photons below a certain weight threshold. This technique uses a roulette constant m to determine whether or not the photon will survive. The photon packet has one chance in m to survive, in which case it will be given a new weight of mW where W is the initial weight (this new weight, on average, conserves energy). All other times, the photon weight is set to 0 and the photon is terminated. This is expressed mathematically below:

${\displaystyle W={\begin{cases}mW&\xi \leq 1/m\\0&\xi >1/m\end{cases}}}$