Introducción
Desde su descubrimiento, en 1895, los rayos X han jugado un papel importante en la ciencia moderna, especialmente en el diagnóstico y tratamiento de varias enfermedades en medicina y odontología, lo que a su vez ha repercutido significativamente en el mejoramiento de la salud de muchas personas alrededor del mundo (Pernicka & McLean, 2007).
Rayos X terapéuticos de kilovoltaje (kV) son rutinariamente empleados para el tratamiento curativo de algunos tumores, tales como el carcinoma basocelular y el escamo-celular de piel, bien como para el tratamiento paliativo de metástasis óseas (Bilge, 2004).
El espectro de energía es una parte esencial para la caracterización completa de un equipo de rayos X (Malezan, et al., 2015), puesto que tiene un gran impacto sobre la calidad y la dosis entregada por el haz de rayos X (Nickoloff & Berman, 1993; Pamplona & Costa, 2010; Abbene, et al., 2012). De aquí que el espectro y la calidad del haz de rayos X sean parámetros imprescindibles para el estudio de las propiedades dosimétricas en radiodiagnóstico (Chen, et al., 2012) y radioterapia (De la Vega, et al., 2008).
Se define al espectro de rayos X como la distribución energética de la radiación producida por un equipo de rayos X (Nickoloff & Berman, 1993). El espectro de rayos X, además, suele ser requerido para el modelamiento matemático y la optimización de sistema de imágenes en radiodiagnóstico (Mainardi & Bonzi, 2008) y radioterapia (Li, et al., 2008).
El espectro de energía de rayos X se puede medir de manera directa utilizando un espectrofotómetro (Mainardi & Bonzi, 2008, Chen, et al., 2012). Con todo, dado que la medición espectrométrica del haz de rayos X es demasiado costosa, exige una alta cualificación y requiere mucho tiempo, su aplicación en la rutina clínica no es común (Baird, 1981; Pamplona & Costa, 2010; Chen, et al., 2012).
Como alternativa, varios autores han reportado métodos indirectos de determinación del espectro de energía de un haz de rayos X a partir de la curva de atenuación del haz (Dossel & Schlegel, 2010) y del uso de transformadas de Laplace (Delgado, 1999; Delgado, 2007; Mainardi & Bonzi, 2008).
En su artículo pionero, Silberstein fue el primero en proponer el uso de la transformada de Laplace para la reconstrucción de espectros de rayos X a partir de una curva de atenuación (Silberstein, 1993). Poco más tarde, en 1936, Bell testeó la eficiencia del método de Silberstein usando curvas experimentales (Bell, 1936). Décadas después, en 1970, Twidell desarrolló un código computacional para reconstruir espectros de rayos X de manera iterativa, por minimización del error entre datos de atenuación experimentales y numéricos (Twidell, 1970). En la década de los 80, fueron publicados varios trabajos para la reconstrucción de espectros de rayos X a partir de datos de atenuación, utilizando modelos matemáticos de tres o más parámetros, basados en diferentes métodos de solución (Baird, 1981; Archer & Wagner, 1982; Kramer & Von Seggern, 1983; Kramer, et al., 1983; Archer, et al., 1985; Archer & Wagner, 1988a, 1988b).
En la última década, diferentes estudios han explorado los aspectos experimentales y numéricos de la reconstrucción del espectro de rayos X a partir de la aplicación de distintas metodologías (Delgado, 2009; Glover & Chantler, 2009; Menin, et al., 2014; Malezan, et al., 2015; Panthi, 2018). Entre los trabajos que usaron el método de mínimos cuadrados para determinar el espectro de rayos X, se destaca el de Archer & Wagner, 1988b. En ese estudio, se plantea una modificación al modelo analítico de atenuación para remover artefactos relacionados con soluciones no realísticas del espectro. No obstante, no queda claro si la introducción de esa modificación era necesaria para suprimir los artefactos generados o si más bien bastaba con imponer restricciones al método de mínimos cuadrados.
Por otra parte, también ha sido reportada la derivación de espectros de energía de rayos X a partir del uso de códigos Monte Carlo (Ay, et al., 2004; Bonifácio, et al., 2005) y programas computacionales tales como IPEM78 (Ay, et al., 2004) y SpekCalc (Poludniowski, et al., 2009).
La capa semirreductora (CSR) es un parámetro comúnmente utilizado para describir la calidad de un haz de rayos X (Nickoloff & Berman, 1993; ChenAy, et al., 2012; Al, 2017). La capa semirreductora se define como la espesura de un material absorbedor que reduce a la mitad la intensidad original del haz (Abbene, et al., 2012). En haces de kV, la energía efectiva puede ser derivada a partir de la capa semirreductora, convirtiendo esta última a coeficiente de atenuación linear o másico (Dossel & Schlegel, 2010; Chen, et al., 2012). La energía efectiva es definida como la energía de un haz polienergético con igual capa semirreductora que un haz monoenergético (Chen, et al, 2012, Correa, et al., 2012). Luego, la energía efectiva también puede ser usada para evaluar la calidad de un haz de rayos X (Correa, et al., 2012).
En el presente artículo, se propone reconstruir el espectro de energía de un haz de rayos X terapéutico, aplicando el método indirecto propuesto por Archer & Wagner, 1988a, fundamentado en la utilización de un par de transformadas de Laplace y en la medición experimental de la curva de atenuación del haz. Para resolver este método se utilizará, como novedad, el algoritmo multistart junto con la función lsqnonlin de MATLAB. La validación de la metodología de reconstrucción será realizada por comparación de los valores de la capa semirreductora de la curva de atenuación y del espectro, en lugar de hacerla cotejando el espectro reconstruido con el real, como hecho en artículos anteriores.
Materiales y Métodos
Modelamiento Matemático
Para un mismo medio material, bajo la condición de equilibrio de partículas cargadas (EPC), la dosis absorbida es equivalente al kerma colisional (Attix, 2008; Ma, 2014). Por tanto, si el medio material es el aire, se tiene que (Bos, 2011),
donde D ar es la dosis absorbida en aire, K col,ar es el kerma colisional en aire, E max es la energía máxima de los fotones del haz de energía nominal E, Φ(E), es la distribución energética de fluencia, es el coeficiente másico de absorción de energía y p es la densidad del material absorbedor.
El kerma colisional en aire también se puede escribir como (Smith, 2000),
Kar= MN k (2)
en que M es la lectura de la cámara de ionización corregida por todas las variables de influencia y N es el coeficiente de calibración en términos de kerma en aire. Reordenando los términos, a partir de las Ec. (1) y Ec. (2), se tiene que
Si un haz estrecho es atenuado por un material absorbedor de espesor, x, la distribución energética de fluencia viene a ser,
donde µ m es el coeficiente másico de atenuación del material absorbedor.
De esta manera, la lectura de la cámara de ionización como función del espesor del material, M(x), es,
donde el kerma en aire del material es, justamente,
La transmisión relativa del haz a través del material absorbedor de espesor, T(x), se define como la razón entre el kerma total del haz atenuado en aire, Kar (x), y el kerma total del haz no atenuado en aire, Kar (0),
Siendo lo mismo que,
Por otra parte, sea,
Entonces, a partir de Ec. (5), Ec. (8) y Ec. (9), la transmisión relativa del haz se torna,
en que F(E) dE es la fracción de M(0) debida a los fotones de energía entre E y E + dE.
Cambiando la variable de integración de dE para dµ m , en la Ec. (9), resulta (Archer & Wagner, 1988a)
Para simplificar, considere que,
donde f (µm) es el prespectro relacionado con el espectro de energía del haz, F(E) (Archer & Wagner, 1988a),
Luego, a partir de Ec. (11) y Ec. (12), se obtiene,
en que si µ m es una función diferenciable y descendiente con respecto a E, siempre que E ⟶ 0, entonces µ m ⟶ ∞, o, si E ⟶ ∞, entonces µ m ⟶ 0 (Pamplona & Costa, 2010).
Como la Ec. (13) concuerda con la definición de transformada de Laplace (Dodson, 2002), la curva de transmisión de un haz de rayos X puede ser rescrita así,
T(x) = L [f (μm)]. (14)
La Ec. (14) es la representación matemática exacta del proceso de atenuación espectral. Si la curva de transmisión es conocida, la transformada inversa de Laplace permite derivar el espectro de energía del haz de rayos X, F(E), tal que,
o también,
En este artículo se estudia un haz de rayos X de 80 kV producido en un ánodo de tungsteno. Para este material, la energía total de los electrones de la capa K es -69.5 keV, por lo que la emisión de rayos X característicos solo ocurriría para un valor de voltaje de tubo superior a este valor. La mayor parte de los rayos X emitidos son de frenado por causa de que aumentan fuertemente a medida que se aumenta del voltaje aplicado, lo que no sucede con los rayos X característicos. Sin embargo, la cantidad de rayos X característicos se incrementa para electrones con energía por encima de la energía de enlace de la capa K. La estabilidad de este proceso se alcanza para energías cinéticas más elevadas (Dyson & Dyson, 2005). La radiación característica de la capa L, no es generalmente importante y no será tenida en cuenta en este trabajo.
Para valores de voltaje de tubo superiores a la energía de enlace de la capa K, la curva de transmisión puede ser expresada como la superposición de su componente de frenado (espectro continuo), Tf (x) y su componente característica (espectro discreto), T c (x), esto es (Archer & Wagner, 1988a),
T(x) = Tf (x) + Te (x). (17)
Si r es la fracción de kerma en aire debida a la radiación de frenado en el espectro no atenuado, consecuentemente 1 - r será la fracción de kerma en aire debida a la radiación característica. Con esto, los términos TF (x) y Te (x) son dados por,
Donde Ff (µm) es la fracción de frenado del espectro de energía y,
donde C¡ es la abundancia relativa del i-ésimo haz de rayos X característicos y µm,i es el coeficiente másico de atenuación del i-ésimo haz de rayos X existente para un haz de energía nominal E.
Usando la propiedad asociativa de la transformada inversa de Laplace, de las Ec. (14) y Ec. (17), se tendría,
El significado físico de esta formulación está en que el último término es una suma de las funciones delta de Dirac, ô, que representa a los rayos X característicos, ya que
Usando la identidad,
El espectro de la radiación característica, F C (E), se puede escribir como,
con lo que el espectro de energía del haz vendría a ser,
F(E) = rF F (E) + (1 - r)F c (E), (24)
donde Fe(E) será entonces el espectro de energía de los rayos característicos. Las Ec. (17) y Ec. (24) constituyen la descripción matemática del espectro de rayos X teniendo en cuenta sus componentes de radiación de frenado y de radiación característica. (Archer & Wagner, 1988).
En este trabajo, se asume que la curva de transmisión viene dada por el modelo propuesto por Archer & Wagner (1988a), es decir,
donde r,a,b, ϑ son parámetros de ajuste a ser determinados y es el coeficiente másico de atenuación del absorbedor a la energía máxima del haz.
Asimismo, el espectro de energía del haz de rayos X, F(E), viene expresado matemáticamente como (Archer & Wagner, 1988a),
donde Γ (ϑ) es la función gamma y Γ ϑ -1/2 (1/2 (a - b)(µm - µm,0)) es una función de Bessel modificada.
Metodología de reconstrucción espectral
La metodología seguida para el cálculo del espectro de rayos X es la siguiente:
1. Cálculo de la curva de transmisión. En primer lugar, se toma la lectura de la cámara de ionización sin material absorbedor. Seguidamente, se interponen placas de aluminio de 1 mm de espesor, progresivamente, de forma que la espesura va aumentando cada vez que se toma una nueva lectura. Se hace uso de la Ec. (8) para encontrar la curva de transmisión.
Resumiendo: la primera lectura, al no tener placa de aluminio (x = 0 mm), se correspondería con T(0); la segunda, para x = 1 mm (1 placa), con T(1); la tercera, para x = 2*1 mm = 2 mm (2 placas), con T(2); la cuarta, para x = 3*1mm = 3 mm (3 placas), con T(3) y así sucesivamente hasta llegar a x = 5 mm, con T(5). Luego, la curva T(x) se expresa en función de x, con x = 1,2, ...,5 mm. Como regularmente el espesor es indicado en términos de la densidad del absorbedor, x se convierte en x = px.
2. Cálculo los parámetros de ajuste r, a, b, ϑ de la Ec. (17). Para reproducir el coeficiente másico de atenuación de un haz polienergético, un polinomio de quinto grado es utilizado. Hecho esto, se procede a hallar (Archer, et al., 1985). Las constantes del polinomio fueron determinadas utilizando la función fit de MATLAB. Para ánodos de tungsteno, los valores de la abundancia relativa y de los coeficientes másicos de atenuación para rayos X característico se consultaron en la tabla 1 de Archer & Wagner (1988). Alternativamente, los valores de los coeficientes másicos de atenuación se pueden encontrar en las tablas del NIST (National Institute of Standards and Technology, en inglés).
Una vez establecidas las abundancias relativas y los coeficientes másicos de atenuación de los fotones de frenado y característicos en aluminio, se emplea algún método de optimización conducente a minimizar el funcional, ξ(T),
donde T(x) es obtenida experimentalmente a través de la Ec. (8) y los datos del numeral 1. Como tentativa inicial de solución se tomaron valores aleatorios, para cada parámetro, usando la función rand de MATLAB. Las restricciones impuestas a los parámetros fueron a > 0; b > 0; ϑ > 0 y 0 < r ≤ 1.
3. Reconstrucción del espectro de rayos X. Después de que los parámetros de ajuste de la Ec. (27) han sido descifrados, se substituyen en la Ec. (18) para así dar con el espectro. La minimización del funcional, ξ(T), fue llevada a cabo mediante el algoritmo multistart para la función lsqnonlin, perteneciente a la caja de herramientas de optimización de MATLAB. El número de puntos de partida del algoritmo multistart fue establecido en 50.
La validación de la reconstrucción espectral se realizó a partir de la comparación del valor de la capa semirreductora obtenida a partir del espectro de rayos X reconstruido con el de la capa semirreductora extraída de la curva de transmisión (Santos, et al., 2016). Esta forma de validación es más simple, rápida y práctica que aquella basada en la comparación del espectro reconstruido y el real, empleada en trabajos anteriores (Archer & Wagner, 1988a, 1988b), dado que no requiere del conocimiento del espectro real.
Asimismo, la medición del espectro real es un proceso complejo y problemático ya que: i) necesita disponer de un equipo especializado y altamente costoso como es el espectrómetro Compton; ii) exige preparar y/o contar con un operador experimentado en el manejo de tal equipo; iii) la caracterización del equipo y la medición del espectro toma mucho tiempo y iv) no es clínicamente práctico pues tocaría utilizar por mucho tiempo, a fin de realizar las mediciones, los equipos dedicados al tratamiento de los pacientes.
La solución ideal al problema de conseguir el espectro real sería que el fabricante lo proporcionara junto con el equipo o de alguna otra forma. Empero, esta información no suele ser suministrada con facilidad, incluso cuándo solicitada.
La capa semirreductora del espectro se calcula usando la ecuación,
donde K ar0 , K ar1 , y K ar2 , son los valores de kerma en aire medidos con la cámara de ionización sin material absorbedor (x 0 ) y para valores de espesor del material absorbedor de x1 (medida experimental del espesor inmediatamente menor a la CSR obtenida por la curva de transmisión) y x 2 (medida experimental del espesor inmediatamente mayor a la CSR obtenida por la curva de transmisión), respectivamente (Santos, et al., 2016).
De la discretization de la Ec. (6) y de la Ec. (18), se encuentra el valor del kerma en aire en función del espesor del absorbedor, o sea,
donde i = [1, 2,...n] es el i-ésimo haz de energia E i , que conforma el espectro de rayos X, n es el índice del haz de máxima energía de los fotones que conforma el espectro de rayos X; j = [0,1, 2, ...5] es el j-ésimo valor del espesor del material absorbedory ΔE ¡ es el intervalo de energía escogido.
Para este trabajo, se utilizó un equipo de ortovoltaje Siemens Stabilipan II, perteneciente al Hospital das Clínicas de Ribeirão Preto (SP), Brasil, operando a 80 kV y 20 mA. Una cámara de ionización Radcal 2086 fue empleada para realizar las lecturas. Las lecturas fueron registradas a una distancia fuente-centro de la cámara de ionización de 40 cm.
Los cálculos involucrados en el proceso de reconstrucción fueron realizados con MATLAB® 2015a (Microsoft Windows 10 Home, CPU: 2,3 GHz, RAM: 4 Gb). El tiempo de cálculo del espectro fue de 8 s, aproximadamente.
Resultados y discusión
En la figura 1 se presenta la curva de transmisión normalizada del haz de rayos X terapéutico de 80 kV El patrón de la curva de transmisión es el esperado de la literatura, esto es, un decaimiento exponencial progresivo que parte desde el valor de transmisión sin material absorbedor.
Los valores de los parámetros de ajuste de la ecuación (27) encontrados en el proceso de minimización, indispensables para generar el espectro de energía del haz de rayos X, son mostrados en la tabla 1.
Energía nominal del haz, E = 80 kV; distancia fuente-centro de cámara de ionización de 40 cm.
La figura 2 muestra el espectro de energía de rayos X reconstruido para un haz de energía nominal de 80 kV Su apariencia coincide con los datos de la literatura, es decir, una zona de corta extensión y baja energía seguida por un pico ancho que concentra la mayoría de los fotones y acabando con una zona que decae exponencialmente, con líneas verticales propias de los rayos X característicos, hasta la energía nominal del haz. La energía pico del espectro reconstruido se localizó en 29,5 kV.
La capa semirreductora del espectro, como mencionado en la sección anterior, se calcula mediante las ecuaciones (28) y (29). Los valores de x 1 y x 2 son extraídos de las medidas experimentales de la curva de transmisión de la figura 1, siendo estos x 1 = 0,58 g/cm2 y x 2 = 0,74 g/cm2. Una vez hecho todos los cálculos, considerando ΔE ¡ = 2 kV, se obtuvo el valor de la capa semirreductora del espectro.
Entretanto, para calcular el valor de la capa semirreductora de la curva de transmisión, se evalúa el valor del espesor de aluminio para el cual T=0,5.
Los valores de la capa semirreductora obtenida a partir de la curva de transmisión y del espectro reconstruido se encuentran consignados en la tabla 2.
La validación del método, realizada comparando el valor de CSR F(E) con el de CSR T(x) vía el error relativo porcentual entre CSR F(E) y CSR T(x) , resultó en 2,65%.
Dado el pequeño error relativo porcentual resultante y el corto tiempo de cálculo del espectro, se colige que el procedimiento empleado para la reconstrucción del espectro de energía de un haz de rayos X terapéutico a partir de su curva de transmisión y la sinergia multistart-lsqnonlin mostró ser eficaz y rápido.
Adicionalmente, vale la pena hacer mención de que este método de reconstrucción espectral no presentó, en ningún caso, soluciones no realistas del espectro como las reportadas por Archer & Wagner, 1988b, debido a que el valor de ϑ, siempre fue mayor que 0,6. Esto relieva el buen funcionamiento del algoritmo escogido para solucionar el problema.
Por otro lado, como no fue posible comparar el espectro reconstruido con el espectro real del equipo, el espectro reconstruido en este trabajo es un "espectro estimado" (Baird, 1981) que, por los resultados del método de validación, posee el mismo efecto clínico que el real, pero del cual no se tiene certeza si es próximo de este. Luego, una perspectiva futura de este trabajo sería hacer la validación de la reconstrucción tanto por comparación de CSR, como por comparación de espectros para, así, evaluar que tanto el espectro reconstruido se aproxima del espectro real del equipo.
Conclusión
El método de reconstrucción espectral de un haz de rayos X a partir de su curva de atenuación y de la sinergia multistart-lsqnonlin de MATLAB mostró resultados satisfactorios en términos de eficacia y rapidez. Así, este método podría ser aplicado con confianza para calcular magnitudes dosimétricas o radiométricas a partir del espectro reconstruido.
El análisis de la curva de atenuación en comparación con los métodos espectroscópicos es una buena alternativa para caracterizar el haz, ya que estas curvas pueden ser obtenidas de manera simple en relación con la medición directa del espectro. Para esto, solo basta disponer de un equipo de rayos X y un conjunto electrómetro-cámara de ionización, elementos de rutina en cualquier servicio de radioterapia.
Para generalizar el uso de este método de reconstrucción se propone, como perspectiva futura, considerar un análisis que abarque haces de rayos X de diferente energía nominal y/o algún método de optimización híbrido que resulte en una mejor aproximación entre las capas semirreductoras asociadas a los espectros. Otra posibilidad futura del trabajo es la validación por comparación de espectros, previa determinación experimental, o por simulación Monte Carlo, del espectro real.