Resolución numérica de un modelo de frontera libre para el crecimiento tumoral

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Resoluci´on num´erica de un modelo de frontera libre para el crecimiento tumoral Carmen Calzada1, Gema Camacho2, Enrique Fern´ andez-Cara2, Mercedes Mar´ ın1 1Dpto. de Inform´atica y An´alisis Num´erico, Universidad C´ordoba, Campus de Rabanales, Ed. C2-3, E-14071 C´ordoba. E-mail: [email protected], [email protected]. 2Dpto. E.D.A.N., Universidad de Sevilla, Aptdo. 1160, E-41080 Sevilla. E-mails: [email protected], [email protected]. Palabras clave: Crecimiento tumoral, frontera libre, dominios ficticios, level set, elementos finitos Resumen En este trabajo presentamos un esquema de resoluci´on num´erica de un modelo de frontera libre que simula el crecimiento de un tumor entre la fase avascular y vascular. Para la resoluci´on se ha utilizado un esquema que permite desacoplar las diferentes inc´ognitas en cada etapa de tiempo. Los problemas de Dirichlet que aparecen se resuelven mediante un m´etodo de dominios ficticios para no tener que remallar el dominio en cada etapa. Los nuevos dominios se calculan utilizando un m´etodo de tipo “level set”. La ecuaci´on asociada se resuelve utilizando t´ecnicas de estabilizaci´on de tipo GLS. Para la discretizaci´on espacial se emplean elementos finitos P2, a fin de poder calcular la curvatura de la frontera de forma m´as precisa. 1. Introducci´on Uno de los retos m´as grandes de la ciencia moderna es entender la din´amica del crecimiento de tumores cancerosos. En las ´ultimas d´ecadas, se ha desarrollado una gran variedad de modelos, basados en EDP’s, que intentan simular el crecimiento de este tipo de tumores [1]. Estos modelos est´an basados en leyes de conservaci´on de la masa y en procesos de reacci´on-difusi´on que tienen lugar dentro del tumor. En principio, seg´un la fase que se quiera estudiar (avascular, angiog´enesis, vascular), el modelo ser´a diferente, ya que los procesos biol´ogicos que tienen lugar y las variables que intervienen en cada una de ellas son de naturaleza distinta. Dada la complejidad de estos modelos, a´un son pocos los trabajos dedicados a su resoluci´on num´erica, sobre todo cuando se trata con problemas 1 C. Calzada, G. Camacho, E. Fern´andez-Cara, M. Mar´ın tridimensionales en geometr´ıas arbitrarias. M´as escasos a´un son los trabajos que utilizan elementos finitos para la discretizaci´on espacial. En general, los modelos que simulan la fase avascular vienen descritos por problemas de frontera libre en los que interesa estudiar c´omo cambia la forma del tumor con el tiempo. En este trabajo presentamos un esquema de resoluci´on num´erica de un modelo de frontera libre propuesto en [4] (ver tambi´en [6] y [8]) que simula el crecimiento de un tumor s´olido (carcinoma) entre la fase avascular y vascular, sin tener en cuenta el proceso de angiog´enesis intermedio. El tumor es tratado como un fluido incompresible, no teni´endose en cuenta la elasticidad del tejido. Las fuerzas de adhesi´on c´elula-c´elula se modelan mediante una tensi´on superficial en el tejido tumoral exterior. El crecimiento de la masa del tumor viene gobernado por un equilibrio entre la mitosis celular y la apoptosis (muerte programada de las c´elulas). La velocidad de la mitosis depende de la concentraci´on de nutrientes dentro del tumor, donde se supone que la concentraci´on de capilares sangu´ıneos es uniforme. Se considera que el tumor no tiene zona necr´otica formada por c´elulas muertas y que no existen inhibidores qu´ımicos en los tejidos externos. Estas condiciones se dan en el caso de tumores de peque˜no tama˜no o cuando la concentraci´on de nutrientes en la sangre o en los tejidos exteriores es alta. No obstante, el modelo que utilizamos puede tambi´en servir para predecir el comportamiento en otras situaciones [6]. 2. Descripci´on del modelo Denotamos ω(t) el dominio ocupado por el tumor en un instante tyγ(t) la frontera que lo separa de los tejidos sanos. Denotamos σ=σ(x, t) la concentraci´on de nutrientes dentro del tumor (por ejemplo, ox´ıgeno o glucosa). ´ Esta viene gobernada por una ecuaci´on de reacci´on-difusi´on y, puesto que la velocidad con la que se difunden los nutrientes es del orden de minutos mientras que la velocidad con la que proliferan las c´elulas es del orden de d´ıas, podemos considerar que se parte de un estado estacionario. Suponemos que los vasos sangu´ıneos existentes aportan nutrientes con una velocidad lineal y que las c´elulas tumorales los absorben de igual forma. Esto da lugar a la ecuaci´on para σen ω(t): Dσ∆σ−γs(σ−σs)−δσσ= 0,(1) con Dσ, γsyδσconstantes positivas y σs> σ la concentraci´on de nutrientes en sangre. El tumor se modela como un fluido incompresible. El correspondiente campo de velocidades u=u(x, t) verifica la ecuaci´on de continuidad en ω(t): ∇ · u=γTσ−δT,(2) donde el t´ermino γTσ(γT>0) se debe a la proliferaci´on (mitosis) de las c´elulas tumorales yδT>0 determina la muerte (apoptosis) de las mismas. El campo de velocidades en ω(t) est´a relacionado con el gradiente de la presi´on por la ley de Darcy u=−wT∇p, (3) siendo wT>0 el coeficiente de movilidad de las c´elulas tumorales (que suponemos constante) y p=p(x, t) la presi´on dentro del tumor. Combinando las ecuaciones (2) y (3), se 2 Resoluci´on num´erica de un modelo de crecimiento tumoral llega a que ∆p=−γT wT σ+δT wT en ω(t).(4) Se supone que la concentraci´on de nutrientes sobre la frontera del tumor coincide con la de los tejidos que lo rodean, que consideraremos uniforme: σ=σγ(t)sobre γ(t).(5) La presi´on sobre la frontera satisface la relaci´on de Young-Laplace: p=aκ sobre γ(t),(6) donde aes la tensi´on superficial (correspondiente a las fuerzas de adherencia c´elula-c´elula) yκes la curvatura local de γ(t). La velocidad normal con la que se mueve la frontera del tumor viene dada por Vn=u·n=−wT∇p·nsobre γ(t),(7) siendo nel vector normal exterior unitario a γ(t). Una vez adimensionalizadas las ecuaciones (ver [8]), llegamos al siguiente problema a resolver: σ−∆σ= 0 en ω(t),(8) σ= 1 sobre γ(t),(9) −∆p=µ(σ−˜σ) en ω(t),(10) p=κsobre γ(t),(11) ∂p ∂n =−Vnsobre γ(t),(12) para t∈(0, T). Por simplicidad, hemos denotado las nuevas variables adimensionalizadas y los nuevos dominios de la misma forma que los antiguos. La formulaci´on utilizada permite describir la evoluci´on del tumor mediante dos par´ametros: µ, que representa la raz´on entre la mitosis y los mecanismos de relajaci´on (movilidad celular y adhesi´on c´elula-c´elula) y ˜σ, que describe el equilibrio entre la apoptosis y la mitosis. Ambos par´ametros incluyen tambi´en los efectos de la vascularizaci´on. As´ı, hay tres grados de vascularizaci´on que pueden identificarse en t´erminos de µy ˜σ: baja vascularizaci´on (µ≥0, ˜σ > 0), en el cual domina la difusi´on como ocurre en experimentos in vitro, vascularizaci´on moderada (µ≥0, ˜σ≤0) y alta vascularizaci´on (µ < 0, ˜σ < 0 ´o ˜σ > 0), que corresponden a reg´ımenes observados en experimentos in vivo (para m´as detalles, ver [6]). 3. Resoluci´on num´erica Para la resoluci´on de (8)–(12) se ha utilizado un esquema que permite desacoplar las diferentes inc´ognitas en cada etapa de tiempo. Se divide el intervalo temporal [0, T] en Msubintervalos de amplitud ky se parte del dominio inicial ω0=ω(0) y de su frontera γ0=∂ω0. Supuesto conocido el dominio ωmen un instante mk para m≥0, se calculan σm,pmyωm+1 (aproximaciones de σ(·, mk), p(·, mk) y ω((m+ 1)k)) como sigue: 3 C. Calzada, G. Camacho, E. Fern´andez-Cara, M. Mar´ın Nutrientes. Se resuelve el problema de Dirichlet σm−∆σm= 0 en ωm, σm= 1 sobre γm.(13) Presi´on. Se resuelve el problema de Dirichlet −∆pm=µ(σm−˜σ) en ωm, pm=κmsobre γm,(14) donde κmes la curvatura de γm. Nuevo dominio. Se calcula la velocidad normal Vm nen el instante mk: Vm n=−∇pm·n.(15) A partir de ella se obtiene el nuevo dominio ωm+1, se hace m=m+ 1 y se repite el proceso. 3.1. El m´etodo de dominios ficticios Los problemas (13) y (14) se resuelven utilizando un m´etodo de dominios ficticios con multiplicadores de Lagrange distribuidos en volumen (ver [7]). Ambos problemas se pueden escribir en la forma siguiente: dado ω⊂Rd(d≥1) de frontera γ, hallar u∈H1(ω) tal que αu −β∆u=fen ω, u=gsobre γ, (16) con f∈H−1(ω), g∈H1/2(γ), α≥0 y β > 0. Necesitamos resolver muchos problemas de este tipo en dominios ωque van cambiando. Interesa un m´etodo que permita resolverlos sin necesidad de modificar el mallado. Para ello, buscamos un nuevo dominio Ω de frontera Γ que incluya a todos los anteriores y formulamos en ´el un problema variacional equivalente a (16). Denotamos ωc= (Ω \¯ω) y suponemos definidas ˜g∈H1(ωc) y ˜ f∈L2(Ω), tales que ˜g|Γ= 0,˜g|γ=gy˜ f|ω=f. Definimos los espacios V0=v∈H1 0(Ω) : v= 0 en ωc, Vg=v∈H1 0(Ω) : v= ˜gen ωc. El problema (16) es equivalente a encontrar ˜u∈Vgtal que ZΩ (α˜uv +β∇˜u· ∇v)dx =ZΩ ˜ fvdx ∀v∈V0,(17) en el sentido de que la soluci´on ude (16) es la restricci´on a ωde la soluci´on ˜ude (17). La relaci´on ˜u= ˜gen ωcse impone mediante un multiplicador de Lagrange λ∈Λ, con Λ = {v∈H1(ωc) : v|Γ= 0}, 4 Resoluci´on num´erica de un modelo de crecimiento tumoral dotado del producto escalar s(µ, µ0) = αZωc µµ0dx +βZωc ∇µ· ∇µ0dx ∀µ, µ0∈Λ.(18) De esta forma, se sustituye el problema (17) por el de buscar el punto de silla del Lagrangiano L:H1 0(Ω) ×Λ7→ Rdefinido como L(v, µ) = 1 2ZΩαv2+β|∇v|2dx −ZΩ ˜ fvdx +s(µ, v|ωc−˜g).(19) Este punto de silla se aproxima utilizando un m´etodo de Uzawa con paso ´optimo que resulta ser bastante efectivo. Para la discretizaci´on espacial se utilizan elementos finitos P2–Lagrange, a fin de poder calcular la curvatura de la frontera de forma m´as precisa. Una ventaja de utilizar dominios ficticios es que disponemos de la presi´on pen todo el dominio Ω y, por tanto, la velocidad normal Vntambi´en se tiene ampliada. Esto es muy conveniente para calcular el dominio ocupado por el tumor en la siguiente etapa temporal. 3.2. El m´etodo “level set” Para calcular el nuevo dominio ω(t) en cada paso de tiempo, se utiliza un m´etodo de tipo “level set”. Estas t´ecnicas fueron introducidas por Osher y Sethian en un contexto general en [10] y permiten seguir la evoluci´on de una frontera m´ovil. La idea b´asica es definir una funci´on escalar ψ=ψ(x, t) (funci´on de nivel) en todo el dominio computacional Ω, de tal manera que su valor en cada punto x∈Ω indica si se est´a dentro o fuera del dominio ω(t): ψ(x, t)   <0 si x∈ω(t), = 0 si x∈γ(t), >0 si x/∈ω(t). (20) La frontera γ(t) viene dada, por tanto, como el conjunto de nivel cero de ψ. Esta funci´on satisface en Ω la ecuaci´on ψt+Vn|∇ψ|= 0,(21) donde Vnes la velocidad normal de la frontera γ(t) ampliada a todo Ω. Como valor inicial ψ(x,0) se toma la funci´on distancia con signo a la frontera γ(0) (negativa dentro y positiva fuera). Esta formulaci´on posee varias ventajas sobre los m´etodos expl´ıcitos de seguimiento de frentes. As´ı, el dominio ocupado por el tumor en cada etapa de tiempo se identifica simplemente por el signo de la funci´on de nivel. Las propiedades geom´etricas de la frontera son f´aciles de determinar. Por ejemplo, en cualquier punto de la misma el vector normal ny la curvatura κvienen dados por n=∇ψ |∇ψ|, κ =∇ · ∇ψ |∇ψ|.(22) Adem´as, la misma formulaci´on es v´alida independientemente de la dimensi´on espacial y los cambios topol´ogicos en la evoluci´on de la frontera se pueden tratar de forma natural. El conjunto de nivel cero de ψno se limita necesariamente a una ´unica curva y puede romperse o mezclarse conforme tavanza. 5 C. Calzada, G. Camacho, E. Fern´andez-Cara, M. Mar´ın Al ser (21) una ecuaci´on hiperb´olica, para su resoluci´on es necesario aplicar un m´etodo de elementos finitos estabilizado. En nuestro caso, hemos utilizado un m´etodo de tipo GLS (aproximaci´on de Galerkin con estabilizaci´on por m´ınimos cuadrados) [2], [3]. Se emplea una discretizaci´on simult´anea espacio-tiempo [13]. La formulaci´on por elementos finitos se escribe para una sucesi´on de bandas espacio-tiempo Qm= Ω×Imcon Im= [tm +, tm+1 −]. La discretizaci´on de una banda Qmse lleva a cabo dividi´endola en elementos QT m=T×Im con T∈ Th, siendo Thuna triangulaci´on de Ω. Definimos el espacio de interpolaci´on Vh m=nvh∈C0(Qm) : vh|T×Im∈Pk(T)×Pk0(Im)∀T∈ Tho.(23) Las funciones de interpolaci´on son continuas dentro de cada franja, pero pueden ser discontinuas en tiempo, al pasar de una franja temporal a otra. Denotamos (·)− my (·)+ mlos valores de la funci´on en tmen la banda de abajo y en la de arriba, respectivamente. La formulaci´on variacional estabilizada para (21) se escribe entonces as´ı: dada (φh)− m, hallar φh∈Vh mtal que, para todo wh∈Vh m, se tiene: ZImZΩ wh(φh t+Vh n|∇φh|)dx dt +ZΩ (wh)+ m((φh)+ m−(φh)− m)dx +X T∈Th ZImZT (wh t+Vh n|∇wh|)τ(φh t+Vh n|∇φh|)dx dt = 0,(24) siendo τ > 0 un par´ametro de estabilizaci´on convenientemente elegido (ver [3]). Para la discretizaci´on en tiempo se han considerado tanto funciones constantes como lineales a trozos que dan lugar a esquemas expl´ıcitos. Para la discretizaci´on en espacio se han considerado funciones lineales y cuadr´aticas. Al haber utilizado esquemas expl´ıcitos de discretizaci´on en tiempo, ha de cumplirse una condici´on CFL, en concreto que ∆t < h m´axΩ|Vn|.(25) Esto significa en la pr´actica que en un paso de tiempo un punto de la frontera del tumor no puede atravesar m´as de dos elementos espaciales. Como apunt´abamos en la secci´on anterior, al haber utilizado el m´etodo de dominios ficticios para la resoluci´on de los problemas de Dirichlet, disponemos de la velocidad normal Vnampliada de forma natural a todo Ω, lo que evita tener que realizar esta ampliaci´on en cada etapa temporal. El c´alculo de esta velocidad ampliada es uno de los inconvenientes que presentan normalmente los m´etodos “level set”. A pesar de que, tradicionalmente, se suelen considerar esquemas de alto orden en tiempo, en la pr´actica, para los problemas que estamos resolviendo, se observa que el m´etodo “level set” es menos sensible a la precisi´on temporal (una vez que se ha elegido el paso de tiempo adecuado para garantizar la convergencia) que a la precisi´on espacial. Por tanto, un esquema en tiempo de primer o segundo orden ofrece a menudo resultados seguros, reduciendo as´ı el coste computacional (ver [8]). Sin embargo, las soluciones num´ericas obtenidas s´ı se muestran sensibles a la precisi´on espacial, sobre todo en un r´egimen de baja vascularizaci´on. Hay que hacer notar que, aunque la condici´on inicial ψ(·,0) sea una funci´on regular, puede perder esta regularidad cuando pasan varios pasos de tiempo, apareciendo problemas num´ericos. Puesto que s´olo interesa encontrar el contorno que define la frontera, es 6 Resoluci´on num´erica de un modelo de crecimiento tumoral Figura 1: Crecimiento en t= 0, t= 0,5 y t= 1. Figura 2: Concentraci´on de nutrientes (izquierda) y presi´on (derecha) en t= 1. posible regularizar la funci´on manteniendo dicho contorno. En la literatura pueden encontrarse varias t´ecnicas de regularizaci´on (ver [12],[11]). Estas t´ecnicas se suelen utilizar en cada paso de tiempo, lo cual puede resultar bastante costoso. Sin embargo, los resultados num´ericos que presentamos se han conseguido sin necesidad de regularizar la funci´on de nivel calculada lo cual es interesante desde el punto de vista computacional. 4. Resultados Num´ericos Se ha considerado como frontera inicial del tumor una circunferencia perturbada, definida por las ecuaciones param´etricas: (x(α), y(α)) = (2,1+0,5 cos(2α))(cos(α),sin(α)), α ∈[0,2π].(26) Para comparar con resultados existentes, se han tomado como par´ametros del modelo µ= 20 y ˜σ= 0,5, que corresponde a un r´egimen de baja vascularizaci´on. Como paso de discretizaci´on espacial se ha tomado h= 0,5 y como paso de tiempo k= 0,01. En la figura 1 aparece la evoluci´on del tumor en diferentes etapas de tiempo. El comportamiento coincide con el que aparece en [6]. En la figura 2 se presentan la concentraci´on de nutrientes y la presi´on en el interior del tumor en el instante t= 1. 5. Conclusiones Se presenta un m´etodo num´erico de resoluci´on de un modelo de crecimiento de un tumor s´olido entre la fase avascular y vascular. El m´etodo resulta ser bastante efectivo, al 7 C. Calzada, G. Camacho, E. Fern´andez-Cara, M. Mar´ın combinar t´ecnicas de dominios ficticios con t´ecnicas de tipo “level set”, ya que, gracias a las primeras, evitamos tener que realizar la ampliaci´on de la velocidad normal a todo el dominio computacional. Tambi´en se evita la regularizaci´on de la funci´on de nivel en cada paso de tiempo. En este modelo simplificado no se han tenido en cuenta los mecanismos asociados a la angiog´enesis. Estos mecanismos y otros adicionales ser´an incorporados al modelo en trabajos posteriores. Se tiene previsto utilizar otros m´etodos de estabilizaci´on y de discretizaci´on temporal para el c´alculo de la funci´on de nivel, an´alogos a los que aparecen en [9]. Tambi´en se incorporar´an t´ecnicas de regularizaci´on cada cierto tiempo, a fin de comparar con los resultados obtenidos en este trabajo. Agradecimientos Este trabajo ha sido financiado por los proyectos EXC/2005/FQM-520 de la Junta de Andaluc´ıa y MTM2006-07932 del Ministerio de Educaci´on y Ciencia. Los autores agradecen a Jordi Blasco, de la UPC, los comentarios y sugerencias realizados. Referencias [1] Araujo R.P., McElwain D.L.S., A history of the study of tumor growth: the contribution of mathematical modeling, Bull. Math. Biol., 66, 1039-1091 (2004). [2] Barth, T.J. , Sethian, J.A. , Implementation of Hamilton-Jacobi and Level-Set Equations on Triangulated Domains, von Karman Institute Lectures Series, Computational Fluid Mechanics (1998). [3] Barth, T.J. , Sethian, J.A. , Numerical Schemes for the Hamilton-Jacobi and Level-Set Equations on Triangulated Domains, J. Comp. Phys., 145(1), 1–40 (1998). [4] Bazaliy, B.V. , Friedman, A. , A Free Boundary Problem for a Elliptic-Parabolic System: Application to a Model of Tumor Growth. Comm. Partial Differential Equations, 28 ,3-4, 517–560, (2003). [5] Brooks A., Hughes T., Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering 32, 199–259 (1982). [6] Cristini V., Lowengrub J., Nie Q. Nonlinear simulation of tumor growth, J. Math. Biol 46, 191–224 (2003). [7] Glowinski, R, Numerical methods for fluids (Part 3), Handbook of Numerical Analysis, Vol. IX Ed. North-Holland, Amsterdam (2003). [8] Hogea, C.S., Murray, B.T., Sethian, J.A., Implementation of the level set method for continuum mechanics based tumor growth models, FDMP 1(2), 109–130 (2005). [9] Malidi, A., Dufour, S., N’dri, D., A study of time integration schemes for the numerical modelling of free surface flows, Int. J. Numer. Meth. Fluids 48(10),1123–1147 (2005). [10] Osher, S., Sethian, J.A., Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulation, J. Comp. Phys. 79, 12–49 (1988). [11] Osher, S., Fedkiw, R., Level Set Methods and Dynamic Implicit Surfaces, Springer–Verlag, New York(2002). [12] Sethian, J.A. Level set methods and fast marching methods. Evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Second edition. Cambridge University Press, Cambridge, (1999). [13] Tezduyar, T.E., Finite Element Methods for Flow Problems with Moving Boundaries and Interfaces, Archives of Computational Methods in Engineering, 8, 83–130 (2001). 8