Recepción: 11 Diciembre 2017
Aprobación: 12 Abril 2018
Resumen: Se presenta un modelo simplificado de la piel humana adulta como un medio turbio de tres capas para estudiar el comportamiento de la Absorbancia en este tejido debido a la variación de bilirrubina, hemoglobina y melanina mediante una simulación basada en el método de Monte Carlo e implementando un algoritmo en paralelo en el lenguaje de programación python3 para la ejecución del mismo con las librerías subprocess y multiprocessing. Se probo el algoritmo en paralelo tomando 100 archivos iniciales y ejecutándolos de manera secuencial y en paralelo. La implementación consiguió una mejora en el rendimiento de 3,71 y eficiencia de 0,93 para cuatro procesadores. Se confirmó que la melanina es uno de los Cromóforos más influyentes en la absorción de la piel, perdiéndose información óptica de la región inferior a la epidermis para una fracción de volumen de melanosomas mayor a 0,033. Se observó que para concentraciones de bilirrubina mayor a 0,02 g/l hay cambios en la concavidad de la reflexión difusa entre 450 y 500 nm. Se determino que la densidad de vasos sanguíneos y la concentración de hemoglobina tienen efectos similares en la Absorbancia y ambas afectan predominantemente en el espectro entre 400 y 600 nm.
Palabras clave: Monte Carlo, paralelo, absorbancia, reflexión difusa.
Abstract: A simplified adult skin model that proposed it as a three layer turbid medium is presented for studying the Absorbance of skin as a function of melanin, bilirubin and hemoglobin concentration. The model was tested in a simulation based on the Monte Carlo method with the implementation of and algorithm in python3 for running the simulation in parallel with the libraries subprocess and multiprocessing. With data from 100 initial files, the simulation was run sequentially and in parallel. For four processors, the speedup was 3,71 and the efficiency 0,93. It was confirmed that melanin plays and important rol as a Chromophore, losing information from the dermis and hypodermis for Melanosom volumen fraction above 0.033. It was also found that for bilirubin concentrations above 0.02 g/l, there are changes in the shape of the diffuse reflectance curve in the 450 - 500 nm range. The blood vessel density and the hemoglobin concentration were found to affect the Absorbance of the skin in the same way and both affect it in the 400 - 600 nm range mainly.
Keywords: Monte Carlo, parallel, absorbance, diffuse reflectance.
1. Introducción
El transporte de energía en un medio a través de luz esta descrito por la Ecuación de Transferencia Radiativa. Esta es una ecuación integro-diferencial cuya solución analítica no es trivial. En su lugar existe la alternativa de utilizar métodos numéricos para simular el comportamiento de la luz en un medio turbio. Uno de los métodos de elección para la simulación del transporte de fotones es el método de Monte Carlo [1].
El método de Monte Carlo surgió durante el periodo de posterior a la Segunda Guerra Mundial, finalizando la década de los 40. Luego de la creación de la EINAC, acrónimo de Electronic Numerical Integrator And Computer (Computador e Integrador Numérico Electrónico), Stanislaw Ulam, quien había sido parte de la Facultad de Matemáticas de la Universidad de California, tuvo la visión aprovechar la capacidad de computo de la EINAC para utilizar técnicas estadísticas [2]. John Von Neumann, intrigado por las ideas de Ulam, penso en un acercamiento estadístico al problema de la difusión de neutrones en material de fisión, donde los efectos de absorción, esparcimiento, fisión y escape están regidos por el azar. Este problema describe la esencia del método. Hay muchas definiciones para el método MC, pero la que se elije es que: Un método de Monte Carlo (MC) se puede definir como un método numérico en el que se utilizan deliberadamente números aleatorios en un cálculo que tiene la estructura de un proceso estocástico[3, p 2], un proceso cuya evolución este definida por una serie de eventos aleatorios. El resultado se obtiene normalmente realizando una serie de repeticiones del calculo y al final obteniendo un promedio de los resultados. El promedio, según la Ley de Números Grandes, se acercara más al valor esperado a medida que se hagan mayor numero de repeticiones. Esto también produce simulaciones más largas, por lo que se fija usualmente un error aceptable y solo se realizan las repeticiones necesarias. Para solventar el problema del tiempo también se utiliza en ocasiones programación en paralelo, donde se divide el trabajo de computo en varios núcleos de un mismo procesador.
El método MC (Monte Carlo) se ha popularizado y su uso se ha expandido a diversas áreas como la óptica y la física médica. En particular en el área de física médica uno de los métodos de diagnostico alternativo no invasivo existente para enfermedades de la piel, la espectofotometría de reflexión difusa, mide la respuesta óptica espectral de este tejido con el fin de extraer datos de las curvas de reflexión difusa de la misma. La piel puede ser modelada como un medio turbio de múltiples capas: la epidermis formada principalmente por keratinocitos que contienen melanosomas y estos a su vez contienen melanina, el pigmento y cromóforo predominante de la piel, para proteger sus núcleos celulares; la dermis que contiene fibras de colágeno que le dan fortaleza y elasticidad a la piel, además de contener vasos sanguíneos, donde la hemoglobina de los hematocritos y la bilirrubina, un compuesto amarillento que surge del proceso normal de descomposición de la hemoglobina que en altos niveles puede indicar problemas hepáticos, son dos cromóforos importantes; y la capa de tejido subcutaneo o hipodermis, constituida principalmente por grasa.
En este trabajo presenta un modelo simplificado de la piel, para estudiar el efecto de la hemoglobina, la bilirrubina y la melanina en la Absorbancia de la piel adulta y se evalua la efectividad de una implementación de la programación en paralelo en este tipo de simulaciones para aumentar la eficiencia del mismo.
2. Metodología
2.1. Un modelo simplificado de la piel
El modelo de la piel utilizado presenta como un medio de tres capas planas e infinitas horizontalmente: la epidermis, dermis y tejido subcutáneo o hipodermis, tomando como medio superior aire.
La respuesta base del tejido de epidermis y la dermis, sin melanina ni sangre, puede ser aproximada de la data de Iyad Saidi [4] (se logra la expresión utilizando una esfera de integración y tomando en cuenta la contribución extra de ambos pigmentos en las muestras), como
donde λ es la longitud de onda. La data de Saidi viene de muestras in vitro de piel neonatal en el rango entre 450 - 700 nm.
El principal agente absorbente de la epidermis es la melanina que se asume esta dispuesta de manera uniforme en el volumen de la capa y para englobar el resto de la contribución del tejido se tiene µa,base. La absorción de esta capa
(1)donde fme es la fracción volumen de los melanosomas que contienen melanina en la epidermis, los valores de esta fueron tomados de datos publicados [5] y µa,me (λ) es el coeficiente de absorción espectral de los melanosomas, cuya expresión, de acuerdo con data experimental publicada [6], se puede aproximar como
(2)En el caso de la dermis se toman en cuenta dos cromóforos importantes: la hemoglobina y la bilirrubina. En la sangre hay dos tipos de hemoglobina: oxigenada y desoxigenada. Cada una tiene respuestas espectrales diferentes. Los coeficientes de absorción de los cromóforos se pueden despejar de la Ley de Beer-Lambert, suponiendo estos como pigmentos disueltos en un medio con una determinada concentración
donde P peso molar [g/mol], C es la concentración en [g/l] y ε el coeficiente de extinción molar [cm−1/(mol/l)]. Cada una de estas expresiones debe ser escalada por un factor f similar al de la ecuación (2) y los coeficientes de extinción molar ε λ fueron tomados de la data experimental (ver sec 2a de [5]). Por parte de la hemoglobina el balance de la hemoglobina oxigenada y desoxigenada esta dado por la oxigenación en sangre S y
(3)Así el coeficiente de absorción de la dermis para una unidad de volumen sera la suma de la contribución de los constituyentes que absorben luz que en este caso serían la hemoglobina y la bilirrubina con datos de la fracción de volumen de los vasos sanguineos según datos publicados [7].
El esparcimiento en la piel es causado principalmente por la keratina y el colageno en terminos de la contribución del esparcimiento según Mie y según Rayleigh [8]. La expresión del factor de anisotropía g fue tomada de propiedades ópticas publicadas para la epidermis y la dermis [9].
Por último la hipodermis o capa subcutanea se considera principalmente esparsora. El coeficiente de absorción y el factor de anisotropía se toman de datos experimentales como 600 nm-1 y 0,8, respectivamente y ambos constantes en el espectro [10].
2.2. Monte Carlo en multiples capas
El programa Monte Carlo en multiples capas o MCML por sus siglas en ingles (Monte Carlo Multi- layered), es un programa escrito en el lenguaje ANSI C para el modelado de transporte de fotones en un medio de múltiples capas, donde las propiedades ópticas (µs, µa, g y n) y el grosor de cada capa (d) son conocidos [1]. El código describe un haz de fotones de una única longitud de onda e infinitesimalmente delgado que incide perpendicularmente en el medio. El simulador lanza un paquete de fotones con intensidad uno, de manera perpendicular a la superficie del medio. Calcula el camino libre medio antes de que ocurra una interacción, determina el resultado de las interacciones de absorción, esparcimiento y en caso de un cambio de capa refracción de manera aleatoria, hasta que el fotón sale del medio o pierde suficiente intensidad para ser desechado. Luego se sigue con el siguiente fotón hasta que se lancen el número de fotones previamente establecido. Para más detalles del programa ver [1].
Los parámetros simulados son la absorbancia, transmitancia, reflectancia difusa y especular. El MCML guarda el escape de fotones en la parte superior como Rd (r, α) e inferior como T (r, α), donde r es el radio y α el ángulo con respecto a la normal del medio. El programa toma un archivo de entrada en el que se pueden indicar múltiples simulaciones independientes en un mismo archivo de entrada, generando incluso archivos de salida separados para cada simulación, por lo que puede generar data para un estudio espectral con un solo archivo.
Indiferentemente de si se indican datos para múltiples simulaciones en un archivo o estos se separan en archivos distintos, el simulador ejecuta un solo set a la vez, comprendiendo como set los datos necesarios para realizar una simulación (µs, µa, g, n y d para cada capa). Se ejecuta una simulación a la vez, se guardan los resultados en un archivo de salida y se repite para cada archivo (ver Figura 1), con la variación de que el usuario debe ser quien ejecute el simulador para cada archivo en caso de suministrar los sets de datos en archivos de entrada separados. Nótese además que para una sola instancia de parámetros del modelo, habrán tantos sets como longitudes de onda λ.
Figura 1
Ilustración de la manera en la que se ejecuta el MCML de manera secuencial.
2.3. El MCML en paralelo
Se decidió que el punto en el que se realiza el calculo en paralelo sera en la ejecución de cada set de datos en lugar de en cada fotón. Para la implementación en el MCML se modifico la salida de datos en consola que no se consideraron necesarias, pero la simulación en si no se modifico. En su lugar se escribió un programa en el lenguaje python3, de nombre pyMCML, para facilitar y automatizar la generación de archivos iniciales, ejecución del simulador y realización de gráficas. En el pyMCML se utilizan dos librerías para ejecutar el programa en paralelo subprocess y multiprocessing.
La librería subprocess brinda una serie de funciones que permiten ordenar la ejecución de un programa (como por ejemplo el MCML). Por otro lado, multiprocessing permite dividir la ejecución de distintas funciones entre varios núcleos de la maquina. Combinando estas dos funcionalidades se puede ejecutar cualquier programa en paralelo.
La implementación amerita la creación de una función que llamaremos sim(A), donde A son los datos necesarios para la ejecución del programa que se desea ejecutar en paralelo. En este caso A es la ruta de un archivo de entrada y sim, utilizara subprocess para que se ejecute el MCML con este archivo. Luego el número completo de archivos iniciales N = km (k es el número de longitudes de onda y m el número de variaciones de parámetros en el modelo) se divide equitativamente entre el número de procesos independientes que se de desean ejecutar n, quedando n grupos de archivos de entrada G con N n elementos. En caso de que N no sea divisible entre n, número restante de archivos se distribuye entre los grupos y el número de elementos no es igual en cada grupo. Luego, con multiprocessing se definen n procesos para que cada uno ejecute a sim sobre cada A de uno de los grupos G. En la Figura 2 se observa una ilustración del proceso. De esta forma solo se deben preparar los archivos iniciales de manera conveniente para ser utilizados en este algoritmo, colocando un set en cada archivo inicial.
Figura 2
Ilustración elemental de la manera en la que se ejecuta el MCML en paralelo al utilizar pyMCML.
Se ejecutaran 100 archivos iniciales con sets diferentes de manera secuencial y en paralelo y se evaluara la mejora en el rendimiento
donde n es el número de procesadores utilizados y T el tiempo para ese número de procesadores; y la eficiencia
2.4. Simulaciones
Las simulaciones fueron realizadas con una computadora con procesador Intel® Core™i5-3570K de 4 núcleos de 3,40 GHz de 64 bits y 4 Gb de RAM con Ubuntu 16.06. Todas las simulaciones fueron realizadas con 100000 fotones, en el rango de 400-700nm, con una resolución de 4nm para un total de 76 longitudes de onda y con valores de control a menos que se indique lo contrario de fsngr = 0,002, fme= 0,033 (Fototipo III[5]), Cbi = 0,01 g/l y Che = 150 g/l. En la primera simulación se evaluo el efecto de la melanina variando la fracción de melanosomas para todos los Fototipos según Narea [5] fme = {0,015; 0,031; 0,033; 0,049; 0,076} . Para la segunda simulación se probo el efecto de la bilirrubina, variando Cbi = {0,003; 0,01; 0,02; 0,03; 0,04} g/l. Para la tercera y cuarta se evalúan los efectos de la hemoglobina, primero variando el volumen de los vasos sanguíneos fsngr = {0,0005; 0,001; 0,002; 0,003; 0,004} y luego de la concentración en sangre Che = {75, 100, 150, 200, 300} g/l. Con cinco valores por parámetro en cada simulación se generaron 380 archivos iniciales para cada una.
3. Resultados
Al ejecutar los archivos de manera secuencial y en paralelo se observó que los tiempos de computo no son constantes a lo largo del espectro, tomando más tiempo en las longitudes donde la absorción es más débil y por lo tanto con reflexión más fuerte. En estas longitudes de onda los fotones tienen más interacciones débiles con el medio, lo que les permite eventualmente salir de este luego de esparcimientos múltiples, por lo que el tiempo de computo se relaciona con la cantidad de interacciones calculadas antes de pasar al fotón siguiente. En cuanto a la eficiencia, con 107 fotones, en el método secuencial el tiempo de ejecución promedio por fotón fue de 197 µs/fotón con una eficiencia de y para el cálculo en paralelo el promedio fue de 53 µs/fotón utilizando los cuatro núcleos de la computadora. Con el algoritmo en paralelo hubo una mejoría de S(4) = 3,71 en el rendimiento y se tuvo una eficiencia de E(4) = 0,93.
Se muestra el coeficiente de absorción de la epidermis de la primera simulación en la Figura 3. La forma exponencial de la curva es consecuencia del modelo utilizado y de la ecuación (2). La absorción varia proporcionalmente junto con la concentración de melanosomas como consecuencia de ecuación (1).
En la Figura 4 se presentan los coeficientes de absorción de la dermis para la segunda, tercera y cuarta simulación. El aumento en la absorción en la región de 400 a 600 nm observado en las Figuras 4a y 4b es similar y equivalente según las ecuación (3). En la Figura 4c hay un cambio de concavidad entre los 450 y 500 nm a partir de Cbi= 0,02 g/l es cual se considera la máxima concentración saludable, además del aumento proporcional de la absorción. Antes de los 450 nm el cambio no es tan pronunciado como en la variación de hemoglobina. Se denota también que la escala de la absorción por la melanina es mayor que la de otros cromóforos.
Figura 4
Coeficiente de absorción de la Dermis para variaciones de (a) fsngr , (b) Che y (c) Cbi.
En la Figuras 5, 6, 7 y 8 se presentan las gráficas de reflectancia difusa obtenidas de las simulación. En la Figura 5 se ve que la fracción de melanosomas y produce una reflexión difusa más lineal, perdiendo los valles presentes entre 500 y 600 nm producto de la absorción de la dermis. Para los fm = 0,049 y 0,076 no se observa el primer valle entre 400 y 450 nm que se observa en las demás.
En la Figuras 5, 6, 7 y 8 se presentan las gráficas de reflectancia difusa obtenidas de las simulación. En la Figura 5 se ve que la fracción de melanosomas y produce una reflexión difusa más lineal, perdiendo los valles presentes entre 500 y 600 nm producto de la absorción de la dermis. Para los fm = 0,049 y 0,076 no se observa el primer valle entre 400 y 450 nm que se observa en las demás. La intensidad de la luz reflejada es aproximadamente la mitad para los valores mayores a La melanina absorbe la mayor parte de la luz incidente en todo el espectro, impidiendo que esta llegue la dermis.
Figura 5
Refletancia difusa de piel simulada para varios fme.
En la Figura 6 en la región de 400 - 500 nm la curva de reflectancia cambia una tendencia lineal para un concentración normal de bilirrubina, por una concavidad que aumenta con la concentración. La variación se limita a la misma zona donde se produce el cambio de concavidad en la Figura 4c.
Figura 6
Refletancia difusa de piel simulada para varios Cbi.
En las Figura 7 y 8 se observa que la variación de fsngr y Che respectivamente. Ambas tienen comportamientos parecidos entre si, al igual que los cambios en µa,dermis. Se esperaría mayor absorción y menor reflexión al aumentar la Che, ya que al aumentar fsngr , disminuye la contribución del tejido base. Esto no se aprecia en los resultados obtenidos. El cambio más prominente esta entre los 500 y 600 nm. Por debajo de 500 nm, a pesar del incremento en el coeficiente de absorción, no se ve un cambio pronunciado. Esta también es la zona de mayor absorción de la epidermis y por lo que se pierde información de la dermis aún para un fme = 0,033 que es el equivalente a un Fototipo III. A pesar de esto se logra apreciar el valle entre 400 y 450 nm correspondiente al pico de absorción de la dermis.
Figura 7
Refletancia difusa de piel simulada para varios fsngr.
Figura 8
Refletancia difusa de piel simulada para varios Che.
4. Conclusiones
Con una mejoría en el rendimiento de S(4) = 3,71 y E(4) = 0,93 con cuatro núcleos de procesamiento, se desea seguir evaluando el alcance esta implementación del simulador en paralelo. Entre los desafios que se presentan, se encuentra la porción del programa que no se puede paralelizar. Según la Ley de Amdahl para calcular SAm(n, p) de acuerdo al número de nucleos y la porción paralelizable, con un programa donde se pueda paralelizar una fracción 0,98 (98 %) del código, se obtendría SAm(4, 0,98) = 3,77 y EAm(4, 0,98) = 0,94. Estando en el mismo orden que los resultados, tomando 98 % como referencia y duplicando el número de núcleos se obtendría EAm (8, 0,98) = 0,88, EAm(16, 0,98) = 0,77 y EAm (n, 0,98) < 0,5 para n > 51. Esta tendencia asintótica no incluye las limitaciones del hardware en por ejemplo: El acceso de memoria debido a la arquitectura del procesador. Por ello, al aumentar el número de nucleos, estos valores presentarian el tope superior aproximado de la eficiencia que se espera obtener. Aunque aún con una escalabilidad menor a la predicha por Amdahl, el método de implementación puede ser aceptable para otros tipos de simulación que tengan un esquema de entrada y salida de datos similar al del MCML por la facilidad de su implementación. Se planea buscar un sistema con más procesadores para evaluar la escalabilidad del algoritmo y las limitaciones que pueda tener python.
En cuanto al modelo, se puede acotar que la melanina es un compuesto altamente absorbente que absorbe la mayor parte de la luz, como se referencia en la literatura. Para concentraciones fme > 0,033 se reduce la cantidad de luz difusa reflejada a tal grado que empieza a perder información óptica de la dermis. Este efecto debe tomarse en cuenta al realizar cualquier medición en sujetos de piel obscura. Para la bilirrubina hay cambios significativos en la tendencia de la reflexión difusa en la zona de 450 500 nm para Cbi > 0,02 g/l, que se considera el tope para el funcionamiento normal del cuerpo, donde la región pasa de tener una tendencia lineal a ser concava.
Al comparar los efecto del aumento de la fracción de volumen de vasos sanguineos y de la concentración de hemoglobina en sangre, se determino que bajo este modelo el tejido base, limpio de pigmentos, no afecta de manera importante la absorción de la piel, por lo que se considera que se puede despreciar. El aumento de hemoglobina reduce la reflexión en la zona entre 400 y 600 nm.
5. Referencias
[1] Lihong Wang, Steven L. Jacques, and Liqiong Zheng. MCML–Monte Carlo modeling of light transport in multi-layered tissues. Computer Methods and Programs in Biomedicine, 47(2):131–146, jul 1995.
[2] N. Metropolis. The beginning of the Monte Carlo method. Los Alamos Science, (5):125–130, 1987.
[3] Malvin H. Kalos and Paula A. Whitlock. Monte Carlo Methods Volume 1: Basics. John Wiley & Sons, Inc., 1986.
[4] Iyad Salam Saidi. Transcutaneous optical measurement of hyperbilirubinemia in neonates. PhD thesis, Rice University, USA, 1992.
[5] F. Narea, S. Vivas, and A. Muñoz. Retrieving the absorption coefficient of epidermis in human skin. Optica Pura y Aplicada, 48(3):207–214, sep 2015.
[6] Steven L. Jacques and Daniel J. McAuliffe. The melanosome: Threshold temperature for explosive vaporization and internal absorption coefficient during pulsed laser irradiation. Photochemistry and Photobiology, 53(6):769–775, jun 1991.
[7] J.A. Delgado Atencio, S.L. Jacques, and S. Vázquez y Montiel. Monte Carlo Modeling of Light Propagation in Neonatal Skin. In Charles J. Mode, editor, Applications of Monte Carlo Methods in Biology, Medicine and Other Fields of Science, chapter 17, pages 297–314. InTech, feb 2011.
[8] Steven L. Jacques. Optical properties of biological tissues: a review. Physics in Medicine and Biology, 58(11):R37–R61, may 2013.
[9] M.J.C. Van Gemert, S.L. Jacques, H.J.C.M. Sterenborg, and W.M. Star. Skin optics. IEEE Transactions on Biomedical Engineering, 36(12):1146–1154, 1989.
[10] Tsumura Norimichi, Kawabuchi Miki, Haneishi Hideaki, and Miyake Yoichi. Mapping pigmentation in human skin by multivisible-spectral imaging by inverse optical scattering technique. Journal of Imaging Science and Technology, 45(5):444–450, sep 2001.
Enlace alternativo
http://servicio.bc.uc.edu.ve/ingenieria/revista/v25n2/art04.pdf (pdf)
