Cálculo de Dosis Monte Carlo en Paciente y Electrones

Cálculo de Dosis en Paciente con Monte Carlo

El cálculo de dosis en el paciente es el punto culminante de toda simulación Monte Carlo (MC) en radioterapia. A diferencia de los modelos analíticos simplificados, el MC rastrea partícula por partícula — y eso lo cambia todo. Los principales códigos disponibles hoy son EGSnrc, MCNP, GEANT4 y PENELOPE, todos basados en la técnica de historia condensada propuesta por Berger en 1963, donde miles de interacciones de electrones se agrupan en pasos macroscópicos.

Existen dos clases de implementación. En la clase I (usada por MCNP/ETRAN), todas las colisiones se agrupan y las partículas secundarias por encima de un umbral energético se crean a posteriori. En la clase II (EGSnrc, GEANT4, PENELOPE), las colisiones «duras» o catastróficas se rastrean individualmente mientras las «blandas» siguen el esquema de clase I. Los artefactos de step-size y los problemas de cruce de fronteras ya están bien comprendidos y resueltos — lo que mejoró significativamente la precisión de los códigos MC modernos. Para una visión completa del tema, consulte nuestra guía completa sobre Monte Carlo en Radioterapia.

Incertidumbres Estadísticas y Geometrías de Scoring

Toda simulación MC introduce ruido estadístico inevitable en la dosis calculada. La fluctuación alrededor del valor medio en cada voxel — caracterizada por la varianza $\sigma^2$ — afecta directamente las líneas de isodosis y los histogramas dosis-volumen (DVHs). En la simulación ideal, la varianza cero sería el objetivo, pero eso requeriría tiempo infinito de computación.

El rendimiento de un cálculo MC se mide por la eficiencia $\varepsilon$, definida como:

$$\varepsilon = \frac{1}{\sigma^2 T}$$

Donde:

• $\sigma^2$ = estimación de la varianza verdadera de la cantidad de interés
• $T$ = tiempo de CPU necesario para alcanzar esa varianza

Para mejorar la eficiencia, podemos reducir $\sigma^2$ para un dado $T$ (técnicas de reducción de varianza) o reducir $T$ para un dado $N$ sin alterar la varianza.

Geometrías de Scoring: Dosels, Kugels y Órganos Segmentados

No todo scoring necesita ocurrir en voxels cúbicos regulares. El motor PEREGRINE, por ejemplo, utiliza dosels — esferas superpuestas independientes de la malla de transporte de material. La simulación monitorea la desviación estándar en el dosel con mayor dosis y termina cuando alcanza el nivel especificado por el usuario.

Otro enfoque son los kugels, utilizados por el motor Macro Monte Carlo (MMC). Los electrones se transportan en pasos macroscópicos a través de esferas de diferentes materiales y radios, con resultados precalculados almacenados en tablas lookup. El concepto, introducido por Neuenschwander et al. (1995), permite cálculos rápidos de planificación para haces de electrones.

Cuando los voxels individuales se combinan en volúmenes mayores — como órganos segmentados — se pierde resolución espacial, pero se gana en velocidad y precisión en la predicción de dosis por órgano. ¿La limitación? El contorno de volúmenes depende significativamente del clínico que lo realiza. Por eso, el scoring por órgano funciona mejor en ambientes controlados como el phantom NCAT, basado en el dataset Visible Human, que modela anatomía 3D con superficies NURBS e incorpora movimiento respiratorio y cardíaco 4D.

Efectos del Tamaño de Voxel y Varianza Latente

El tamaño del voxel de scoring impacta directamente el tiempo de cálculo, pues determina el tiempo gastado en verificaciones de cruce de fronteras. Cygler et al. (2004) demostraron que el número de historias necesarias por unidad de área es linealmente proporcional a la masa de los voxels de scoring.

Existe también la varianza latente, concepto acuñado por Sempau et al. (2001). Se trata de la incertidumbre residual proveniente de las fluctuaciones estadísticas en el phase space generado por la simulación del cabezal del acelerador. Incluso reutilizando partículas del phase space infinitas veces, la incertidumbre converge a este valor finito e irreducible. Si la varianza latente es significativa en el presupuesto total de incertidumbre, más partículas independientes de phase space deben generarse.

Métodos de Estimación de Incertidumbre

El método batch divide las historias en $N$ lotes (típicamente 10) y estima la incertidumbre como:

$$s_{\bar{X}} = \sqrt{\frac{\sum_{i=1}^{N}(X_i – \bar{X})^2}{N(N-1)}}$$

Walters et al. (2002) identificaron tres problemas: (1) pocos batches generan fluctuaciones en la propia incertidumbre, (2) la agrupación arbitraria ignora correlaciones entre partículas, y (3) el método añade una dimensión extra a las cantidades de scoring.

El método historia-por-historia, implementado según el enfoque de Salvat, resuelve estos problemas manteniendo registros acumulados de $\sum X_i^2$ y $\sum X_i$ durante la simulación, calculando la incertidumbre final vía:

$$s_{\bar{X}} = \sqrt{\frac{1}{N(N-1)} \left[ \sum_{i=1}^{N} X_i^2 – \frac{\left(\sum_{i=1}^{N} X_i\right)^2}{N} \right]}$$

Este método elimina la necesidad de almacenar datos en batches y es computacionalmente más eficiente.

Métodos de Denoising y Suavizado

Reducir fluctuaciones estadísticas simulando más historias es el camino puro, pero impracticable sin técnicas de reducción de varianza poderosas. La alternativa es eliminar el ruido a posteriori con algoritmos de denoising — análogos a la restauración de imagen — acelerando el cálculo MC.

Kawrakow (2002) estableció cinco criterios de precisión que se convirtieron en estándar para evaluar algoritmos de suavizado:

1. Inspección visual de líneas de isodosis (diferencias mínimas respecto al benchmark)
2. Diferencia de área entre DVHs suavizados y de referencia
3. Diferencia máxima de dosis
4. RMSD (root-mean-square difference) entre distribuciones
5. Test x%/y mm — fracción de voxels que difieren más de x% sin punto en el benchmark más cercano que y mm (recomendación actual: 2-3%/2 mm)

Denoising de DVHs

Sempau y Bielajew (2000) propusieron tratar el DVH calculado como la convolución del DVH verdadero (varianza cero) con ruido, aplicando deconvolución para recuperar la señal original. En la práctica, las decisiones basadas en DVHs podrían tomarse rápidamente, viabilizando el uso de MC en todas las etapas de la planificación inversa.

Jiang et al. (2000) adoptaron un enfoque similar, tratando el DVH MC como imagen borrosa y aplicando minimización por mínimos cuadrados iterativos para restaurarlo.

Denoising de Distribuciones 3D de Dosis

El denoising de DVHs no resuelve completamente el problema, ya que las distribuciones de dosis se evalúan de múltiples formas (isodosis, TCP/NTCP). Diversos métodos atacan la distribución 3D directamente:

Filtrado digital (Deasy, 2000) — Primero en proponer denoising 3D, demostró que los filtros digitales mejoran la usabilidad visual y la confiabilidad clínica de distribuciones MC para electrones.

Wavelets (Deasy et al., 2002) — Coeficientes wavelet por debajo de un umbral se anulan, suprimiendo ruido con poca introducción de bias. Aceleración de 2x o más.

Savitzky-Golay 3D (Kawrakow, 2002) — Generalización 3D con ventana adaptativa basada en la incertidumbre local. Reduce historias necesarias por factores de 2 a 20. Extremadamente útil en la fase iterativa de planificación.

Difusión anisotrópica adaptativa (Miao et al., 2003) — Parámetros de filtrado se modifican según el ruido local. Reduce ruido 2-5x (factor de hasta 20x en tiempo de simulación) preservando gradientes de dosis.

IRON (Fippel y Nüsslin, 2003) — Algoritmo iterativo de reducción de ruido que minimiza derivadas parciales segundas de la dosis. Ganancia de 2 a 10x en tiempo de cálculo.

Filtro híbrido mediana-media adaptativo (El Naqa et al., 2005) — Combina operador mediana en regiones de fuerte gradiente con operador media en regiones homogéneas.

Fuente: Monte Carlo Techniques in Radiation Therapy (2nd ed., CRC Press, 2022), Capítulo 11

Conversión CT a Medio: Métodos Estequiométricos

La conversión de números Hounsfield (HU) de la CT en composición material y densidad másica es uno de los eslabones más críticos para la precisión del cálculo MC. Un error aquí compromete todo lo que sigue.

Los primeros algoritmos usaban solo 6 materiales (aire, pulmón, grasa, agua, músculo, hueso). du Plessis et al. (1998) investigaron la influencia de la composición tisular y concluyeron que 7 subconjuntos básicos de tejido eran suficientes para 1% de precisión en dosis con haces de fotones MV. Crearon 57 tipos de tejido variando la densidad física para pulmón y hueso cortical, cubriendo un rango de 3.000 HU.

Conversión Estequiométrica de Schneider

La conversión estequiométrica, propuesta por Schneider et al. (1996) inicialmente para protones, representó un avance significativo. El procedimiento: materiales de composición conocida se escanean por CT para medir los HU correspondientes. Los valores medidos se ajustan mediante la parametrización de Rutherford et al. (1976):

$$\sigma_i(E) = Z_i \left[ K_{KN}(E) + Z_i^{1.86} K_{sca}(E) + Z_i^{3.62} K_{ph}(E) \right]$$

Donde $K_{KN}$ es el coeficiente Klein-Nishina, $K_{sca}$ el coeficiente de dispersión Rayleigh y $K_{ph}$ el coeficiente fotoeléctrico. El coeficiente de atenuación lineal $\mu$ se obtiene vía:

$$\mu(E) = \rho \frac{N_A}{A} \sum_{i=1}^{n} \frac{w_i}{A_i} \sigma_i(E)$$

Schneider et al. (2000) expandieron el método a 71 tejidos humanos, creando 4 secciones en la escala CT donde densidad y pesos elementales se interpolan entre dos tejidos de referencia. La densidad del tejido compuesto es:

$$\rho = \frac{\rho_1 \rho_2}{W_1(\rho_2 – \rho_1) + \rho_1}$$

Donde $W_1$ y $W_2 = 1 – W_1$ representan las proporciones de los dos tejidos de referencia.

Subconjuntos Dosimétricamente Equivalentes (14 bins)

Vanderstraeten et al. (2007) generalizaron el método, creando 14 subconjuntos de tejido dosimétricamente equivalentes (de los cuales 10 eran hueso), garantizando que tejidos en la interfaz entre dos bins no difirieran más de 1% en dosis. La comparación con el esquema convencional de 5 bins mostró diferencias de hasta 5% en las interfaces — particularmente entre tejido adiposo y hueso.

CT Dual-Energy para Conversión Mejorada

La CT de doble energía escanea con dos voltajes diferentes para estimar mejor los números atómicos efectivos $Z$ y las densidades electrónicas relativas $\rho’_e$. Bazalova et al. (2008) encontraron errores medios de extracción de $\rho’_e$ y $Z$ de 1,8% y 2,8%, respectivamente.

Fuente: Monte Carlo Techniques in Radiation Therapy (2nd ed., CRC Press, 2022), Capítulo 11

En la práctica, la CT dual-energy es más ventajosa para haces de baja energía donde la dependencia con $Z$ es mayor.

Registro Deformable de Imagen y Dosis 4D

El movimiento respiratorio en radioterapia de pulmón representa un desafío real para la precisión dosimétrica. El registro deformable de imagen (DIR) permite mapear distribuciones de dosis entre fases respiratorias, viabilizando el cálculo MC 4D.

Keall et al. (2004) realizaron uno de los primeros cálculos MC pseudo-4D: dosis calculada en cada una de las $N=8$ fases de la CT 4D, con posterior mapeo a la fase de referencia vía DIR. Para más detalles sobre entrega dinámica de haz y simulaciones 4D, vea nuestro artículo sobre entrega dinámica de haz y Monte Carlo 4D.

Enfoques de Dose Warping

Heath y Seuntjens (2006) desarrollaron el defDOSXYZnrc, que deforma la grilla de dosis on the fly durante el rastreo de partículas. La desventaja: el tiempo de cálculo aumenta por factor 2 o más.

Siebers y Zhong (2008) propusieron el Método de Transferencia de Energía (ETM), separando transporte de radiación de la deposición de energía. Cuando dos voxels $A_1$ y $A_2$ se fusionan en $X_1$ tras el registro, la interpolación simple calcula:

$$d_{INTERP}(X_1) = \frac{d(A_1) + d(A_2)}{2}$$

Pero lo correcto es considerar la transferencia de energía:

$$d_{ETM}(X_1) = \frac{E(A_1) + E(A_2)}{m(A_1) + m(A_2)}$$

Esta corrección es especialmente importante en regiones de gradiente de dosis, donde la energía por unidad de masa se redistribuye de forma no-lineal.

3D vs. 4D: ¿Cuánto Importa el Movimiento?

Flampouri et al. (2005) estudiaron seis pacientes con IMRT pulmonar y concluyeron que la planificación convencional era suficiente para tumores con movimiento menor a 12 mm. Para movimientos mayores o artefactos severos de CT, el 4D MC se vuelve necesario. Hallazgo destacado: los artefactos de reconstrucción CT tienen mayor efecto que el propio movimiento tumoral.

Seco et al. (2008) compararon 3D FB con 4D FB y 4D gated, introduciendo el índice omega $\Omega(\rho)$ para analizar diferencias de dosis como función de la densidad tisular:

$$\Omega(\rho) = \frac{\int f(\mathbf{x}) \, \delta(\rho(\mathbf{x}) – \rho) \, d\mathbf{x}}{\int \delta(\rho(\mathbf{x}) – \rho) \, d\mathbf{x}}$$

Diferencias de 3-5 Gy entre inhalación y free-breathing se observaron tanto en regiones de baja densidad (pulmón) como de alta densidad (hueso). El MC 4D mejora predicciones no solo en el tumor, sino en los tejidos circundantes que componen el PTV expandido.

En la intercomparación de métodos de dose warping por Heath et al. (2008), los tres métodos probados no mostraron diferencias clínicamente significativas. Sin embargo, cuando el movimiento no era contabilizado, la dosis en el volumen blanco era subestimada en hasta 16%.

Planificación Inversa con Monte Carlo

En la planificación IMRT, la matriz $D_{ij}$ convierte fluencias en distribución de dosis:

$$D(x_i) = \sum_{j \in \text{Beamlet}} D_{ij} \, x_j$$

Donde $D_{ij}$ es la influencia dosimétrica y $x_j$ es la intensidad de cada beamlet. Esta matriz puede calcularse por pencil beam (PB), superposición-convolución o MC. El problema: restricciones de entrega del MLC y dependencias del factor de output con tamaño de campo generalmente no se contabilizan en la generación de los $D_{ij}$, llevando a distribuciones subóptimas.

Jeraj y Keall (1999) fueron pioneros en planificación inversa MC usando MCNP para generar fluencia inicial y EGS4 para calcular dosis en el paciente. Siebers y Kawrakow (2007) propusieron un método híbrido donde la predicción inicial usa PB e iteraciones subsecuentes corrigen con MC, obteniendo ganancia de 2,5x sobre optimización MC pura.

Enfoque Broad Beam (Nohadani et al., 2009)

La innovación de Nohadani et al. fue usar un phase space abierto tras las mandíbulas (que no cambian durante la entrega IMRT en linacs Varian) dividido en píxeles PB-like. Tres ventajas claras: (1) contabiliza heterogeneidades en el paciente, (2) el factor de output es el de campo abierto estándar sin correcciones adicionales, y (3) variaciones de output por cambio de campo durante entrega por MLC son pequeñas.

La comparación entre DVHs de planes IMRT optimizados con MC-PB versus PB convencional reveló algo importante: el plan basado en PB aparenta ser mejor en términos de cobertura del PTV y sparing de órganos en riesgo — pero se desvía significativamente de la dosis real depositada. Los errores de cálculo de dosis no pueden negligirse en la planificación clínica.

Aplicaciones Clínicas en Haces de Electrones

Para haces de electrones, el MC demuestra ventajas claras sobre algoritmos analíticos como pencil beam. En heterogeneidades (hueso, aire, pulmón), el pencil beam falla en predecir dispersión lateral, backscatter y efectos de cavidades — situaciones donde el MC sobresale al rastrear cada partícula a través de todas las interacciones físicas relevantes.

El commissioning de un TPS basado en MC para electrones requiere atención a: verificación de cálculos con datos medidos, normalización y prescripción de dosis, incertidumbre estadística versus tamaño de voxel y la distinción entre dosis-al-medio y dosis-al-agua. Conozca más sobre Monte Carlo para electrones y braquiterapia en nuestro artículo dedicado.

La evolución de los TPSs comerciales de MC para electrones trajo tiempos de cálculo clínicamente aceptables — aproximadamente 4 minutos para cálculos con incertidumbre estadística de 2% y voxel de 0,25 cm.

Consideraciones Finales

El cálculo de dosis Monte Carlo en paciente evolucionó de herramienta académica a componente viable de sistemas de planificación clínica. Los avances en denoising (reducciones de 2-20x en tiempo de cálculo), conversión estequiométrica CT (de 5 bins a 14 bins dosimétricamente equivalentes) y planificación inversa MC posibilitan cálculos más precisos sin sacrificar la practicidad.

Para aplicaciones en electrones, donde las heterogeneidades tisulares desafían algoritmos convencionales, el MC ofrece la solución definitiva — siempre que sea acompañado de commissioning riguroso. Con la creciente integración de registros deformables 4D, el MC camina hacia capturar no solo la anatomía estática, sino toda la dinámica respiratoria del paciente.

Para profundizar en los fundamentos del Monte Carlo o explorar cómo el método se aplica a fotones en aplicaciones clínicas, consulte los artículos dedicados de nuestra serie.

Fuente principal: Monte Carlo Techniques in Radiation Therapy (2nd ed., CRC Press, 2022), Capítulos 11-12.