Método de Monte Carlo para el transporte de fotones

Aplicación de modelado

El modelado de la propagación de fotones con métodos de Monte Carlo es un enfoque flexible pero riguroso para simular el transporte de fotones. En el método, las reglas locales del transporte de fotones se expresan como distribuciones de probabilidad que describen el tamaño del paso del movimiento de fotones entre los sitios de interacción fotón-materia y los ángulos de desviación en la trayectoria de un fotón cuando ocurre un evento de dispersión. Esto es equivalente a modelar el transporte de fotones analíticamente mediante la ecuación de transferencia radiativa (RTE), que describe el movimiento de fotones utilizando una ecuación diferencial. Sin embargo, las soluciones de forma cerrada de la RTE a menudo no son posibles; para algunas geometrías, la aproximación de difusión se puede utilizar para simplificar la RTE, aunque esto, a su vez, introduce muchas imprecisiones, especialmente cerca de fuentes y límites. En contraste, las simulaciones de Monte Carlo se pueden hacer arbitrariamente precisas aumentando el número de fotones trazados. Por ejemplo, vea la película, donde una simulación de Monte Carlo de un haz de lápiz incidente en un medio semi-infinito modela tanto el flujo de fotones balístico inicial como la propagación difusa posterior.

El método de Monte Carlo es necesariamente estadístico y, por lo tanto, requiere un tiempo de cálculo significativo para lograr precisión. Además, las simulaciones de Monte Carlo pueden realizar un seguimiento de múltiples cantidades físicas simultáneamente, con cualquier resolución espacial y temporal deseada. Esta flexibilidad hace que el modelado de Monte Carlo sea una herramienta poderosa. Por lo tanto, si bien son ineficientes desde el punto de vista computacional, los métodos de Monte Carlo a menudo se consideran el estándar para las mediciones simuladas del transporte de fotones para muchas aplicaciones biomédicas.

Simulación de Monte Carlo de un haz de lápiz incidente en un medio de dispersión semiinfinito.

Aplicaciones biomédicas de los métodos de Monte Carlo

Imágenes biomédicas

Las propiedades ópticas del tejido biológico ofrecen un enfoque para la obtención de imágenes biomédicas. Existen muchos contrastes endógenos, incluida la absorción de la sangre y la melanina y la dispersión de las células nerviosas y los núcleos de las células cancerosas. Además, las sondas fluorescentes pueden dirigirse a muchos tejidos diferentes. Las técnicas de microscopía (incluida la confocal , la de dos fotones y la tomografía de coherencia óptica ) tienen la capacidad de obtener imágenes de estas propiedades con alta resolución espacial, pero, como dependen de los fotones balísticos, su profundidad de penetración está limitada a unos pocos milímetros. Obtener imágenes más profundas en los tejidos, donde los fotones se han dispersado de forma múltiple, requiere una comprensión más profunda del comportamiento estadístico de grandes cantidades de fotones en un entorno de este tipo. Los métodos de Monte Carlo proporcionan un marco flexible que se ha utilizado con diferentes técnicas para reconstruir las propiedades ópticas en las profundidades del tejido. Aquí se presenta una breve introducción a algunas de estas técnicas.

• Tomografía fotoacústica En la tomografía fotoacústica, se absorbe la luz láser difusa, lo que genera un aumento de la temperatura local. Esta variación de la temperatura local genera, a su vez, ondas ultrasónicas a través de la expansión termoelástica que se detectan mediante un transductor ultrasónico. En la práctica, se varían diversos parámetros de configuración (es decir, longitud de onda de la luz, apertura numérica del transductor) y, como resultado, el modelado de Monte Carlo es una herramienta valiosa para predecir la respuesta del tejido antes de los métodos experimentales.
• La tomografía óptica difusa (DOT) es una técnica de obtención de imágenes que utiliza una serie de fuentes de luz y detectores de infrarrojo cercano para medir las propiedades ópticas de los tejidos biológicos. Se pueden medir diversos contrastes, incluida la absorción debida a la oxihemoglobina y la desoxihemoglobina (para la neuroimagen funcional o la detección del cáncer) y la concentración de sondas fluorescentes. Para reconstruir una imagen, se debe conocer la forma en que la luz viaja desde una fuente determinada hasta un detector determinado y cómo la medición depende de la distribución y los cambios en las propiedades ópticas (conocido como el modelo directo). Debido a la naturaleza altamente dispersante del tejido biológico, dichas trayectorias son complicadas y las funciones de sensibilidad son difusas. El modelo directo se genera a menudo utilizando métodos de Monte Carlo.

Radioterapia

El objetivo de la radioterapia es suministrar energía, generalmente en forma de radiación ionizante, al tejido canceroso sin afectar al tejido normal circundante. El modelado de Monte Carlo se emplea habitualmente en radioterapia para determinar la dosis periférica que experimentará el paciente debido a la dispersión, tanto del tejido del paciente como de la colimación anterior en el acelerador lineal.

Terapia fotodinámica

En la terapia fotodinámica (TFD), se utiliza luz para activar los agentes quimioterapéuticos. Debido a la naturaleza de la TFD, resulta útil utilizar métodos de Monte Carlo para modelar la dispersión y la absorción en el tejido a fin de garantizar que se administren los niveles adecuados de luz para activar los agentes quimioterapéuticos.

Implementación del transporte de fotones en un medio dispersor

Aquí se presenta un modelo de un método de Monte Carlo de fotones en un medio infinito homogéneo. Sin embargo, el modelo se puede extender fácilmente para medios multicapa. Para un medio no homogéneo, se deben considerar los límites. Además, para un medio semi-infinito (en el que los fotones se consideran perdidos si salen del límite superior), se debe tener una consideración especial. Para obtener más información, visite los enlaces en la parte inferior de la página. Resolveremos el problema utilizando una fuente puntual infinitamente pequeña (representada analíticamente como una función delta de Dirac en el espacio y el tiempo). Las respuestas a geometrías de fuente arbitrarias se pueden construir utilizando el método de funciones de Green (o convolución , si existe suficiente simetría espacial). Los parámetros requeridos son el coeficiente de absorción , el coeficiente de dispersión y la función de fase de dispersión. (Si se consideran los límites, también se debe proporcionar el índice de refracción para cada medio). Las respuestas resueltas en el tiempo se encuentran haciendo un seguimiento del tiempo total transcurrido del vuelo del fotón utilizando la longitud del camino óptico . Las respuestas a fuentes con perfiles de tiempo arbitrarios pueden luego modelarse a través de convolución en el tiempo.

En nuestro modelo simplificado, utilizamos la siguiente técnica de reducción de varianza para reducir el tiempo de cálculo. En lugar de propagar fotones individualmente, creamos un paquete de fotones con un peso específico (generalmente inicializado como unidad). A medida que el fotón interactúa en el medio turbio, depositará peso debido a la absorción y el peso restante se dispersará a otras partes del medio. Se puede registrar cualquier cantidad de variables a lo largo del proceso, según el interés de una aplicación particular. Cada paquete de fotones pasará repetidamente por los siguientes pasos numerados hasta que se termine, se refleje o se transmita. El proceso se muestra en el esquema de la derecha. Se puede lanzar y modelar cualquier cantidad de paquetes de fotones, hasta que las mediciones simuladas resultantes tengan la relación señal-ruido deseada. Tenga en cuenta que, como el modelado de Monte Carlo es un proceso estadístico que involucra números aleatorios, utilizaremos la variable ξ en todo momento como un número pseudoaleatorio para muchos cálculos.

Esquema para modelar el flujo de fotones en un medio de dispersión y absorción infinitas con simulaciones de Monte Carlo.

Paso 1: Lanzamiento de un paquete de fotones

En nuestro modelo, ignoramos la reflectancia especular inicial asociada con la entrada a un medio que no tiene un índice de refracción adecuado. Con esto en mente, simplemente necesitamos establecer la posición inicial del paquete de fotones, así como la dirección inicial. Es conveniente utilizar un sistema de coordenadas global. Usaremos tres coordenadas cartesianas para determinar la posición, junto con tres cosenos de dirección para determinar la dirección de propagación. Las condiciones iniciales de inicio variarán según la aplicación, sin embargo, para un haz en forma de lápiz inicializado en el origen, podemos establecer la posición inicial y los cosenos de dirección de la siguiente manera (las fuentes isotrópicas se pueden modelar fácilmente al aleatorizar la dirección inicial de cada paquete):

${\displaystyle {\begin{aligned}x&=0\\{\text{Position: }}y&=0\\z&=0\\\\\mu _{x}&=0\\{\text{Direction cosines: }}\mu _{y}&=0\\\mu _{z}&=1\end{aligned}}}$

Paso 2: Selección del tamaño del paso y movimiento del paquete de fotones

El tamaño del paso, s , es la distancia que recorre el paquete de fotones entre los sitios de interacción. Existen diversos métodos para seleccionar el tamaño del paso. A continuación, se muestra una forma básica de selección del tamaño del paso de fotones (derivada mediante el método de distribución inversa y la ley de Beer-Lambert ) que utilizamos para nuestro modelo homogéneo:

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

donde es un número aleatorio y es el coeficiente de interacción total (es decir, la suma de los coeficientes de absorción y dispersión). ${\displaystyle \xi }$ ${\displaystyle {\mu _{t}}}$

Una vez que se selecciona un tamaño de paso, el paquete de fotones se propaga a lo largo de una distancia s en una dirección definida por los cosenos directores. Esto se logra fácilmente simplemente actualizando las coordenadas de la siguiente manera:

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

Paso 3: Absorción y dispersión

En cada sitio de interacción se absorbe una parte del peso del fotón. Esta fracción del peso se determina de la siguiente manera:

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

donde es el coeficiente de absorción. ${\displaystyle {\mu _{a}}}$

La fracción de peso se puede registrar en una matriz si la distribución de absorción es de interés para el estudio en particular. El peso del paquete de fotones se debe actualizar de la siguiente manera:

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

Después de la absorción, el paquete de fotones se dispersa. El promedio ponderado del coseno del ángulo de dispersión del fotón se conoce como anisotropía de dispersión ( g ), que tiene un valor entre −1 y 1. Si la anisotropía óptica es 0, esto generalmente indica que la dispersión es isotrópica. Si g se acerca a un valor de 1, esto indica que la dispersión se produce principalmente en la dirección hacia delante. Para determinar la nueva dirección del paquete de fotones (y, por lo tanto, los cosenos de dirección de los fotones), necesitamos conocer la función de fase de dispersión. A menudo se utiliza la función de fase de Henyey-Greenstein. Luego, el ángulo de dispersión, θ, se determina utilizando la siguiente fórmula.

${\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}}}$

Además, se supone que el ángulo polar φ se distribuye uniformemente entre 0 y . Con base en esta suposición, podemos establecer: ${\displaystyle 2\pi }$

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

A partir de estos ángulos y de los cosenos directores originales, podemos encontrar un nuevo conjunto de cosenos directores. La nueva dirección de propagación se puede representar en el sistema de coordenadas global de la siguiente manera:

${\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{aligned}}}$

Para un caso especial

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

usar

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

o

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

usar

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

Código C:

/*********************** Indicatriz **********************Nuevos cosenos de dirección después de la dispersión por el ángulo theta, fi.* mux nuevo=(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 nuevo= - sqrt(1-muz^2)*sin(theta)*cos(fi)+muz*cos(theta)*------------------------------------------------- --------*Aporte:* muxs,muys,muzs - coseno de dirección antes de la colisión* mutheta, fi - coseno del ángulo polar y del ángulo azimutal*------------------------------------------------- --------*Producción:* muxd,muyd,muzd - coseno de dirección después de la colisión*------------------------------------------------- --------*/Indicatriz vacía (doble muxs, doble muys, doble muzs, doble mutheta, doble fi, doble *muxd, doble *muyd, doble *muzd){ doble costheta = mutheta; doble sinteta = sqrt(1.0-costheta*costheta); // pecado(theta) doble sinfi = sin(fi); doble cosfi = cos(fi); si (muzs == 1.0) { *muxd = sintheta*cosfi; *muyd = sintheta*sinfi; *muzd = costotheta; } elseif (muzs == -1.0) { *muxd = sintheta*cosfi; *muyd = -sintheta*sinfi; *muzd = -costheta; } demás { doble denominación = sqrt(1.0-muzs*muzs); doble 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; }}

Paso 4: Terminación del fotón

Si un paquete de fotones ha experimentado muchas interacciones, para la mayoría de las aplicaciones el peso que queda en el paquete tiene poca importancia. Como resultado, es necesario determinar un medio para terminar los paquetes de fotones de peso suficientemente pequeño. Un método simple utilizaría un umbral, y si el peso del paquete de fotones está por debajo del umbral, el paquete se considera muerto. El método mencionado anteriormente es limitado ya que no conserva la energía. Para mantener constante la energía total, a menudo se emplea una técnica de ruleta rusa para fotones por debajo de un cierto umbral de peso. Esta técnica utiliza una constante de ruleta m para determinar si el fotón sobrevivirá o no. El paquete de fotones tiene una probabilidad en m de sobrevivir, en cuyo caso se le dará un nuevo peso de mW donde W es el peso inicial (este nuevo peso, en promedio, conserva la energía). En todos los demás casos, el peso del fotón se establece en 0 y el fotón se termina. Esto se expresa matemáticamente a continuación:

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

Unidades de procesamiento gráfico (GPU) y simulaciones rápidas de Monte Carlo del transporte de fotones

La simulación de Monte Carlo de la migración de fotones en medios turbios es un problema altamente paralelizable, en el que una gran cantidad de fotones se propagan de forma independiente, pero de acuerdo con reglas idénticas y diferentes secuencias de números aleatorios. La naturaleza paralela de este tipo especial de simulación de Monte Carlo lo hace muy adecuado para su ejecución en una unidad de procesamiento gráfico (GPU). El lanzamiento de las GPU programables inició este desarrollo y desde 2008 ha habido algunos informes sobre el uso de GPU para la simulación de Monte Carlo de alta velocidad de la migración de fotones. ^{[1] [2] [3] [4]}

Este enfoque básico se puede paralelizar mediante el uso de varias GPU conectadas entre sí. Un ejemplo es el "GPU Cluster MCML", que se puede descargar desde el sitio web de los autores (Monte Carlo Simulation of Light Transport in Multi-layered Turbid Media Based on GPU Clusters): http://bmp.hust.edu.cn/GPU_Cluster/GPU_Cluster_MCML.HTM

Véase también

• Ecuación de transferencia radiativa y teoría de difusión para el transporte de fotones en tejido biológico
• Método de Monte Carlo
• Convolución para respuestas ópticas de haz ancho en medios de dispersión
• Métodos de Monte Carlo para el transporte de electrones

Enlaces a otros recursos de Monte Carlo

• Laboratorio de imágenes ópticas de la Universidad de Washington en St. Louis (MCML)
• Centro de láser médico de Oregón
• Investigación sobre migración de fotones en Monte Carlo en la Universidad de Lund, Suecia. Aceleración por GPU de simulaciones de Monte Carlo y Monte Carlo escalable. Código fuente abierto para descargar.
• Monte Carlo basado en la nube para el transporte de luz en un medio de dispersión turbio La herramienta es de uso gratuito en investigaciones y actividades no comerciales.
• Transporte de luz en el tejido como ejemplo de simulación de Monte Carlo (con código fuente C++).

Referencias

• Wang, LH; Wu Hsin-I (2007). Óptica biomédica: principios e imágenes . Wiley.
• L.-H. Wang; SL Jacques; L.-Q. Zheng (1995). "MCML—Modelado de Monte Carlo del transporte de luz en tejidos multicapa". Métodos y programas informáticos en biomedicina . 47 (2): 131–146. doi :10.1016/0169-2607(95)01640-F. PMID  7587160.
• L.-H. Wang; SL Jacques; L.-Q. Zheng (1997). "Convolución de convolución para respuestas a un haz de fotones de diámetro finito incidente en tejidos multicapa" (PDF) . Métodos y programas informáticos en biomedicina . 54 (3): 141–150. doi :10.1016/S0169-2607(97)00021-7. PMID  9421660.
• SL Jacques; L.-H. Wang (1995). "Modelado de Monte Carlo del transporte de luz en tejidos" (PDF) . En AJ Welch; MJC van Gemert (eds.). Respuesta térmica óptica del tejido irradiado con láser . Nueva York: Plenum Press. págs. 73–100.
• L.-H. Wang; SL Jacques (1994). "Posiciones radiales y angulares optimizadas en el modelado de Monte Carlo" (PDF) . Física Médica . 21 (7): 1081–1083. Bibcode :1994MedPh..21.1081W. doi :10.1118/1.597351. PMID  7968840.

Referencias en línea

1. E. Alerstam; T. Svensson; S. Andersson-Engels (2008). "Computación paralela con unidades de procesamiento gráfico para simulación Monte Carlo de alta velocidad de migración de fotones" (PDF) . J. Biomed. Opt . 13 (6): 060504. Bibcode :2008JBO....13f0504A. doi : 10.1117/1.3041496 . PMID  19123645.
2. Q. Fang; DA Boas (2009). "Simulación de Monte Carlo de la migración de fotones en medios turbios tridimensionales acelerada por unidades de procesamiento gráfico". Opt. Express . 17 (22): 20178–20190. Bibcode :2009OExpr..1720178F. doi :10.1364/oe.17.020178. PMC 2863034 . PMID  19997242. 
3. N. Ren; J. Liang; X. Qu; J. Li; B. Lu; J. Tian (2010). "Simulación de Monte Carlo basada en GPU para propagación de luz en tejidos heterogéneos complejos". Opt. Express . 18 (7): 6811–6823. Bibcode :2010OExpr..18.6811R. doi : 10.1364/oe.18.006811 . PMID  20389700.
4. Alexander Doronin; Igor Meglinski (2011). "Herramienta computacional Monte Carlo orientada a objetos en línea para las necesidades de la óptica biomédica". Biomed. Opt. Express . 2 (9): 2461–2469. doi :10.1364/boe.2.002461. PMC 3184856 . PMID  21991540.