HTML

ePub

PDF

1. INTRODUCCIÓN

Cada vez más, dos ciencias que a priori podrían parecer distantes como la mecánica y la medicina están cada vez más cerca. Investigadores en ambos ámbitos llevan varios años tratando de simular numéricamente el comportamiento mecánico de órganos humanos para planificar con exactitud biopsias, para mejorar el diagnóstico de lesiones, para la preparación de intervenciones u otros procedimientos clínicos, para simuladores quirúrgicos utilizados en el aprendizaje, etc. [1, 2, 3, 4]. El desarrollo de algoritmos avanzados es necesario para obtener una simulación realista y con el nivel de precisión adecuado para garantizar la seguridad del paciente, el correcto aprendizaje de los futuros cirujanos o la planificación intervenciones quirúrgicas. Ciertos autores han conseguido ya simulaciones numéricas de órganos humanos con suficiente precisión [5, 6], basadas en modelos de comportamiento de material hiperelásticos lineales como el propuesto por Kirchhoff-Saint Venant o no lineales como neoHokean. Otros autores, sin embargo, consideran que es posible hallar modelos de comportamiento mecánico del material que se ajusten con mayor precisión al comportamiento de los órganos humanos tales como el modelo propuesto por Mooney-Rivlin de orden 3 [7].

Sea cual sea el modelo de material utilizado, es fundamental conocer cuáles son las constantes elásticas del modelo que mejor ajustan la simulación al comportamiento real del órgano de un paciente en concreto. La caracterización con precisión de tejidos vivos es uno de los grandes retos para la simulación numérica si requerimos cierto nivel de fiabilidad ya que difícilmente es posible interactuar directamente con el órgano. Para hallar dichas constantes normalmente los investigadores se sirven de un proceso de optimización. Ciertos autores [7] proponen el uso de rutinas de optimización basadas en algoritmos genéticos para hallar el valor de dichas constantes a partir imágenes médicas del hígado de un paciente en dos estados de deformación diferente. Este enfoque pretende evitar la medida invasiva de la respuesta mecánica del órgano y permite, además, que la determinación de estas constantes elásticas del modelo constitutivo sea específica de paciente y pueda realizarse in vivo.

Basándonos en [7], en este trabajo se propone obtener las constantes elásticas del modelo constitutivo de un hígado de un paciente anónimo a partir de dos imágenes médicas del hígado recogidas mediante Tomografías Axiales Computerizadas (TAC). Una de ellas correspondiente al estado de inhalación total y otra al estado de exhalación total. Estas imágenes fueron proporcionadas por el Servido de Radiodiagnóstico del Hospital Universitari i Politècnic La Fe de València y fueron segmentadas (aislando el hígado del resto de órganos del abdomen) usando el software comercial Simpleware®, el cual también proporcionó la malla de elementos usada para la simulación numérica realizada mediante el Método de los Elementos Finitos (MEF). El modelo constitutivo que se asume para modelar el comportamiento mecánico del hígado en este trabajo es un modelo de Mooney-Rivlin de tercer orden. En el procedimiento propuesto, se define una función objetivo que es la distancia euclídea entre ciertos puntos de control del hígado deformado real y los mismos puntos en la simulación numérica para ciertos valores de las constantes elásticas del modelo constitutivo. El conjunto de constantes que definen el modelo constitutivo del material son aquellas que resulten del proceso de optimización que minimice la función objetivo. En trabajos previos [7] fueron usados algoritmos genéticos para el proceso de optimización. Por el contrario, se plantea en este trabajo el uso de algoritmos de optimización basados en el gradiente para minimizar la función objetivo. A diferencia de los algoritmos genéticos, cuyo coste computacional es notoriamente mayor, los algoritmos basados en el gradiente tienen el inconveniente de que pueden converger a mínimos locales. Dicha problemática no ha tenido lugar en las simulaciones realizadas hasta el momento. No obstante, técnicas para paliar esta indeseable característica del método existen en la literatura [8] y podrían ser aplicables en caso de ser necesario.

Las técnicas propuestas se han verificado a través de un experimento sintético cuyo punto de partida fue la imagen del hígado de un paciente anónimo en exhalación total, la cual se asumió como el estado inicial del hígado sin deformar, a partir del cual se construyó la geometría del modelo. La imagen del hígado en inhalación total, estado deformado, se usó para obtener las condiciones de contorno de Dirichlet. Dichas condiciones de contorno se han obtenido mediante un método de registro que proporcionó el desplazamiento que habría que aplicar a los nodos de la superficie externa del hígado sin deformar para alcanzar el estado deformado. Para validar el método se escoge un valor de referencia para las constantes del modelo constitutivo y se realiza una simulación MEF, obteniendo así la posición de todos los puntos en el estado deformado. Posteriormente se lanza una búsqueda ciega de dichas constantes con el método de optimización propuesto.

En el presente trabajo, después de la introducción, se establece el modelo constitutivo y la formulación del problema elástico en la sección 2. La sección 3 trata la formulación del problema de optimización que se validará en los ejemplos numéricos de la sección 4. Finalmente, tras una discusión del método propuesto, en la sección 5, en la sección 6 se destacarán los aspectos más importantes del trabajo.

2. PROBLEMA ELÁSTICO Y MODELO CONSTITUTIVO

Como se ha mencionado anteriormente, el modelado del comportamiento mecánico de un órgano humano de tejido blando, como puede ser el hígado en este caso, pasa por resolver el problema de elasticidad. En este caso, debido a los grandes desplazamientos que se sufren durante el proceso respiratorio, es necesario considerar el problema elástico considerando grandes desplazamientos. Por otro lado, el modelo constitutivo lineal de Kirchhoff-Saint Venant no es capaz de representar adecuadamente el comportamiento de un hígado vivo [7], por ello se requiere de modelos constitutivos no lineales y más complejos. En este trabajo consideramos en modelo propuesto por Mooney-Rivlin que se detalla en los siguientes apartados.

El problema elástico en grandes desplazamientos con un modelo constitutivo no lineal se puede enunciar del siguiente modo: Encontrar el campo de desplazamientos u que satisfaga , entendiendo como el trabajo virtual de las fuerzas internas y el trabajo virtual de las fuerzas externas. Para resolver este problema y aplicar el MEF se considera el uso de la formulación lagrangiana total.

Hay que tener en cuenta que en el problema planteado no se considera la acción de las fuerzas gravitatorias y que todos los contornos externos están sujetos a condiciones de contorno de Dirichlet. Por ello, el trabajo de las fuerzas virtuales externas es nulo en cualquier caso. Debido a la consideración de grandes desplazamientos y la no linealidad del material, el problema elástico es no lineal y para resolverlo se utiliza el método de Newton-Raphson. Para ello es necesario evaluar la matriz de rigidez tangente cuyo proceso de obtención se detalla a continuación para el presente caso.

2.1. Deformación de Green

Para medir la deformación de Green observemos la figura 1 donde en la parte izquierda se representa la configuración inicial del sólido sin deformar, mientras que a la derecha se representa el sólido deformado. Como se observa, cada punto X en la configuración inicial se encuentra en la posición x en el sólido deformado, siendo u el desplazamiento sufrido por dicho punto.

En primer lugar, se definen las coordenadas finales de cada punto como una función de las coordenadas iniciales . El desplazamiento que experimenta cada punto se puede escribir como:

Con esto, el gradiente de deformación se puede escribir del siguiente modo:

siendo F la matriz jacobiana de la transformación entre la configuración inicial y final, . El tensor de deformaciones de Green se define como:

Una vez establecido el tensor de deformaciones para gran- des desplazamientos referido a la configuración inicial se pasa en el siguiente apartado a detallar las diferentes medidas de tensiones.

2.2. Medidas de tensión

En el problema elástico lineal, la medida de tensiones se realiza a través del tensor de tensiones de Cauchy, , también conocidas como tensiones verdaderas ya que son las tensiones de cada punto en la configuración deformada. No obstante, en el caso de grandes desplazamientos en el que nos encontramos nos interesa referir el tensor de tensiones a la configuración inicial y no a la configuración deformada. Para este paso, en primer lugar definimos el primer tensor de tensiones de Piola-Kirchoff, P. Este tensor de tensiones es no simétrico ya que relaciona las fuerzas en la configuración deformada con las áreas en la configuración inicial. La relación entre y P es la siguiente:

donde = det F. Finalmente se define también el segundo tensor de tensiones de Piola-Kirchoff S como:

El segundo tensor de tensiones de Piola-Kirchoff es simétrico y relaciona las fuerzas y las áreas sobre la configuración inicial no deformada. La relación entre el tensor de tensiones de Cauchy y el segundo tensor de tensiones de Piola-Kirchoff es la siguiente:

2.3. Formulación lagrangiana total

Como se ha mencionado anteriormente, el problema elástico se resuelve utilizando el método de los trabajos virtuales en el que se plantea el problema de forma débil integrando en la configuración inicial. Ello da lugar a la formulación lagrangiana total. Debido a que en este caso no se consideran fuerzas externas aplicadas, solo se considera el trabajo virtual de las fuerzas internas, que expresado en la configuración inicial, toma la siguiente forma:

para cualquier desplazamiento virtual . Utilizando la relación entre P y S(5), se obtiene que:

Además, observando la definición del tensor de Green asociado al desplazamiento virtual:

se puede escribir el trabajo virtual de las fuerzas internas en función del segundo tensor de tensiones de Piola-Kirchoff y del tensor de Green de los desplazamientos virtuales en la configuración inicial como:

Con esto, la matriz de rigidez tangente se obtiene a partir de (10) considerando la variación del trabajo virtual de las fuerzas internas con respecto a los desplazamientos virtuales y linealizando la ecuación constitutiva, :

2.4. Modelo hiperelástico de Mooney-Rivlin

La relación entre el estado de deformaciones y el estado tensional al que está sometido un material se denomina modelo constitutivo del material. Como se ha justificado en puntos anteriores, para el cometido de este trabajo, el modelo constitutivo de Mooney-Rivlin de tercer orden proporciona resultados de suficiente precisión para el caso de estudio considerado. El modelo de Mooney-Rivlin define la densidad de energía elástica de deformación como:

donde y son las constantes del modelo constitutivo, siendo la correspondiente al módulo de incompresibilidad. son los invariantes reducidos del tensor de deformaciones de Green-Lagrange C, que se definen como:

siendo los invariantes del tensor de deformaciones de Green-Lagrange C, definidos como:

Con ello, el segundo tensor de tensiones de Piola-Kirchhoff se puede escribir como:

Las derivadas de los invariantes reducidos con respecto al tensor de deformaciones de Green se pueden evaluar como:

siendo:

La matriz de la ecuación constitutiva para el modelo de Mooney-Rivlin se define como:

siendo:

con:

3. ESTIMACIÓN DE PROPIEDADES ELÁSTICAS. OPTIMIZACIÓN

A partir del modelo en grandes deformaciones con el material propuesto en la sección 2 del presente trabajo se puede evaluar el comportamiento mecánico de un hígado con precisión [?, 6, 7]. No obstante, las constantes elásticas del modelo constitutivo no son ni conocidas ni accesibles de manera directa. Por ello se establece, en este trabajo, un procedimiento para evaluar las constantes elásticas del modelo de material de Mooney-Rivlin usando un proceso de optimización basado en el gradiente a partir de la información proporcionada por imágenes médicas. Dichas imágenes son obtenidas en dos posiciones extremas del hígado en un proceso de deformación natural del mismo como es el de respiración. Se eligen, por tanto, los estados de inhalación y exhalación total. El proceso consta de los siguientes pasos: i) segmentación y mallado del hígado en el estado de exhalación, ii) segmentación del hígado el en estado de inhalación (deformado), iii) definición de las condiciones de contorno de Dirichlet de los nodos del contorno externo del hígado mediante el registro de los contornos del hígado usando el algoritmo Coherent Point Drift [9], iv) identificación de puntos de control en ambas geometrías.

Una vez realizado este análisis previo, se establece la función objetivo (21) como la suma de la distancia euclídea entre las posiciones de los N puntos de control reales (o de referencia) en el hígado deformado con las posiciones de los puntos de control predichas mediante el modelo numérico para ciertos valores de las propiedades elásticas del material, :

Notar que para los casos presentados en este trabajo, el valor del módulo de incompresibilidad se toma a partir de la literatura [7] ya que se puede evaluar con precisión en ensayos in vitro. Para hallar las constantes elásticas restantes se resuelve el siguiente problema de optimización:

Para resolver el problema de optimización planteado en (22) se pueden considerar distintos tipos de aproximaciones como algoritmos genéticos o algoritmos basados en el gradiente. Los algoritmos genéticos [7] son capaces de evitar en la mayoría de las situaciones caer en mínimos locales, no obstante, su coste computacional es elevado ya que requieren la realización de gran cantidad de simulaciones. Es decir, resolver gran cantidad de veces el problema elástico en grandes deformaciones planteado en la sección 2 para diversos valores de con el fin de conocer el valor de la función objetivo. Por otro lado, los métodos basados en el gradiente son mucho más eficientes a nivel computacional ya que permiten obtener un óptimo con pocas iteraciones, es decir, con pocas evaluaciones de la función objetivo. Por el contrario, pueden converger a mínimos locales ya que no son capaces de conseguir mínimos globales de manera garantizada a no ser que se usen técnicas avanzadas. A pesar de ello, en los ejemplos numéricos utilizados se ha llegado siempre al mínimo global sin necesidad de utilizar técnicas mejoradas basadas en el método del gradiente.

Para aplicar el método del gradiente es necesario plantear un análisis de sensibilidades de la función objetivo con respecto a las variables a optimizar:

Utilizando la relación x = X + u, la ecuación (23) se puede reescribir como:

Para el presente problema elástico con el modelo de material de Mooney-Rivlin, el análisis de sensibilidades conduce a que:

donde se evalúa en la posición de equilibrio según la ecuación (11). Como se puede observar, el análisis de sensibilidades resulta en un problema lineal ya que se realiza en torno a la posición de equilibrio. Además, teniendo en cuenta (15) y (18) se observa que la matriz tangente depende linealmente de las propiedades elásticas del material. Por ello, se puede evaluar del mismo modo que pero imponiendo y . De modo similar se pueden hallar el resto de sensibilidades.

4. EJEMPLOS NUMÉRICOS

Para validar la técnica propuesta se utilizan dos imágenes de un paciente anónimo, obtenidas a través de un TAC, de su hígado en dos estados de deformación distintos: inhalación total figura 2 y exhalación total figura 3. A partir de la imagen de exhalación total, que se toma como geometría de referencia, se obtiene la geometría en el estado inicial del hígado y su malla mediante el software Simpleware®, tras un proceso de segmentación. La imagen en el estado de inhalación total se utiliza, en primer lugar, para hallar las condiciones de contorno. Dichas condiciones de contorno son los desplazamientos aplicados a la superficie externa del hígado para pasar de la geometría en estado de exhalación total a su geometría en estado de inhalación total. Para conocer dichos desplazamientos se utiliza el algoritmo de registro Coherent Point Drift [9] que permite conocer, con suficiente nivel de precisión, el punto en la geometría final que corresponde a cada punto en la geometría inicial, es decir, el campo de desplazamientos de la superficie. Para la generación del modelo mediante el MEF, se genera una malla tridimensional de elementos tetraédricos lineales. Dicha malla es proporcionada por el software Simpleware® y consta de 49.080 nodos con 216.601 elementos.

Como primer test, se ha resuelto el problema asignan- do unos valores de las constantes elásticas de material obtenidas en bibliografía [7], = 0,018788, = 0,016730 y = 0,038696, fijando el valor = 0,2058. De este modo se obtiene una solución de referencia al problema elástico. Para comprobar la validez del método propuesto, se ha lanzado una batería de pruebas con diferente número de puntos de control. Notar que la posición final de dichos puntos de control es conocida al haber resuelto el problema con los parámetros elásticos de referencia. La prueba consiste ahora en encontrar dichos parámetros mediante el algoritmo propuesto. Para ello fijamos un intervalo de búsqueda cuyo intervalo es: = [0,001, 0,1], = [ 0,1, 0,001] y = [0,001, 0,1]. Los resultados obtenidos de muestran en el cuadro 1.

Como se observa en el cuadro 1, independientemente del número de puntos de control, la precisión alcanzada es la misma. Cabe mencionar que el modelo de material propuesto para este caso no cumple la condición de "poli-convexidad" [10, 11], por ello, al realizar una búsqueda de los parámetros del material en un proceso de optimización basado en el gradiente se pueden encontrar combinaciones de que den lugar a resultados no físicos. Dichas combinaciones deben descartarse durante el proceso de optimización. Finalmente, el método propuesto se ha aplicado a un caso más real en el que los puntos de control no se distribuyen de manera aleatoria dentro del hígado, sino que son aquellos que corresponden a una superficie interna conocida. En este caso, son aquellos pertenecientes a la superficie de un tumor interno, ver figura 4. Aplicando el mismo procedimiento que en los casos anteriores, los resultados obtenidos son = 0,018788, = 0,016730 y = 0,038696, obteniendo de nuevo los valores de las constantes del modelo constitutivo usadas como referencia.

5. DISCUSIÓN

El trabajo presentado en este artículo muestra que es posible obtener una estimación del valor de las constantes elásticas de un modelo constitutivo como el propuesto por Mooney-Rivlin, a partir de la información extraída de imágenes médicas en dos estados, inhalación total y exhalación total. La caracterización del comportamiento mecánico del hígado, en este caso, permitirá conocer las posiciones y movimientos de, por ejemplo, pequeños tumores o lesiones durante el proceso de respiración, con el fin de optimizar su tratamiento y disminuir el riesgo para los pacientes.

En este trabajo hemos supuesto que podemos asumir que el hígado durante el proceso respiratorio sigue un proceso cuasiestático, es decir, se han despreciado los efectos inerciales. Además hemos supuesto que el hígado en exhalación total no se encuentra deformado y que la de- formación se produce por el proceso respiratorio. Esto es una aproximación al proceso real ya que en el estado de exhalación total el hígado se encuentra sometido a las fuerzas gravitatorias y a las interacciones con todos los otros órganos que lo rodean. No obstante, en primera aproximación, no se han tenido en cuenta estas acciones ya que su modelado numérico con cierta precisión se encuentra hoy en día fuera de nuestro alcance. No obstante, a pesar de las simplificaciones adoptadas, la relevancia de los resultados obtenidos es muy alta. Éstos son lo suficientemente precisos como para poder seguir lesiones internas del hígado durante el proceso respiratorio en un paciente in vivo.

6. CONCLUSIONES

La caracterización de las propiedades elástica de tejido blando, como un hígado, es el objetivo principal de este trabajo. Para ello, se asume que el comportamiento del un hígado sigue el modelo constitutivo de Mooney.Rivlin de orden 3, tal y como se indica en la sección 2.4. Para hallar las constantes elásticas, se propone un método de optimización basado en el gradiente que minimice cierta función objetivo. En este caso es la distancia euclídea entre ciertos puntos de control en la posición deformada para cierta combinación de los parámetros del modelo constitutivo y la posición de esos puntos en un caso real o de referencia, como se ha realizado en los ejemplos numéricos. Los resultados muestran que con no demasiados puntos de control o con aquellos relacionados con ciertas estructuras internas conocidas es suficiente para alcanzar el objetivo perseguido en este trabajo.

REFERENCIAS

[1] K. K. Brock, S. J. Hollister, L. A. Dawson, and J. M. Balter, “Technical note: creating a four-dimensional model of the liver using finite element analysis.” Medical physics, vol. 29, no. 7, pp. 1403–5, jul 2002. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/12148719

[2] B. D. de Senneville, M. Ries, G. Maclair, and C. Moonen, “MR-guided thermotherapy of abdominal organs using a robust PCA-based motion descriptor.” IEEE transactions on medical imaging, vol. 30, no. 11, pp. 1987–95, nov 2011. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/21724501

[3] A. Mena, D. Bel, I. Alfaro, D. González, E. Cueto, and F. Chinesta, “Towards a pancreatic surgery simulator based on model order reduction,” Advanced Modeling and Simulation in Engineering Sciences, vol. 2, no. 1, p. 31, dec 2015. [Online]. Available: http://www.amses-journal.com/content/2/1/31

[4] E. Cueto and F. Chinesta, “Real time simulation for computational surgery: a review,” Advanced Modeling and Simulation in Engineering Sciences, vol. 1, no. 1, p. 11, 2014. [Online]. Available: http://amses-journal.springeropen.com/articles/10.1186/2213-7467-1-11

[5] S. Niroomandi, D. González, I. Alfaro, F. Bordeu, A. Leygue, E. Cueto, and F. Chinesta, “Real-time simulation of biological soft tissues: a PGD approach.” International Journal for Numerical Methods in Biomedical Engineering, vol. 29, no. 5, pp. 586–600, may 2013. [Online]. Available: http://www.ncbi.nlm.nih.gov/pubmed/23495247

[6] S. Niroomandi, I. Alfaro, D. González, E. Cueto, and F. Chinesta, “Model order reduction in hyperelasticity: a proper generalized decomposition approach,” International Journal for Numerical Methods in Engineering, pp. n/a–n/a, jul 2013. [Online]. Available: http://doi.wiley.com/10.1002/nme.4531

[7] F. M. Martínez, “Determining the Biomechanical Behaviour of the Liver using Medical Image Analysis and Evolutionary Computation,” Ph.D. dissertation, Universitat Politècnica de Va- lència, 2014.

[8] J. Snyman, Practical Mathematical Optimization, ser. Applied Optimization. New York: Springer-Verlag, 2005, vol. 97. [Online]. Available: http://link.springer.com/10.1007/b105200

[9] A. Myronenko and X. Song, “Point Set Registration: Coherent Point Drift,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 32, no. 12, pp. 2262–2275, 2010. [Online]. Available: http://ieeexplore.ieee.org/xpl/login.jsp?tp={\&}arnumber=5432191{\&}url=http{\%}3A{\%}2F{\%}2Fieeexplore.ieee.org{\%}2Fxpls{\%}2Fabs{\_}all.jsp{\%}3Farnumber{\%}3D5432191

[10] S. Hartmann, “Parameter estimation of hyperelasticity relations of generalized polynomial-type with constraint conditions,” International Journal of Solids and Structures, vol. 38, no. 44-45, pp. 7999–8018, 2001.

[11] S. Hartmann and P. Neff, “Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for near- incompressibility,” International Journal of Solids and Structures, vol. 40, no. 11, pp. 2767–2791, 2003.

IR