Algoritmos de cálculo de dosis para haces de
fotones.
SAFIM 2014. Buenos Aires, Argentina
Armando Alaminos Bouza
MEVIS Informática Médica LTDA., São Paulo, Brasil
www.mevis.com.br
Formulación determinista de la distribución de dosis :
Ecuación del Transporte de Boltzmann
Resolver la ecuación del transporte directamente era un problema
imposible de enfrentar hasta hace poco tiempo
(Ecuación 2 es la Boltzmann Fokker-Plack, que trata transporte de electrons)
¿ Como solucionar la formulación de Boltzmann del transporte ?
- Método estocástico : Monte Carlo (MC).
Resuelve el problema en forma indirecta al seguir historias de un gran
número de particulas generadas aleatoriamente.
- Método determinista : Grid Based Boltzmann Solver (GBBS ).
Método de solución directa que discretiza las variables, volume,
ángulos, energia.
Dato curioso : ¡ Ambos métodos nacieron en “Los Alamos National
Laboratory” !
El desarrollo de los GBBS fue esfuerzo para encontrar alternativas más
rápidas que la solución con MC.
¿ Como solucionar la formulación de Boltzmann del transporte ?
- Método estocástico : Monte Carlo (MC). Muy lento para
fines clínicos o demanda un poder de cálculo con costo
muy alto para radioterapia. Esto ha comenzado a cambiar
en los últimos años. Pero aún así, es muy lento.
- Método determinista : No se conocia un método de solución
con resultados suficientemente realista.
¿ Como resolver el problema ? Nacen diversos métodos
aproximados.
Su precision ha mejorado en la medida que aumentó el poder
de cómputo del hardware accessible.
Si observamos la historia de este proceso en los últimos 40 años
vemos una evolución de métodos basados en
“datos/correciones” a métodos basados en “modelos físicos”.
Los métodos basados en “datos medidos y correciones” no los vamos a
discutir hoy. Aún tienen valor para algunas aplicaciones, pero han perdido
mucho terreno ante los métodos basados en modelos físicos, además los test
de aceptación actuales tienen tendencia a fallar con ellos.
Dos familias de algoritmos son muy frecuentes en los sistemas de
planificación de radioterapia modernos (2004-2014) :
- Integración de Pencil Beams (¿ expresión en Español ? )
- Integración de kernels de depósito de energía, conocidos como “Energy
Deposition Kernels” ( EDK)
Sea un haz de fotones monoenergéticos.
En principio, si conocemos la distribución del flujo de fotones a la entrada de un medio y conocemos el coeficiente de atenuación para todo el
volume del medio podemos calcular la densidad del flujo de fotones en
cualquier lugar del medio y además la cantidad de energía que los
fotones liberan en cualquier punto ( r ). La energía total liberada se
identifica como TERMA (total energy released to unit mass – Anhesjo
1987).
Donde: TE(r) es el TERMA en r, µ(E,r) es el coeficiente de
atenuación lineal para la energía E en el punto r, ψE(r0)
densidad de flujo de energía en el plano r0 .
(Ahnesjö 1989)
Ahora nos bastaría conocer como la energía liberada al medio por los
fotones (el TERMA) se deposita alrededor del punto de liberación.
Si resolvemos este problema tendremos un algoritmo muy robusto y flexible.
Varios autores publicaron al respecto en los años 80 :
- Dean 1980
- Mackie 1985
- Boyer y Mok 1985
- Mohan y Chui 1986 - Differential Pencil Beam
- Ahnesjö 1987 y 1898.
Primer problema a resolver : Crear los “mapas” de distribución de la dosis alrededor del punto de liberación (o interacción primaria).
EDK : IntroducciónEl EDK es un mapa de la distribución de energia alrededor de un punto
de interacción primaria de uma partícula (en este caso fóton).
Los EDK son casi imposibles de generar analíticamente o en forma
experimental. En la práctica siempre se emplea un código de MC para
generar una biblioteca de EDKs
Los EDK contienen fracción de energia depositada por unidad de
masa, normalmente normalizado con la energia entregada en el
origen. Generalmente la geometria empleada es una esfera (agua) de
600 mm de radio, dividida en casquetes concéntricos y conos.
Energy Deposition Kernel para fotones monocromáticos en agua.
Cubo de 1 metro de lado. Curvas de 1E-5, 1E-6, 1E-7, 1E-8, 1E-9.
1 Mev 4.0 MeV 15.0 MeV
EDK creados con EGSnrc, modificado para 64 esferas y 90 conos. Secciones de interacción
según NIST XCOM, ECUT=EA=0.512 MeV y PCUT=PA = 0.001 keV, mínimo número de historias
iniciales entre 4.0E6..,8.0E6 para tener error < 1 % em región de dispersión primaria.
(desarrollado con colaboración de Dr. Ernesto Mainegra, NRC, Canada)
dosis por electrones primarios dosis total
Mapas de depósito de energía en EDK de 4 MeV. Cubo de agua de 1 metro de lado
Curvas de 1E-5, 1E-6, 1E-7, 5.0E-5, 1E-8, 1E-9.
EDKnrc, NIST XCOM, ECUT=EA=0.512 MeV y PCUT=PA = 0.001 keV
(desarrollado con colaboración de Dr. Ernesto Mainegra, NRC, Canada)
En principio, puede implementarse la convolución de los EDKs y el TERMA
en forma directa, integrando la contribución de todos los voxels dispersores
en forma explícita, pero este abordaje produce tiempos de cálculo
inaceptables para todos los fines prácticos. Además, los EDK presentan una singularidad en r = 0, creando inestabiblidad en la integración
numérica.
Una solución empleada inicialmente fue la utilización del teorema de la convolución que permite resolver la convolución como el producto punto
de las transformadas de Fourier del EDK y el TERMA.
Pero la solución via teorema de la convolución solo es válida para el caso
de un EDK invariante para todo el volumen, esta restricción impone serias
limitaciones. En primer lugar, la invariabilidad de los EDKs entra en
contradicción con los medios heterogéneos, además de no modelar
correctamente la inclinación de los EDKs en haces divergentes ni el endurecimiento del espectro con la profundidad.
La solución aproximada de la convolución del TERMA y los EDKs que menos
restricciones impone fue publicada en 1989 por Anders Ahnesjo.
En lugar de integrar la contribución de todos los voxels dispersores al punto
de cálculo se propone dividir el volumen dispersor en un número finito de
conos.
En cada cono integramos solamente considerando los valores de TERMA y
propiedades del tejido presente únicamente en el eje del cono.
La precisión de esta solución depende de la cantidad de conos
utilizados, pues algunas heterogeneidades pequeñas pueden ser ignoradas
cuando ningún eje de cono las intercepta. Afortunadamente la densidad
de ejes de conos en regiones próximas al punto central de depósito de dosis
es alta resultando menos probable que se omita una heterogeneidad de
importancia clínica.
El método descrito por Ahnesjo se denomina Convolución con Conos Colapsados - collapsed cones convolution (CCC)
EDK polienergético analítico de Ahnesjö, publicado en 1989 :
Como se espera, contiene una singularidad en el origen ( r = 0), pero su
Integral en una esfera desde el origen al infinito, es finita, al igual que la
integral en un cono. El kernel para un cono de ángulo solido Ω en la
dirección θm , Φn :
¡ Con esto eliminamos la singularidad ! Además, disminuimos la
complejidad de la integración numérica de O(N7) para O(N3 x M).
Donde N3 son los voxels presentes y M los conos empleados para
cada kernel.
X =
Curvas de isoTERMA Curvas de isodosis
La expresión analítica de EDK de Ahnesjö es usada em implementaciones
convensionales de CC, pero la necesidad de esta parametrización puede
resultar muy limitante para modelar cualquier tipo de haz (Lu, Mackie 2005).
Se puede implementar el método CC accesando diretamente las tablas de
EDK generadas com MC.
Métodos para implementar el CC empleando las tablas deben ser estables
frente a la fuerte variación del EDK diferencial dentro de la dimensión de un
voxel para localizaciones próximas al origen. Hoy se prefiere crear kernels pre
integrados que se denominan “cumulative” o “cumulative-cumulative” kernel.
De este modo se logra reproducibilidad de los resultados para voxels de
diferentes tamaños.
6 MeV , CCC-CAT3D vs EGSnrc 1.25 MeV, CCC-CAT3D vs EGSnrc
Para los haces polienergéticos se puede pensar en dos soluciones :
1. Método de components. Sería la implamentación más rigurosa, pero
genera una convolución para cada component del espectro tornandolo
muy lento.
2. Método de EDK polienergético. Creado como suma ponderada de EDK
monoenergéticos, considerando fracción en el espectro y coeficiente de
atenuación lineal de cada componente.
Espectro de 1.25, 6.0 MeV, CCC-CAT3D (método 2) vs EGSnrc
¿ Como se tratan las heterogeneidades ?
• El TERMA se calcula con el camino
radiológico efectivo entre la Fuente y el
punto de interacción. El TERMA se procesa
energía por energía del espectro, esto es
muy preciso.
• El EDK se escala según el camino radiológico
efectivo entre el voxel emisor y el colector de
energía. Casi todos los TPS consideran que el
EDK se escala según la densidad electrónica
relativa al agua, esto resulta muy bueno para
pulmón pero es menos preciso para huesos y
protesis de metales.
EDK de 1.5 MeV para angulos de 0 a
3.75 y 86.5 a 90 grados. Para agua,
hueso y titanio. Huang et al. 2013
¿ Como implementar la convolución ?
a) Punto de vista del emisor de energía. Usado por Helax y Oncentra
b) Punto de vista del colector de energía. Usado por Pinnacle, Xio, CAT3D.
La variante con el punto de vista del colector es muy
atractiva en las arquitecturas que permiten procesamiento
paralelo como las CPU multinúcleos, y el procesamiento
con GPU o Xeon Phi.
A juzgar por los resultado de pruebas empleando de 1 a 12
hilos (threads) de procesamiento, estimo que con
multithread masivo, como en un Xeon Phi, los tiempos de
convergencia para el CCC puedan reducirse
significativamente resultando en respuestas casi
instantáneas.
Solución con Pencil Beam.
• Kernel puede obtenerse por MC o por
deconvolución de resultados experimentales.
• Solución con FFF impone restricciones.
• Integración o superposición es más flexible.
• Corrección de heterogeneidades por camino
radiológico efectivo longitudinal.
• Correción por heterogeneidad lateral escalando
el kernel.
Positivo :
• Rápido
• Modelo adecuado para optimizar IMRT
Negativo :
• Falla en interfase de heterogeneidades.
• Ignora el “build-down” en salida del haz.
• Tratamiento de heterogeneidades es incompleto.Generación de PB con MC
- Método determinista : Grid Based Boltzmann Solver (GBBS ).
Attila fue el primer GBBS
comercial de propósito
general. Desarrollado por
Transpire INC.
El Acuros es un código
especializado para RT,
desarrollado a partir del Attila,
pero aproximadamente 10
veces más rápido.
Desarrollado por Transpire
INC.
Comparación de resultados. MC , Acuros , AAA, CCC
Aire Pulmón Hueso
Journal Med.
Phys. Clin
Engin.Radit.
Oncol. 2012
Dosis a la salida del haz (“exit dose”)
Yub I. Tan, et al 2014. Detalle de “exit dose” CCC / CAT3D
Como era de esperar, PB no modela
Falta de retro dispersion a la salida.
Tan et al.: Evaluation of TPS algorithms , 2014
Gamma en dosis absoluta 3%/3mm, campo 20 cm2, diferentes profundidades
Conclusiones :
• La accesibilidad a poder de cálculo creciente traerá a los TPSs
algoritmos cada vez más complejos y próximos a la realidade física.
• Disponer de varios algoritmos para planificar RT es deseable y el físico
médico debe ser capaz de elegir el más adecuado para cada caso.
• Todos los programas tienen algunos “bugs”, no hay excepción. Portanto, siempre emplee un método de cálculo, independiente de su
TPS principal, para para verificar Unidades Monitoras, antes de tratar al
paciente.
¡ Muchas gracias !
(ver referencias bibliograficas en próximas láminas)
Lecturas recomendadas:
Generación de EDK :• “Generation of photon energy deposition kernels using the EGS Monte Carlo Code”, Mackie T.R,
Bielajew A.F, Rogers D.W.O, Battista J.J. Phys. Med. Biol, 1988, vol 33, No. 1.
• “Calculation of photon energy deposition kernels and electron dose point kernels in water”, Mainegra-
Hing E, Rogers D.W.O., Kawrakov I., Med. Phys. 32, 685 (2005).
• “Investigation of various energy deposition kernel refinements for the convolution/superposition
method”, Huang Jessie, Eklund D, Childress N., Howell R., Mirkovic D., Followill D, Kry S., Med. Phys. 40,
Dec. 2013.
Trabajo original sobre “Collapsed Cone Convolution”:• “Collapsed cone convolution of radiant energy for photons dose calculation in heterogeneous media”,
Anders Ahnesjö. Med. Phys. 16 (4) , Jul/Aug 1989.
Como implementar CCC a partir de EDK tabulados, directos del MC.
• “Accurate convolution/superposition for multi-resolution dose calculation using cumulative tabulated
kernels”, Lu W., Oliveira G., Chen M., Reckwerdt P., Mackie T., Phys. Med. Biol., 50 (2005) 655-680.
Resultado experimentales presentados : • “Accuracy of the Small Field Dosimetry Using the AcurosXB Dose Calculation Algorithm within and
beyond Heterogeneous Media for 6 MV Photon Beams”. S. Stathakis, C.Esquivel, L.Quino, P.Myers,
O.Calvo, P.Mavroidis, A.Gutiérrez, N.Papanikolaou., International Journal of Medical Physics,Clinical
Engineering and Radiation Oncology, 2012, 1, 78-87.
Resultado experimentales presentados :
• “Evaluation of six TPS algorithms in computing entrance and exit doses”. Y. Tan, M.
Metwaly, M Glegg, S. Baggarley, A Elliott., Journal of Applied Clinical Medical
Physics, Vol. 15, no. 3, 2014.
Sobre el Acuros :
• “Validation of a new grid-based Boltzmann equation solver for dose calculation in
radiotherapy with photon beams”, O. Vassiliev, T. Wareing, J McGhee, G. Failla, M.
Salehpour, F. Mourtada., Phys. Med. Biol., 55, 581 (2010).
Top Related