Control de EDPs orientado a la terapia de un tumor cerebral

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–9) Control de EDPs orientado a la terapia de un tumor cerebral R. Echevarr´ ıa,1, A. Doubova1, E. Fern´ andez-Cara1 I. Gayte-Delgado 1 1Dpto. E.D.A.N., Universidad de Sevilla, Aptdo. 1160, E-41080 Sevilla. E-mails: [email protected], [email protected], [email protected], [email protected]. Palabras clave: control, EDP, tumor cerebral, terapia de tumores cerebrales Resumen Consideramos un sistema de ecuaciones en derivadas parciales que modela los efectos de una terapia sobre un tumor cerebral (glioblastoma). La densidad de c´elulas tumorales verifica una EDP parab´olica semilineal, que est´a acoplada con otra EDP similar para un agente citot´oxico. El control es distribuido, su soporte es un peque˜no subdominio del abierto que representa el cerebro y act´ua a trav´es del segundo miembro de la ecuaci´on para la densidad de anticuerpos. Presentamos resultados de control para este modelo, as´ı como algunas experiencias num´ericas. 1. El modelo considerado Sea Ω ⊂RNun abierto conexo acotado y sea T > 0. Denotamos Q= Ω ×(0, T) y Σ = ∂Ω×(0, T). Consideramos el siguiente problema, que describe la evoluci´on en Ω×(0, T) de un tumor cerebral (glioblastoma):              ct−∇·(D(x)∇c) = f(c)−F(c, β) en Q, βt−µ∆β=h(β)−H(c, β) + v1ωen Q, ∂c ∂n = 0,∂β ∂n = 0 sobre Σ, c(0) = c0, β(0) = β0en Ω. (1) Aqu´ı, se interpreta que Ω es el cerebro; c=c(x, t) y β=β(x, t) son respectivamente las concentraciones de c´elulas tumorales y de anticuerpos (agentes citot´oxicos) generados por el organismo; D=D(x) es el coeficiente de difusi´on de c´elulas tumorales, por simplicidad, 1 R. Echevarr´ıa, A. Doubova, E. Fern´andez-Cara, I. Gayte-Delgado se supone que el coeficiente de difusi´on de anticuerpos es constante e igual a µ;fyhson funciones que determinan respectivamente los ritmos de producci´on de cyβ; por otra parte, la manera en que la interacci´on de c´elulas y anticuerpos afecta a sus respecivas evoluciones viene dada por las funciones FyH. De acuerdo con lo expuesto en [1] y [2] (donde se considera un modelo m´as simple de crecimiento del tumor cerebral y sin la acci´on de terapia), supondremos que D(x) = ½Dwsi x∈Ωw, Dgsi x∈Ωg,(2) donde 0 < Dw< Dg(Ωwy Ωgson las zonas del cerebro respectivamente ocupadas por la materia blanca y la materia gris). Denominaremos Λ : D(Λ) ⊂L2(Ω) 7→ L2(Ω) el operador de difusi´on asociado, con D(Λ) = {w∈H1(Ω) : ∇ · (D(x)∇w)∈L2(Ω) },Λw=∇ · (D(x)∇w)∀w∈D(Λ). Supondremos tambi´en que las funciones f,h,FyHest´an dadas por f(c) = a1c, h(β) = a2β, F(c, β) = b1cβ, H(c, β) = b2cβ, (3) donde las constantes a1,a2,b1yb2son positivas (en particular, aceptaremos que la interacci´on de anticuerpos y c´elulas tumorales es instant´anea y despreciaremos efectos de retardo; tampoco tendremos en cuenta efectos de saturaci´on). En (1), los datos iniciales deben verificar c0, β0∈H1(Ω) ∩L∞(Ω), c0, β0≥0.(4) Por otra parte, v=v(x, t) es un control; en la pr´actica, supondremos que v∈L∞(ω×(0, T)).(5) Cada funci´on vdescribe una terapia que est´a siendo aplicada a lo largo del intervalo temporal (0, T)1. Obviamente, se espera que vdetermine un correcto incremento de anticuerpos que, a su vez, haga decrecer los valores de ca trav´es del t´ermino −F(c, β), en el segundo miembro de la primera ecuaci´on de (1). Para el sistema precedente tenemos el resultado que sigue, cuya demostraci´on usa argumentos bien conocidos: Teorema 1 Supongamos que se cumplen las condiciones (2)–(5). Entonces (1) posee exactamente una soluci´on (c, β), con ½c∈L2(0, T;D(Λ)) ∩L∞(Ω ×(0, T)), ct∈L2(Ω ×(0, T)), β∈L2(0, T;H2(Ω)) ∩L∞(Ω ×(0, T)), βt∈L2(Ω ×(0, T)).(6) Adem´as, c≥0. 1En este modelo, la terapia es aplicada exclusivamente en los puntos de ω. De modo an´alogo, se podr´ıa formular un problema en el que la terapia se aplica sobre parte de la frontera. 2 Control de EDPs orientado a la terapia de un tumor cerebral 2. El problema de control ´optimo Consideremos el funcional J:L∞(ω×(0, T)) 7→ Rdefinido por J(v) = a 2ZΩ |c(T)|2dx +b 2ZZ ω×(0,T ) |v|2dx dt, (7) donde ces, junto con β, soluci´on del sistema (1) y aybson dos par´ametros positivos, con a+b= 1. La funci´on Jpuede interpretarse como un promedio de la cantidad de c´elulas cancer´ıgenas que quedan al final del tratamiento y la intensidad de la medicaci´on empleada. Los pesos aybdeterminan la importancia de cada t´ermino. Pretendemos hallar el control vque minimice este promedio. Para ello planteamos el siguiente problema de control ´optimo: ½Minimizar J(v) Sujeto a v∈Uad,(8) donde Uad ⊂L∞(ω×(0, T)) es un convexo no vac´ıo, cerrado para la norma habitual en L2(ω×(0, T)). Desde el punto de vista m´edico, es imprescindible restringir la medicaci´on. Por tanto, entre otras elecciones, una limitaci´on natural consiste en tomar Uad ={v∈L∞(ω×(0, T)) : 0 ≤v≤M}.(9) Se tiene el resultado siguiente, cuya demostraci´on es est´andar: Teorema 2 Existe al menos una soluci´on del problema de control ´optimo (8). 3. El sistema de optimalidad Sea ˆvla soluci´on de (8) que proporciona el Teorema 2, y sea (ˆc, ˆ β) la soluci´on del correspondiente sistema (1), es decir, el estado asociado2. Aplicando un teorema generalizado de multiplicadores de Lagrange3, obtenemos que (ˆv, ˆc, ˆ β), junto con la soluci´on (d, η) del sistema adjunto asociado                    −dt−∇·(D(x)∇d) = f0(ˆc)d−∂F ∂c (ˆc, ˆ β)d−∂H ∂c (ˆc, ˆ β)ηen Q, −ηt−µ∆η=h0(ˆ β)η−∂H ∂β (ˆc, ˆ β)η−∂F ∂β (ˆc, ˆ β)den Q, ∂d ∂n = 0,∂η ∂n = 0 sobre Σ, d(T) = ˆc(T), η(T) = 0 en Ω, (10) deben verificar ZZ ω×(0,T ) (aη +bˆv)(v−ˆv)dx dt ≥0∀v∈Uad,ˆv∈Uad.(11) 2Lo que sigue tambi´en se aplicar´ıa a cualquier soluci´on local de (8). 3Por ejemplo, un resultado de este tipo puede ser deducido a a partir del llamado formalismo de Dubovitskii-Milyutin, v´ease [3]. 3 R. Echevarr´ıa, A. Doubova, E. Fern´andez-Cara, I. Gayte-Delgado Esto ´ultimo equivale a decir que ˆv=P(−a bη),(12) donde P:L2(ω×(0, T))) 7→ Uad es la proyecci´on ortogonal habitual. En resumen, se tiene el resultado siguiente: Teorema 3 Sea (ˆv, ˆc, ˆ β)soluci´on del problema de control ´optimo (8). Entonces, existe (d, η)que verifica (10) y la relaci´on de optimalidad (12). El sistema de optimalidad (1), (10), (12) tiene soluci´on ´unica para Tsuficientemente peque˜no. Por tanto, para valores peque˜nos de T, el problema de control ´optimo (8) tambi´en tiene soluci´on ´unica. Esto justifica y hace interesante un algoritmo num´erico basado en la resoluci´on de (1), (10), (12). Nosotros hemos utilizado el siguiente algoritmo de tipo punto fijo: a) Elegir v0∈Uad. b) Dados k≥0 y vk∈Uad, (b.1) Hallar (ck, βk) soluci´on de (1) poniendo v=vk. (b.2) Hallar (dk, ηk) soluci´on de (10) poniendo en el segundo miembro ˆc=cky ˆ β=βk. (b.3) Calcular vk+1 =P(−a bηk). 4. Resultados num´ericos Wg w Ww O Figura 1: Dominio Ω del problema junto con la distribuci´on de materia blanca (Ωw) y materia gris (Ωg). La zona donde comienza a desarrollarse el tumor es O. La zona donde act´ua el control es ω. Figura 2: La triangulaci´on utilizada, con 3157 tri´angulos y 1789 v´ertices. Presentamos aqu´ı algunas experiencias num´ericas, por el momento bidimensionales. M´as concretamente, consideramos un dominio Ω que representa una secci´on coronal del cerebro. Las distribuciones de materia gris y materia blanca se muestran en la Figura 1. 4 Control de EDPs orientado a la terapia de un tumor cerebral Para la aproximaci´on num´erica se ha utilizado un esquema de Euler semi-impl´ıcito para la discretizaci´on en tiempo y un m´etodo est´andar de elementos finitos P1-Lagrange para la discretizaci´on espacial. La triangulaci´on utilizada se presenta en la Figura 2. 0 12 25 37 50 Figura 3: Isovalores de la funci´on c0, condici´on inicial en el problema (1), que representa la concentraci´on de c´elulas tumorales cuando se comienza a aplicar la terapia. 0 10 20 30 40 Figura 4: La zona rodeada por las iso-l´ıneas representa la parte del tumor “visible” por las t´ecnicas de detecci´on, al comienzo de la terapia. Es conocido, v´ease por ejemplo [1], que las t´ecnicas de detecci´on de tumores (tomograf´ıa computerizada y resonancia magn´etica) presentan un “umbral de detecci´on”, es decir, s´olo son capaces de detectar un tumor cuando la zona en que se supera una cierta densidad de c´elulas tumorales es suficientemente grande. Siguiendo tambi´en a [1], hemos supuesto que la parte del tumor “detectable” es aqu´ella con densidad mayor que 40 mil c´elulas/cm2y di´ametro superior a 3 cm. Hemos considerado un tumor virtual que inicialmente comienza a desarrollarse en una peque˜na zona del dominio, se˜nalada en la Figura 1 con la letra O, en el centro de la cual el valor inicial es de 16 mil c´elulas/cm2. Hemos calculado su evoluci´on, resolviendo la ecuaci´on correspondiente a cen (1), junto con sus condiciones de contorno e inicial, en el intervalo de tiempo de 0 a 740 d´ıas, que es cuando (seg´un lo dicho antes) empieza a ser detectable. De acuerdo con ello, hemos tomado como condici´on inicial c0=c(740) en el problema de control. En la Figura 3 aparecen los iso-valores de la funci´on c0y en la Figura 4 se puede observar la regi´on en donde c0supera el umbral de detecci´on fijado. Con esta c0inicial resolvemos el sistema de optimalidad en el intervalo de tiempo [0,360]. La elecci´on de los par´ametros en (1) ha sido llevada a cabo teniendo en cuenta el an´alisis de los trabajos [1] y [2] (v´eanse tambi´en sus referencias). Las simulaciones num´ericas realizadas corresponden a los valores que aparecen en la Tabla 1. Con objeto de comparar los efectos del control, hemos representado en las Figuras 5 y 6 la funci´on que se obtiene resolviendo (1) con v≡0. Esta funci´on nos dice c´omo ser´ıa la evoluci´on del tumor sin aplicaci´on de terapia alguna 360 d´ıas despu´es de la diagnosis. En las Figuras 7 y 8, se presenta la concentraci´on de c´elulas tumorales correspondiente al control ´optimo calculado cuando el convexo Uad est´a dado por (9) con M= 5 y, en (7), se toma a=0.02 y b=0.98. En la Figura 11 se representa la evoluci´on en el tiempo del valor medio en ωdel control ´optimo. 5 R. Echevarr´ıa, A. Doubova, E. Fern´andez-Cara, I. Gayte-Delgado Par´ametro Valor Unidad Par´ametro Valor Unidad Dg0.0002 cm2/d´ıa a10.0007 1/d´ıa Dw0.001 cm2/d´ıa b10.004 1/d´ıa µ0.5 cm2/d´ıa a20.0001 1/d´ıa b20.0003 1/d´ıa Tabla 1: Valores de los par´ametros usados en las simulaciones. 0 110 210 320 420 Figura 5: Isovalores de la funci´on concentraci´on de c´elulas tumorales 360 d´ıas despu´es de la diagnosis, sin aplicaci´on de terapia. 0 10 20 30 40 Figura 6: Parte “visible” por las t´ecnicas de detecci´on de la misma funci´on. 0 30 59 89 120 Figura 7: Concentraci´on de c´elulas tumorales correspondiente a la aplicaci´on de la terapia soluci´on del problema de control ´optimo. 0 10 20 30 40 Figura 8: Parte “visible” por las t´ecnicas de detecci´on de la misma funci´on. De forma an´aloga, en las Figuras 7, 8 y 11 se muestran los resultados obtenidos para los valores de los par´ametros a=0.08 y b=0.92 y para un control soportado por tres subdominios disjuntos. Por otro lado, desde el punto de vista pr´actico parece m´as realista la aplicaci´on de una terapia “intermitente”, es decir, con per´ıodos de medicaci´on alternados con per´ıodos de descanso. Esto corresponde, en nuestro problema de control, a la siguiente elecci´on del 6 Control de EDPs orientado a la terapia de un tumor cerebral 0 14 27 41 54 Figura 9: Concentraci´on de c´elulas tumorales en el instante final correspondiente a la aplicaci´on de un control ´optimo soportado en tres zonas disjuntas. 0 10 20 30 40 Figura 10: Parte “visible” por las t´ecnicas de detecci´on de la misma funci´on. convexo de controles admisibles: Uad ={v∈L2(ω×(0, T)) : 0 ≤v≤M, v = 0 en [τk, τk+1], k = 1,3, . . . }(13) para una determinada partici´on {τi}del intervalo [0, T]. En las Figuras 13 y 14 se muestran el control y la concentraci´on de c´elulas tumorales en el instante final correspondientes a una restricci´on del tipo anterior con M= 18. Obs´ervese que, aunque la cantidad total de medicaci´on suministrada es similar en ambas experiencias, los valores de cen este caso son m´as favorables que los que se muestran en la Figura 7, lo cual parece confirmar lo observado en la pr´actica m´edica. 5. Comentarios finales La estrategia utilizada para resolver (8) puede ser distinta. En efecto, en vez de resolver el sistema de optimalidad (1), (10), (12), puede resultar m´as conveniente aplicar un m´etodo de tipo gradiente (con proyecci´on) para la determinaci´on del control ´optimo. Por otra parte, pueden ser m´as realistas otras elecciones de Uad donde intervengan restricciones adicionales. Por ejemplo, parece razonable tomar Uad =U0 ad ∩ { v∈L∞(ω×(0, T)) : ZZ ω×(0,T ) v dx dt ≤K},(14) donde U0 ad est´a dado por (13) y K > 0. Obviamente, tambi´en es m´as realista realizar experiencias num´ericas tridimensionales en espacio. La metodolog´ıa es muy parecida, aunque las caracter´ısticas geom´etricas del problema son mucho m´as complicadas. Finalmente, indiquemos que se puede intentar ir un poco m´as all´a en la aplicaci´on de la teor´ıa de control a la terapia de tumores. En efecto, un problema m´as interesante (y 7 R. Echevarr´ıa, A. Doubova, E. Fern´andez-Cara, I. Gayte-Delgado 0 40 80 120 160 200 240 280 320 360 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 DIAS CONTROL Figura 11: Control ´optimo correspondiente a las Figuras 7 y 8. 0 40 80 120 160 200 240 280 320 360 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 DIAS CONTROL Figura 12: Control ´optimo correspondiente a las Figuras 9 y 10, soportado por tres zonas disjuntas. Cada curva corresponde a una de las zonas. 012 48 60 96 108 144 156 192 204 240 252 288 300 336 348 360 0 2 4 6 8 10 12 14 16 18 20 DÍAS CONTROL Figura 13: Control intermitente. 0 21 43 64 85 Figura 14: Iso-valores de la funci´on c(T) asociada al control intermitente. m´as dif´ıcil) aparece cuando el objetivo es, en vez de minimizar Jen Uad, elegir ven Uad de manera que la soluci´on de (1) verifique c(T)∈Bd(T),(15) donde Bd(T) es un conjunto de estados finales deseados fijado. Esto es un problema de controlabilidad con restricciones. Estas cuestiones ser´an tratadas en varios trabajos de futura aparici´on. Agradecimientos El trabajo de los autores ha sido parcialmente financiado por Ministerio de Educaci´on y Ciencia, Proyecto MTM2006-07932 y por Junta de Andaluc´ıa, Proyecto FQM 520. Para la realizaci´on de las gr´aficas de este trabajo se ha utilizado el software: Scilab (Copyright c °1989-2005. INRIA ENPC) y Medit (Pascal Frey, Laboratoire Jacques Louis 8 Control de EDPs orientado a la terapia de un tumor cerebral Lions de l’Universit´e Pierre et Marie CURIE). Referencias [1] K. R. Swanson, E. C. Alvord Jr, and J. D. Murray, A quantitative model for differential motility of gliomas in grey and white matter, Cell Prolif., 33, 317–329, 2000. [2] K. R. Swanson, E. C. Alvord Jr, and J. D. Murray, Quantifying Efficacy of Chemotherapy of Brain Tumors (Gliomas) with Homogeneous and Heterogeneous Drug Delivery, Acta Biotheoretica, 50(4): 223-237, 2002. [3] I.V. Girsanov, Lectures on mathematical theory of extremum problem, Lectures notes in Economics and Mathematical Systems, 67, Springer-Verlag, Berl´ın 1972. 9