Repositorio Institucional de Documentos

Modelos matemáticos en oncología Simulación numérica Marta Gómez Gómez Trabajo de fin del grado de Matemáticas Universidad de Zaragoza Prólogo “En los últimos años las ciencias de la vida se han transformado, y están empezando a ser, en cierto modo, una ciencia exacta. Se ha abierto un mundo nuevo en el estudio de los seres vivos, un universo con una cantidad de datos enorme. Hay que ponerle forma a eso, y las matemáticas son la herramienta para hacerlo”. Fernando Giráldez, catedrático de Biología del Desarrollo de la Universidad de Pompeu Fabra. En principio las matemáticas no precisan de la biología pero sí recíprocamente, siendo éste el origen de la biología matemática, la cual será claramente líder de la ciencia en un futuro previsible. Aunque no esté específicamente definida, la biología matemática es una ciencia de rápido crecimiento y bien reconocida. Los mejores modelos, muestran cómo funciona un proceso y predicen cómo seguirá. Un modelo matemático es algo más que un conjunto de ecuaciones que hay que analizar para que luego un biólogo pueda interpretar sus resultados. Un modelo matemático se encuentra determinado por el problema biológico que se quiere resolver, la pregunta biológica que se quiere responder. No obstante, construir un modelo matemático y analizarlo no es hacer biología. La exploración matemática de un fenómeno biológico, por ejemplo, mediante un sistema de ecuaciones, no es equivalente a la construcción de una teoría biológica. La aplicación de las matemáticas a la biología tiene esencialmente dos caminos: por un lado la aplicación rutinaria de técnicas conocidas; por otro, el desarrollo de nuevos métodos necesarios para el análisis de sistemas biológicos. La aportación de las matemáticas, en concreto a la lucha contra el cáncer, se hace sobre todo a través de modelos y programas que simulan desde cómo crece un tumor a qué efecto tiene sobre un paciente determinada terapia. Los primeros modelos empezaron a desarrollarse en los años 70 y algunos se emplean hoy ya rutinariamente en los hospitales. Se han desarrollo modelos provenientes de diversos campos como la matemática aplicada, la estadística, la mecánica de fluidos, la ciencia computacional, etc. En particular, los modelos de ecuaciones en derivadas parciales, son herramientas muy usadas en el estudio del crecimiento de tumores y la forma en que se difunden sobre los tejidos que los rodean. Sin embargo, hoy en día se carece de un modelo que proporcione una predicción y caracterización del comportamiento para el crecimiento de tumores cancerosos en sus múltiples formas y para cualquier tipo de población, teniendo en cuenta que los modelos existentes funcionan bajo condiciones ideales y con poblaciones específicas. Así pues, a lo largo de este trabajo se tratará de introducir un primer modelo matemático sencillo, con el fin de producir un acercamiento al lector con todo este campo. Se busca una ecuación que modele en concreto tumores cerebrales, dando lugar a las ecuaciones de reacción-difusión. Posteriormente, se hará un estudio de este tipo de ecuaciones de forma matemática y de métodos numéricos para su resolución, como el método de diferencias finitas. I II Capítulo 0. Prólogo Finalmente, veremos una posible aplicación de todo el estudio, mediante la simulación numérica con Matlab de casos con pacientes virtuales donde nuestro fin será averiguar cuál es la mejor forma de distribuir dosis de radioterapia en función de distintos aspectos o características. No obstante, no se trata de que los matemáticos o los ordenadores vayan a suplantar a los oncólogos, sino de orientar los pasos de los médicos y ayudarles a probar primero lo que parece más efectivo y por tanto, produce mejores resultados. Modelos matemáticos en oncología Simulación numérica Summary Mathematical biology is a science of rapid growth and well recognized, although not clearly defined. In this field, a good mathematical model is determined by the biological problem to be solved, i.e, biological question you want to answer. The contribution of mathematics to cancer, was made primarily through models or programs that simulate from how a tumor grows to what effect has a therapy given on a patient. In particular, models of partial differential equations are widely used tools in the study of tumor growth and in the way they are spread over the surrounding tissues. Chapter 1: Reaction-diffusion equations in oncology Cancer is one disease in which cells divide without control and can invade other tissues. The cells of our body grow and divide in a controlled way to produce more as needed. When cells grow old or get damaged, they die, and they are replaced by new cells. However, sometimes this orderly process goes wrong. The genetic material (DNA) of the cell may be damaged or altered, which produces mutations that affect normal growth and division of cells. When this happens, cells do not die when they should, and form new cells that the body does not need, creating a mass of tissue which is called tumor. In this work, we seek to analyze a simple model for brain tumors, in particular, gliomas, which arise from the glial cells; that is, cells of the central nervous system, which mainly have the function of support of neurons and are involved in brain information processing in the body. The type of equation that models this type of tumor are reaction-diffusion equations. We suppose a population u(x,t). It is in a set Ω⊆Rnwhit Ωcompact and denote by J(x,t)∈Rn the particles flow in and out of Ω. Thus, the conservation equation, says that the rate of change of the density u(x,t)in Ωis equal to the rate of change of the flow of material through ∂Ωplus material created in Ω:∂ ∂tZΩ u(x,t)dΩ=−Z∂Ω J(x,t)−→ n d∂Ω+ZΩ f(u(x,t))dΩ,(1) where f(u(x,t)) describes the birth rate, death, etc. By using the divergence theorem and substituting in (1) the following equation is finally obtained: ∂u ∂t=O(dOu)+ f(u),(2) where dis the so-called diffusion constant. This type of equation is called reaction-diffusion equation. In particular, we work with the Fisher-Kolmogorov equation, which we can obtain by deriving the following logistic model for the growth of a population: (u0(t) = ρu(t)(A−u(t)), u(0) = f0,f0∈(0,A),(3) III IV Capítulo 0. Summary where u=u(t)is the density of population, ρ>0 is the growth rate and A>0 is the capacity of the environment. We assume that the population is uniformly distributed over an area all time, and thus by deriving (3), the Fisher equation is obtained: ut=duxx +ρu(A−u),(4) where dis the diffusion coefficient and u=u(x,t)is again the density of population. Normally, we study this equation together with boundary conditions of Neumann type, ux(0,t) = ux(L,t) = 0,(5) where Lis the domain length. Therefore, in the modeling of the most basic aspect of gliomas, the Fisher equation describes the spatio-temporal dynamics of a density of cancer cells that can migrate and proliferate. Defining f(u) = ρu(A−u)with ρthe parameter of proliferation, we rewrite (4) as: ut=duxx +f(u),(6) where f(u)models the proliferation (birth and death of tumor cells), and duxx models diffusion. We could model other characteristics of tumor cells, however, equation (6) models two of their most important features: proliferation and diffusion. Chapter 2: Finite difference method The method of finite difference is studied for the solution of the reaction-diffusion equation. This method consists in the approximate solution of PDEs by using a discretization on a mesh. First, the study of this method is done by considering the heat equation: ut=uxx,(7) with initial condition u(x,0) = η(x)and boundary conditions of Dirichlet type, (u(0,t) = g0(t), u(1,t) = g1(t),(8) for t>0, 0 ≤x≤1. Generally, a set of finite difference equations is obtained on a grid with discrete points (xi,tj) where: xi=ih 0≤i≤N, tj=jk 0≤n≤M, where hand kare the spatial and temporal discretization steps respectively. So, denote byUj i≈u(xi,tj) the numerical approximation of uat discrete point xi, in the time tj. We distinguish between explicit and implicit methods. On the one hand, we obtain an explicit method considering the following natural discretization of (7): Uj+1 i−Uj i k=1 h2(Uj i−1−2Uj i+Uj i+1).(9) Modelos matemáticos en oncología Simulación numérica V On the other hand, the most usual implicit method is the method of Crank-Nicolson, which can be written as: Uj+1 i−Uj i k=1 2h2(Uj i−1−2Uj i+Uj i+1+Uj+1 i−1−2Uj+1 i+Uj+1 i+1).(10) We can also consider the method of lines, which discretizes the PDE only in space, and solve the resulting system of ODEs by any of the known methods, such as Euler method. For example, for the heat equation (7), it can be discretized in space in the following way: U0 i(t) = 1 h2(Ui−1(t)−2Ui(t)+Ui+1(t)) para i=1,2,...,N-1.(11) Stability and convergence To study the stability and convergence of the method we start by defining the local truncation error, which is obtained by inserting the exact solution u(x,t)of the PDE in the finite difference equation. In particular, for example for explicit methods, it is obtained that: τ(x,t) = u(x,t+k)−u(x,t) k−1 h2(u(x−h,t)−2u(x,t)+u(x+h,t)), refering to τjor kτjkas the local truncation error. Although we do not know u(x,t), in general, we assume that it is smooth and then use Taylor series expansions. It is said that a method is consistent if τ(x,t)−→ 0 when k,h−→ 0. It is said that the difference scheme is consistent of order (p,q)for the partial differential equation given if: kτjk=O(hp)+O(kq). The explicit methods are second order accurate in space and first order accurate in time, as the local truncation error is O(h2+k). However, the method of Crank-Nicolson is second-order accurate in both space and time, that is O(h2+k2). To study the stability, we assume that Uj+1can be written as: Uj+1=QU j,j≥0,(12) and say that a difference scheme of the form (12) is stable with respect to the norm k.kif there exist positive constants h0and k0and non-negative constants αand βsuch that: kUj+1k ≤ αeβtkU0k,(13) for 0 ≤t = (j+1)k; 0 <h≤h0and 0 <k≤k0. In conclusion, the most common approach for the convergence of a finite difference scheme is through the concepts of consistency, stability and Lax theorem. We can state the Lax theorem. It says that if a system is consistent of order (p,q)in the norm k.k for an initial value well-posed problem and it is stable with respect to norm k.k, then it is convergent of order (p,q)with respect to norm k.k. Autora: Marta Gómez Gómez VI Capítulo 0. Summary Difference scheme for the Fisher equation In this chapter, we also study the method for the Fisher equation. Using an explicit scheme, you can obtain the following discretization: Uj+1 i=Uj i+dk h2(Uj i+1−2Uj i+Uj i−1)+kρ(A−Uj i)Uj i, with boundary conditions of Neumann type. We know some properties of the solution of the Fisher equation, particularly with ρ=d=A=1. We suppose u(x,t), satisfying u,ux,uxx,ut∈C([0,1]×[0,∞]), is the solution of a problem of the form:        ut=uxx +u(1−u), ux(0,t) = ux(1,t) = 0, u(x,0) = f(x). (14) Then, if f(x)satisfies 0 <ε≤f(x)≤1+ε, it is fulfilled that: 0<ε≤u(x,t)≤1+ε, for all x∈[0,1],t≥0. For the Neumann problem (14), it is easy to see that the solution can tend to infinity in a finite time for some specific values of gwith g(u) = u(1−u). Note that if the initial condition fis constant, i.e, f(x) = f0,(15) for all x∈[0,1], then; u(x,t) = v(t),x∈[0,1],t>0, where vis the solution of: v0(t) = g(v),v(0) = f0. Therefore the solution of (14) is given by the solution of an ordinary differential equation which tends to infinity in a finite time for some functions g. IMEX method Difference schemes can be extended in several dimensions, for example, for the Fisher-Kolmogorov equation: ut=d(uxx +uyy)+ρu(A−u), with boundary conditions of Neumann type. To treat the non-linear part of the equation, the IMEX method will be used. This method involves treating implicitly the diffusion term and explicitly the reaction term. So, in each time step, we have to solve a system of linear equations. In particular, to an internal node (i,j)of the grid, the difference equation is: Uj+1 i−Uj i k−dUj+1 i+1−2Uj+1 i+Uj+1 i−1 h2−d Uj+1 i+N+1−2Uj+1 i+Uj+1 i−(N+1) h2−ρUj i(A−Uj i) = 0. Modelos matemáticos en oncología Simulación numérica VII Chapter 3: Numerical simulation Numerical solution: explicit method for heat equation First, we use the explicit finite difference method for solving the heat equation;        ut=uxx,x∈(0,1),t<0, u(0,t) = u(1,t) = 0, u(x,0) = sin(2πx). (16) If solving the heat equation by separate variables, the following solution is obtained: u(x,t) = e−4π2tsin(2πx). The aim of this section is to prove that with a finer mesh, we get better numerical approximations. Therefore, first we compute the numerical aprroximation with h=1 6,k=1 80 and after with h=1 20, k=1 800. Also, one can make a study of the errors given by the method. We begin by calculating the error in the two cases previously studied to corroborate more accurately than the numerical solution in the second case is much closer to the exact solution. This study yields that 0.0013 <0.0112; i.e, the error with the values of hand ktaken in the second case is much smaller than the error for the first case. Numerical solution: explicit method for Fisher-Kolmogorov equation A second aim is to apply the explicit finite difference method to other equations, such as the Fisher-Kolmogorov (FK) equation and to study the behaviour for different values of the diffusion and proliferation parameters. The spatial interval where it is solved is [a,b]choosing aand bappropriately, and remember that the FK equation is normally considered with boundary conditions of Neumann type. So, we choose as integration interval [−L/2,L/2], taking L=6 cm and we use as an initial condition the following function: u0(x) = U0e−x2 σ2. This implies that the cancer cells are distributed in the tissue thereby, where the constant (amplitude) 0 <U0<1 takes values in the range U0∈[10−3,10−1]. The diffusion and proliferation parameters are considered as d∈[1,103]mm2/year and ρ∈[10−1,102]year −1. The time interval we want to explore is t∈[0,T]with T∈[10−1,10]years. We choose the remaining parameter σso that the initial spatial distribution of the tumour is in a sufficiently small region of the interval [−L/2,L/2]. In this way, we will calculate the solution varying independently the values of dand ρ. We have implemented in Matlab the explicit finite difference method for solving FK equation and we are going to discuss the different cases obtained. First, we assume a fixed proliferation parameter given by ρ=10 and we analyze the differences that appear if the diffusion parameter is very small d=1 or very high d=1000. It is noted, on the one hand, that for a proliferation parameter not excessively large, regardless of what the tumor is disseminated, in a period of one year both tumors will occupy maximum tissue. On the other hand, we see how for the same proliferation parameter, with a very large diffusion, the number of cells is expanded in a short period of time, covering the entire line of integration; and with a very small diffusion, it is not the covered even when ureaches its maximum value after a year. Autora: Marta Gómez Gómez Capítulo 1 Ecuaciones de reacción-difusión en oncología 1.1. ¿Qué es el cáncer? Definición 1.1.1. Se denomina como cáncer al término usado para enfermedades en las que las células se dividen sin control y pueden invadir otros tejidos. No se trata únicamente de una enfermedad, sino de muchas, ya que existen más de 100 tipos diferentes. No obstante, se pueden agrupar en las siguientes categorías principales: -Carcinoma: cáncer que empieza en la piel o en tejidos que cubren órganos internos. -Sarcoma: cáncer que empieza en hueso, cartílago, grasa, músculo o vasos sanguíneos. -Leucemia: cáncer que empieza en el tejido en el que se forma la sangre, como la médula ósea. -Linfoma y mieloma: cáncer que empieza en las células del sistema inmunitario. -Cáncer del sistema nervioso central: cáncer que empieza en los tejidos del cerebro y la médula espinal. Todos empiezan en las células. Para entender bien qué es el cáncer, es importante saber lo que sucede cuando las células normales se convierten en cancerosas. Las células de nuestro cuerpo crecen y se dividen de una forma controlada para producir más células según sea necesario. Cuando las células envejecen o se dañan, mueren y son reemplazadas por células nuevas. Sin embargo, algunas veces este proceso ordenado se descontrola. El material genético (ADN) de una célula puede dañarse o alterarse, lo cual produce mutaciones que afectan al crecimiento y a la 1 2Capítulo 1. Ecuaciones de reacción-difusión en oncología división normal de las células. Cuando esto sucede, las células no mueren cuando deberían morir, y se forman células nuevas a pesar de que el cuerpo no las necesita. Las células que sobran forman una masa de tejido a la que se le llama tumor. No obstante, no todos los tumores son cancerosos; existen tumores benignos y malignos. -Tumores benignos: pueden estirparse y en la mayoría de los casos no vuelven a aparecer. Las células no se diseminan a otras partes del cuerpo. -Tumores malignos: las células pueden invadir tejidos cercanos y diseminarse a otras partes del cuerpo. Al proceso en el que el cáncer se disemina a otra parte del cuerpo, se le llama metástasis. 1.1.1. Tumores cerebrales Este trabajo muestra modelos para tumores cerebrales, los cuales abarcan cualquier tumor que se inicie en el cerebro. Estos tumores se pueden originar a partir de las células cerebrales, las membranas alrededor del cerebro, nervios o glándulas. Los tumores pueden destruir directamente células cerebrales o provocarles daño produciendo inflamación, ejerciendo presión sobre otras partes del cerebro e incrementando la presión intracraneal. Los tumores cerebrales pueden ocurrir a cualquier edad, pero muchos de ellos son más comunes en un grupo de edad en particular. Por ejemplo, en los adultos, los gliomas y los meningiomas son los más comunes. Los meningiomas son muy frecuentes y por lo general benignos. Se presentan en el tejido aracnoideo de las meninges1y se adhieren a la duramadre. No obstante, pueden causar serias complicaciones e incluso la muerte debido a su tamaño y localización. Los gliomas, por contra, a pesar de que apenas metastatizan, rara vez se pueden curar. Surgen a partir de las células gliales, es decir, células del sistema nervioso central que desempeñan de forma principal, la función de soporte de las neuronas e intervienen activamente en el procesamiento cerebral de la información en el organismo. Son clasificados de acuerdo a su grado, y del cual depende que se augure un mejor o un peor pronóstico, siendo éste por lo general malo para pacientes con gliomas de alto grado. El tratamiento para tumores cerebrales puede involucrar cirugía, radioterapia y quimioterapia. No obstante, el tratamiento depende en cada caso del tamaño, del tipo de tumor y de la salud del paciente en general; y no siempre tiene como objetivo la cura. En algunas ocasiones en las que se sabe a priori que no tiene cura, únicamente se busca el alivio de los síntomas o la mejora de la actividad cerebral. Así pues, el objetivo es encontrar una ecuación en derivadas parciales que modele en concreto los gliomas y resolverla mediante un método numérico adecuado. 1.2. Ecuaciones de reacción-difusión Una clase importante de ecuaciones en derivadas parciales son las ecuaciones de reacción-difusión, para las cuales las variables independientes son el tiempo t y las variables espaciales x. Las ecuaciones de reacción-difusión implican la combinación de dos procesos diferentes: reacción y difusión. Se comienza definiendo ambos conceptos. 1Las meninges son las membranas del tejido conectivo que cubren todo el sistema nervioso central, siendo la aracnoides la meninge intermedia. Modelos matemáticos en oncología Simulación numérica 1.2. Ecuaciones de reacción-difusión 3 Definición 1.2.1. Se entiende como difusión a la tendencia de las moléculas a moverse desde zonas de alta concentración hacia zonas de baja concentración. Definición 1.2.2. Se entiende por reacción al cambio de estado de las partículas, debido por ejemplo a interacciones o de manera espontánea. 1.2.1. Modelo logístico En primer lugar, se recobra un modelo de logística para el crecimiento de una población, del cual se deducirá posteriormente la ecuación de Fisher-Kolmogorov como caso particular de ecuación de reacción-difusión. Este modelo, manifiesta que el crecimiento de la población frente a los recursos limitados es gobernado por la siguiente ecuación: (u0(t) = ρu(t)(A−u(t)), u(0) = f0,f0∈(0,A),(1.1) donde u=u(t)es la densidad de población, ρ>0 es el índice de crecimiento y A>0 es la llamada capacidad de carga del medio ambiente. El modelo nos muestra que para pequeñas poblaciones, se logra un crecimiento exponencial gobernado por: u0(t)≈ρAu(t). Sin embargo, si uaumenta, el término −ρu2comienza a ser significativo, el crecimiento se ralentiza y la población alcanza poco a poco la capacidad de carga del medio ambiente. El problema (1.1) puede ser resuelto analíticamente, mediante el método de variables separadas. du dt =ρu(A−u),u6=0,A, du u(A−u)=ρdt. Es decir, se ha de resolver: Zdu u(u−A)=Z−ρdt, lo cual da lugar a: lnÇu−A uå=−ρAt +C,con C constante, o lo que es lo mismo, despejando u: u(t) = A 1−Ke−ρAt ,con K constante. Aplicando la condición inicial (u(0) = f0), se obtiene la solución buscada. u(t) = A f0 f0+(A−f0)e−ρAt ,t≥0. Notar que u=Aes la solución asintótica cuando t−→ ∞para cualquier valor inicial f0>0. Autora: Marta Gómez Gómez 4Capítulo 1. Ecuaciones de reacción-difusión en oncología Por ejemplo con A=ρ=1, se puede observar gráficamente la solución para algunos valores de f0en la Figura 1.1, donde se aprecia como todas las hipotéticas poblaciones tienden finalmente a A. Figura 1.1: Solución del modelo logístico del crecimiento de una población para distintos valores de la condición inicial f0 1.2.2. Deducción de las ecuaciones de reacción-difusión Se supone una población u(x,t)en un conjunto Ω⊆Rncon Ωabierto y se denota por J(x,t)∈Rn el flujo de partículas que entran y salen de Ω. La ecuación de conservación, nos dice que la tasa de cambio de la densidad u(x,t)en Ωes igual a la tasa de cambio del flujo del material a través de ∂Ωmás el material creado en Ω. Es decir, de forma esquemática: CAMBIO EN Ω−→ FLUJO A TRAVÉS DE ∂Ω + CAMBIO EN LA TASA DE NACIMIENTO Y MUERTE EN Ω Escrito de forma matemática, esto es: ∂ ∂tZΩ u(x,t)dΩ=−Z∂Ω J(x,t)·−→ n d∂Ω+ZΩ f(u(x,t))dΩ,(1.2) donde f(u(x,t)) describe la tasa de nacimiento, muerte, etc. Notar, que consideramos que al final la tasa media del flujo entra. Teorema 1.2.3. Teorema de Gauss de la divergencia Sea Ωun abierto simple de R2y S =∂Ωsu borde, orientado con la norma exterior unitaria −→ n . Sea F :Ω−→ R2un campo vectorial de clase C1(Ω). Entonces: ZΩ divFdΩ=ZS F·−→ n dS. De esta manera, haciendo uso del teorema de la divergencia y sustituyendo en (1.2), se obtiene: ZΩ (∂u ∂t−f(u)+divJ)dΩ=0. Por la Ley de Fick, el flujo es proporcional al gradiente de la concentración del material, donde la constante de proporcionalidad es el coeficiente de difusión d, es decir, J=−dOu. Observar que el signo negativo es debido al hecho de que va de mayor densidad a menor densidad. Modelos matemáticos en oncología Simulación numérica 1.3. Ecuación de Fisher-Kolmogorov en oncología 5 Así, se llega a la ecuación: ∂u ∂t=O(dOu)+ f(u).(1.3) Definición 1.2.4. A una ecuación de la forma (1.3), se le suele llamar ecuación de reacción-difusión. Considerando condiciones iniciales y de contorno, se tiene por ejemplo, el siguiente problema:          ∂u ∂t=dMu+f(u),x∈Ω, u(x,t) = 0,x∈∂Ω,t>0, u(x,0) = g(x),x∈Ω, (1.4) donde fpuede depender de forma no lineal de u. Se puede ver con más detalle en [M] (Volumen I: capítulo 11). 1.3. Ecuación de Fisher-Kolmogorov en oncología En el modelo logístico (1.4), se supone que la variación espacial de la densidad de la población es de poca importancia para el crecimiento de ésta. Es decir, se asume que la población se distribuye de manera uniforme sobre un área todo el tiempo. No obstante, en poblaciones reales, esta suposición es a menudo bastante dudosa. Por lo tanto, se trabaja con la siguiente ecuación, obtenida de derivar el modelo logístico (1.1): ut=duxx +ρu(A−u),(1.5) donde d es el coeficiente de difusión y u=u(x,t)es nuevamente la densidad de población. Definición 1.3.1. La ecuación del modelo de crecimiento de una población (1.5), se suele llamar ecuación de Fisher. La introducción de un término de difusión, conduce a una ecuación diferencial parcial que en contraste con la ecuación diferencial ordinaria (1.1), no puede ser generalmente resuelta analíticamente. Como ya se ha mencionado, el término duxx modela la difusión de la población. Remarcar el hecho de que, términos similares, surgen en muchas aplicaciones donde se quiere capturar la tendencia de la naturaleza para suavizar las cosas. Normalmente, la ecuación de Fisher (1.5) se estudia en conjunto con condiciones de frontera del tipo Neumann, ux(0,t) = ux(L,t) = 0,(1.6) donde L denota la longitud del dominio. La razón por la que se eligen estas condiciones de frontera, es que se asume que el área es cerrada, así que no hay migración a través del dominio. La ecuación de Fisher-Kolmogorov, constituye uno de los ejemplos más elementales de ecuación de reacción-difusión no lineal y uno de sus usos es la modelización de tumores cerebrales (gliomas malignos), principalmente por la ausencia de metástasis, lo que justifica las condiciones de frontera de tipo Neumann. No obstante, como en todos los tumores, los aspectos biológicos y clínicos de los gliomas, son complejos y los detalles de su crecimiento espacio-temporal todavía no se entienden bien. Para construir estos modelos se tienen que hacer algunas suposiciones previas. El modelo teórico más simple incluye sólo el número total de células del tumor, asumiendo normalmente que el tumor tiene un crecimiento exponencial. Autora: Marta Gómez Gómez 6Capítulo 1. Ecuaciones de reacción-difusión en oncología Estos modelos, no tienen en cuenta la disposición espacial de las células en un lugar anatómico en concreto o la extensión espacial de las células cancerosas; aspectos que son cruciales en la estimación del crecimiento del tumor, ya que determinan la capacidad de invasión y el aparente borde del tumor. Así pues, la ausencia de un modelo simple que explique con exactitud el crecimiento de los gliomas humanos, es lo que hace díficil explicar por qué los resultados tras una extirpación quirúrgica son tan decepcionantes. Por lo tanto, se modelan únicamente los aspectos más básicos de los gliomas, usando la ecuación de Fisher-Kolmogorov para describir la dinámica espacio-temporal de una densidad de células cancerígenas que pueden migrar y proliferar. Definiendo f(u) = ρu(A−u)con ρel parámetro de proliferación, se puede reescribir (1.5) como: ut=duxx +f(u),(1.7) donde f(u)modela la proliferación, es decir, la tasa de nacimiento y muerte de las células tumorales yduxx viene del hecho de que cuando las células tumorales han crecido lo suficiente, migran, se difunden. Notar de nuevo, que se podrían modelizar muchas otras características de las células tumorales. No obstante, la ecuación (1.7) modela dos de las más importantes: proliferación ydifusión. En el siguiente capítulo se estudiará un método numérico para la resolución de dicha ecuación. Modelos matemáticos en oncología Simulación numérica Capítulo 2 Método de diferencias finitas 2.1. Métodos explícitos e implícitos En este capítulo, se va a estudiar el método de diferencias finitas para aproximar ecuaciones en derivadas parciales dependientes del tiempo. Se puede ver un estudio completo del método en [L] (capítulo 13). Definición 2.1.1. El método de diferencias finitas es un método de carácter general que permite la resolución aproximada de ecuaciones diferenciales en derivadas parciales definidas en recintos finitos, mediante la discretización del recinto del plano en el que se quiere resolver la ecuación con una malla, por conveniencia cuadrada. Se realiza el estudio del método de diferencias finitas para la ecuación del calor unidimensional en el dominio Ω= (0,1), la cual es un clásico ejemplo de ecuación parabólica: ut=uxx.(2.1) Junto con esta ecuación, se necesitan condiciones iniciales, u(x,0) = η(x),(2.2) y condiciones de contorno, como por ejemplo, condiciones de tipo Dirichlet, u(0,t) = g0(t),u(1,t) = g1(t),(2.3) para t >0, si 0 ≤x≤1. En la práctica, por lo general, se aplica un conjunto de ecuaciones en diferencias finitas en una cuadrícula discreta con puntos discretos (xi,tn) donde dados N,M∈Z+: xi=ih 0≤i≤N, tj=jk 0≤j≤M, con hyklos pasos de discretización espacial y temporal respectivamente. De esta manera, se denota por Uj i≈u(xi,tj)la aproximación numérica en el punto discreto xien el tiempo tj. 7 8Capítulo 2. Método de diferencias finitas 2.1.1. Método explícito Como un primer ejemplo, se puede considerar la siguiente discretización natural de (2.1) para 1≤i≤N−1, j≥0: Uj+1 i−Uj i k=1 h2(Uj i−1−2Uj i+Uj i+1).(2.4) Este ejemplo se trata de un método explícito ya que se puede calcular cada Uj+1 iexplícitamente en términos de los datos de la etapa de tiempo anterior: Uj+1 i=Uj i+k h2(Uj i−1−2Uj i+Uj i+1).(2.5) La Figura 2.1 muestra la molécula para este método. Se trata de un método de un paso en el tiempo, llamado método de dos niveles en el contexto de las EDPS, ya que consta de la solución en dos niveles de tiempo diferentes. Figura 2.1: Molécula para el método (2.5) No siempre se obtienen buenos resultados con el método explícito, debido a la existencia de una condición de estabilidad que se cumple para elecciones apropiadas de los pasos de discretización espacial y temporal. En particular, para la ecuación del calor unidimensional se tiene la restricción: k (h)2≤1 2, de forma que si se eligen unos pasos de discretización que no la cumplan, puede ocurrir que los errores obtenidos en un paso de tiempo sean más grandes que en el paso anterior. Se consideran por ello métodos implícitos para obtener métodos más estables, ya que no tienen ninguna restricción en los tamaños de los pasos de discretización. No obstante, cuentan con el inconveniente de ser más costosos (en cada paso de tiempo) desde un punto de vista computacional, ya que hay que resolver un sistema de ecuaciones en cada nivel temporal. 2.1.2. Método de Crank-Nikolson Por otro lado, como ejemplo clásico de método implícito se tiene el método de Euler: Uj+1 i−Uj i k=1 h2(Uj+1 i−1−2Uj+1 i+Uj+1 i+1).(2.6) Modelos matemáticos en oncología Simulación numérica 2.1. Métodos explícitos e implícitos 9 No obstante, otro método mucho más usual en la práctica, por sus beneficios en cuanto a condiciones de orden que se estudiarán más adelante es el método de Crank-Nikolson, donde se tiene la siguiente discretización: Uj+1 i−Uj i k=1 2(D2Uj i+D2Uj+1 i) = 1 2h2(Uj i−1−2Uj i+Uj i+1+Uj+1 i−1−2Uj+1 i+Uj+1 i+1),(2.7) que puede ser reescrito como: −rU j+1 i−1+(1+2r)Uj+1 i−rU j+1 i+1=rU j i+1+(1−2r)Uj i+rU j i+1,(2.8) donde r=k 2h2. Su correspondiente molécula es: Figura 2.2: Molécula para el método (2.8) El método de Crank-Nikolson es un método implícito y da lugar a un sistema tridiagonal de ecuaciones a resolver en cada paso de tiempo. En forma matricial, se tiene 1:           (1+2r)−r −r(1+2r)−r −r(1+2r)−r ......... −r(1+2r)−r −r(1+2r)                      Uj+1 1 Uj+1 2 Uj+1 3. . . Uj+1 N−2 Uj+1 N−1            =            r(g0(tj)+g0(tj+1))+(1−2r)Uj 1+rU j 2 rU j 1+(1−2r)Uj 2+rU j 3 rU j 2+(1−2r)Uj 3+rU j 4 . . . rU j N−3+(1−2r)Uj N−2+rU j N−1 rU j N−2+(1−2r)Uj N−1+r(g1(tj)+g1(tj+1))            . Un sistema tridiagonal de (N−1) ecuaciones puede ser resuelto con un coste computacional de O(N)por el algoritmo de Thomas. 1Notar que las condiciones de contorno (2.3) entran en esas ecuaciones. Autora: Marta Gómez Gómez 16 Capítulo 2. Método de diferencias finitas Se sigue del Teorema 2.4.1 que u(x,t)≥ε>0,∀x∈[0,1],t≥0 y como consecuencia se tiene: E0(t)≤ −2εZ1 0(1−u(x,t))2dx =−2εE(t). Por tanto, la desigualdad de Gronwall3implica que: E(t)≤e−2εtE(0), de donde se obtiene el siguiente resultado. Teorema 2.4.2. Sea u(x,t)solución del problema (2.21) con f (x)satisfaciendo 0<ε≤f(x)≤1+ε, ∀x∈[0,1]. Entonces la solución asintótica de u(x,t)se aproxima a u =1en el sentido de que: Z1 0(u(x,t)−1)2dx ≤e−2εtZ1 0(1−f(x))2dx, para t ≥0. Para el problema de Neumann (2.21), no es difícil ver que la solución puede tender a infinito en un tiempo finito para algunos valores concretos de gcon g(u) = u(1−u). Notar que si la condición incial fes constante, por ejemplo, f(x) = f0,(2.22) para todo x∈[0,1], entonces: u(x,t) = v(t),x∈[0,1],t>0, donde ves la solucion de: (v0(t) = g(v), v(0) = f0. Por lo tanto, la solución de (2.21) viene dada por la solución de una ecuación diferencial ordinaria, la cual es conocido que tiende a infinito en un tiempo finito para algunas funciones g. Sea, por ejemplo, g(v) = v3,f0>0, entonces la solución será: v(t) = f0 »1−2t f 2 0 , la cual cumple que: v(t)−→ ∞cuando t−→ 1 2f2 0 . En conclusión, la solución de (2.21) tiende a infinito en un tiempo finito con unas condiciones iniciales que satisfacen (2.22) cuando g(u) = u3yf0>0. 3Sea Iun intervalo de la forma [a,b]con a<b. Si ues diferenciable en Iy satisface u0(t)≤β(t)u(t)entonces se cumple u(t)≤u(a)exp(Rt aβ(s)ds). Modelos matemáticos en oncología Simulación numérica 2.5. Extensión del método a dos dimensiones para la ecuación de Fisher 17 2.5. Extensión del método a dos dimensiones para la ecuación de Fisher Los esquemas en diferencias finitas explicados en las secciones anteriores, se pueden extender fácilmente a varias dimensiones. En particular, en el siguiente capítulo se presentarán resultados para la ecuación de Fisher en dos dimensiones: ut=d(uxx +uyy)+ρu(A−u), con condiciones de frontera de tipo Neumann. Se considera una malla uniforme para un cuadrado y una discretización uniforme en tiempo. Para tratar la parte no lineal de la ecuación se usará el método IMEX, que consiste en tratar implícitamente el término de difusión y explícitamente el de reacción. De esta forma, en cada paso de tiempo, tendremos que resolver un sistema de ecuaciones lineales. En particular, para un nodo interior (i,j)de la malla la ecuación en diferencias es: Uj+1 i−Uj i k−dUj+1 i+1−2Uj+1 i+Uj+1 i−1 h2−d Uj+1 i+N+1−2Uj+1 i+Uj+1 i−(N+1) h2−ρUj i(A−Uj i) = 0, o lo que es lo mismo: −dk h2Uj+1 i+1−dk h2Uj+1 i−1−dk h2Uj+1 i+N+1−dk h2Uj+1 i−(N+1)+(4dk h2+1)Uj+1 i=Uj i+kρUj i(A−Uj i). De esta manera, resulta un sistema de ecuaciones lineales (ya que la parte no lineal de la ecuación se ha discretizado explícitamente) Ax =bdonde bincluye la parte no lineal y la solución del sistema es la solución ubuscada. Este tipo de métodos pueden verse en [HV] (capítulo 4). Autora: Marta Gómez Gómez Capítulo 3 Simulación numérica 3.1. Solución numérica: método explícito 3.1.1. Ecuación del calor Un primer objetivo de este capítulo, es comparar la solucion numérica y analítica de la ecuación del calor; para la cual se ha implementado el método de diferencias finitas explícito presentado en el capítulo anterior. Se considera el siguiente problema parabólico:        ut=uxx,x∈(0,1),t<0, u(0,t) = u(1,t) = 0, u(x,0) = sin(2πx). (3.1) Se puede obtener fácilmente la solución analítica de este problema por el método de separación de variables. Para ello se buscan soluciones u(x,t)de la forma u(x,t) = X(x)T(t). Insertando u(x,t) = X(x)T(t)en (3.1), se obtiene: X(x)T0(t) = X00(x)T(t). Dividiendo por X(x)T(t): T0(t) T(t)=X00(x) X(x)=−λ Considerando las condiciones de contorno del problema (3.1), se calculan los valores propios y las funciones propias del problema de Sturm-Liouville: (X00(x)+λX(x) = 0, X(0) = 0,X(1) = 0.(3.2) Se tiene que para k=1,2,... los valores propios son λk= (kπ)2y las funciones propias Xk(x) = sin(kπx). Resolviendo por otro lado T0(t) + λT(t) = 0 como un problema de primer orden se tiene como solución Tk(t) = e−λkt=e−(kπ)2t. Por lo tanto, se sigue que para cada k=1,2,... uk(x,t) = Xk(x)Tk(t) = e−(kπ)2tsin(kπx)es solución de la ecuación en derivadas y satisface las condiciones de contorno. 19 20 Capítulo 3. Simulación numérica Aplicando la condición inicial f(x) = sin(2πx)y suponiendo que la solución es formalmente una combinación lineal infinita de las soluciones anteriores u(x,t), se llega finalmente a: u(x,t) = e−4π2tsin(2πx). Se ha implementado en Matlab el método de Euler explícito para la resolución del problema (3.1). Llamando a la función1: u = ecalorexpl1(Nx, Mt, L, T) se ha representado la solución numérica y la exacta en la Figura 3.1. Se ha elegido como Tfinal =0.1, h=1 6,k=1 80. Es decir, los valores de los parámetros en la función corresponden con Nx =7, Mt =8, T=0.1 y L=1. Figura 3.1: Solución exacta y numérica del problema (3.1) Aunque la solución numérica no da problemas de convergencia y parece aproximase a la solución exacta, es claro que no es la mejor que se puede encontrar. Se puede comprobar que usando una malla más fina, con h=1 20 yk=1 800 la aproximación numérica es mucho mejor. En este caso, los valores de los parámetros son Nx =21, Mt =80, T=0.1 y L=1. Figura 3.2: Solución exacta y numérica del problema (3.1) 1Ver Anexo-Sección 3.1.1.a (página 33). Modelos matemáticos en oncología Simulación numérica 3.1. Solución numérica: método explícito 21 ¿Se puede elegir hykaleatoriamente? Se eligen, por ejemplo, h=1 23 yk=1 800, es decir, Nx =24, Mt =80, T=0.1 y L=1 como valores de los parámetros y se puede observar en la Figura 3.3 la solución obtenida . Figura 3.3: Solución numérica del problema (3.1) que oscila Claramente para estos valores de hykla solución numérica oscila. ¿Por qué sucede esto? El método de diferencias finitas explícito requiere de una condición de estabilidad. Como se comenta en el capítulo 2, hay una restricción en el paso de discretización temporal en función de la espacial. Por último se puede hacer un estudio de los errores cometidos con el método. Se comienza calculando el error cometido en los dos casos estudiados anteriormente, para corroborar con más exactitud que la solución numérica en el segundo caso se aproxima mucho más a la solución exacta. Es decir, se estudian los errores para h=1 6,k=1 80 yh=1 20,k=1 800. Se añade en el programa correspondiente en Matlab, la orden para el cálculo del error2en norma infinito, calculando previamente la solución exacta en cada nodo, y se obtienen los siguientes resultados: NxMtError 7 8 0.0112 21 80 0.0013 donde como era de esperar 0.0013 <0.0112; es decir, el error cometido con los valores tomados de hyken el segundo caso es mucho menor que el error cometido en el primer caso. 3.1.2. Ecuación de reacción-difusión lineal De forma más genérica, se puede hacer un estudio de los errores cometidos para distintos valores de hyk, por ejemplo con la siguiente ecuación de reacción-difusión lineal:        ut= (1+x2)uxx −3u,x∈(−1,1),t>0, u(x,0) = 1+x2,x∈[−1,1], u(−1,t) = u(1,t) = 2e−t,t≥0. (3.3) 2Ver Anexo-Sección 3.1.1.b (página 34). Autora: Marta Gómez Gómez 22 Capítulo 3. Simulación numérica Se parte de Nx =10, Mt =10 y se van duplicando los pasos de malla tanto en espacio como en tiempo, con el fin de verificar el orden de convergencia del método explícito estudiado en el capítulo 2. HHHHH H Nx Mt10 20 40 80 10 6.6407e−06 3.2952e−06 1.6453e−06 8.2205e−07 20 7.0854e−06 3.5222e−06 1.7561e−06 8.7682e−07 40 8.7659e−06 3.5615e−06 1.7741e−06 8.8544e−07 80 62.8877 4.0132e04 23.0589 8.8885e−07 Si se observan los valores de la diagonal en la tabla, se ve como el valor del error cometido va reduciéndose por dos a medida que dividimos por la mitad hyk; lo cual es debido a que el orden de convergencia del método es O(h2+k). 3.1.3. Ecuación de Fisher-Kolmogorov Un segundo objetivo es aplicar el método de diferencias finitas explícito a otras ecuaciones, como por ejemplo a la ecuación de Fisher-Kolmogorov y estudiar el comportamiento de las soluciones para distintos valores de los parámetros de difusión y proliferación. Como se ha visto en el primer capítulo, la ecuación de Fisher-Kolmogorov (FK) se utiliza en la modelización de tumores cerebrales (gliomas), describiendo el comportamiento de una densidad P=P(x,t)de células cancerígenas que pueden migrar y proliferar, y la cual viene definida por: ∂P ∂t=d∂2P ∂x2+ρÇ1−P e PåP,x∈[a,b],t≥0, donde d,ρye Pson constantes que representan el coeficiente de difusión (migración) celular, la tasa de proliferación y la densidad tisular máxima, respectivamente. Notar que las unidades de Pye P son número de células/longitud, dse mide en longitud2/tiempo y ρes inversamente proporcional al tiempo. El intervalo espacial donde se resuelve es [a,b], escogiendo a y b de manera adecuada y recordar que la ecuación de FK se suplementaba normalmente con condiciones de frontera del tipo Neumann. La ecuación de FK se puede simplificar si definimos una nueva variable u(x,t) = P(x,t)/e P, de manera que 0 ≤u(x,t)≤1 para todo x∈[a,b]yt>0. De esta manera, se considera en nuestra simulación la EDP normalizada: ut=duxx +ρ(1−u)u,x∈[a,b],t>0,(3.4) la cual coincide con la estudiada previamente. El intervalo donde se integra es [−L/2,L/2], tomando L=6 cm y se utiliza como condición inicial la siguiente función: u0(x) = U0e−x2 σ2, lo cual supone que las células cancerígenas se distribuyen en el tejido mediante una distribución de tipo gaussiana y donde la constante (amplitud) 0 <U0<1 se toma en el rango de valores U0∈ [10−3,10−1]. Los parámetros de difusión y proliferación se supone que cumplen d∈[1,103]mm2/año Modelos matemáticos en oncología Simulación numérica 3.1. Solución numérica: método explícito 23 yρ∈[10−1,102]año −1. El intervalo de tiempo que se quiere explorar es t∈[0,T]con T∈[10−1,10] años. El parámetro restante σ, se elige de manera que la distribución espacial inicial del tumor esté confinada en una región suficientemente pequeña del intervalo [−L/2,L/2]. Se pretende ver cómo la solución de nuestro problema varía en función de los parámetros dyρ. Se ha implementado en Matlab el método explícito de diferencias finitas para la resolución del problema (3.4) y se llama a la siguiente función3: u = FKi(Nx, Mt, L, T, D, Rho, U0, sigma) En todos los casos se usan los siguientes valores para el resto de argumentos: Nx =10, U0 =0.1, T=1, L=1, sigma =0.1 y Mt el necesario en cada caso para que converja. En primer lugar se supone fijo el parámetro de proliferación, ρ=10 y se compara qué diferencias hay si el parámetro de difusión es muy pequeño, es decir, d=1 o comienza a ser significativo como por ejemplo d=10. Figura 3.4: Solución para un parámetro de difusión muy pequeña en la izquierda y una difusión grande en la derecha Se observa por un lado, como para un parámetro no excesivamente grande de proliferación, independientemente de lo que el tumor se difunda, en un período de un año, ambos tumores llegan a ocupar el máximo del tejido (suponiendo que este sea u=1 con la ecuación normalizada). Por otro lado, se ve como para un mismo valor del parámetro de proliferación, con una difusión significativa (figura de la derecha) el número de células se expande en un período muy breve de tiempo cubriendo toda la línea de integración y sin embargo, con una difusión muy pequeña (figura de la izquierda) no la cubre ni siquiera cuando u=1, es decir cuando alzanca el máximo posible, al cabo de un año. 3Ver Anexo-Sección 3.1.3 (página 34). Autora: Marta Gómez Gómez 24 Capítulo 3. Simulación numérica Se supone ahora fija la difusión d=10, y se compara qué sucede para distintos valores de proliferación, uno apenas inexistente con ρ=0.1 y otro muy elevado como por ejemplo ρ=100. Figura 3.5: Solución para un parámetro de proliferación muy pequeño en la izquierda y una proliferación elevada en la derecha En este caso, lo primero que se observa es que debido al alto valor del parámetro de difusión, ambos tumores cubren toda la línea de integración casi instantáneamente al comienzo del año. En cuanto a la proliferación, en este ejemplo se ve clara la diferencia entre ambos tumores para un mismo parámetro de difusión. Con un valor de proliferación casi inexistente (figura de la derecha), a pesar de que el tumor se difunde muy rápido, apenas cubre la capacidad del tejido, quedándose a lo largo de todo el año en u=0.007 aproximadamente. Sin embargo, con un valor muy alto de proliferación (figura de la izquierda), no sólo se difunde rápidamente, sino que además en apenas un mes y medio el tumor ya ha cubierto el máximo del tejido siendo éste u=1. 3.2. Solución numérica: método IMEX en dos dimensiones 3.2.1. Ecuación de Fisher-Kolmogorov En esta sección, se va a considerar el problema de Fisher-Kolmogorov en dos dimensiones. Para aproximar la solución de dicho problema se va a usar el método IMEX explicado en el capítulo 2. Se ha implementado en Matlab el método IMEX para resolver numéricamente la ecuación: ut=d(uxx +uyy)+ρu(1−u).(3.5) Se llama a la función4: u = Fk2d(Nx, Mt, L, T, D, Rho, U0, sigma) y se fijan todos los argumentos. Se tiene como objetivo ver la evolución de un único tumor a lo largo de un año con un valor de proliferación ρ=10 y de difusión d=1. Los valores de los argumentos que se usan son: Nx =40, Mt =30, D=1, Rho =10, U0 =0.1, sigma =0.1, L=6 y T=1. Se parte de la misma distribución inicial que en el apartado anterior, donde al comienzo del año las células están muy concentradas en un punto, cubriendo aproximadamente el 5% del tejido. 4Ver Anexo-Sección 3.2.1 (página 35). Modelos matemáticos en oncología Simulación numérica 3.2. Solución numérica: método IMEX en dos dimensiones 25 Figura 3.6: Distribución inicial de las células tumorales Se observa la evolución en cuatro períodos de tiempo distintos, de 3 meses de longitud cada uno. Figura 3.7: Evolución del tumor al cabo de 3, 6, 9 y 12 meses Se observa que en los primeros tres meses apenas ha proliferado (cubre un 6% del tejido aproximadamente) y sin embargo ya empieza a difundirse. Al cabo de medio año, se empiezan a notar los efectos de la proliferación (cubre un 20% del tejido aproximadamente) y continua difundiéndose. A los nueve meses, duplica casi su valor cubriendo aproximadamente la mitad del tejido y ya se ha difundido por casi toda la malla. Finalmente, se ve como tras un año está completamente difundido y su valor es aproximadamente u=0.75 siendo 1 el máximo valor que puede alcanzar. Autora: Marta Gómez Gómez Anexo Simulación numérica Sección 3.1.1.a En primer lugar, se ha programado el método de diferencias finitas explícito para la resolución de la ecuación del calor, con el fin de comparar gráficamente la solución exacta con la numérica. function u = ecalorexpl1(Nx, Mt, L, T) hx = L/(Nx-1); %Paso espacial en x ht = T/(Mt-1); %Paso temporal s = ht/hx^2; %Inicializamos la matriz solucion y vectores posicion y tiempo u=zeros(Nx, Mt); x=zeros(1,Nx); t=zeros(1,Mt); for j=1:Nx x(j) = (j-1)*hx; end for m=1:Mt t(m) = (m-1)*ht; end %Imponemos la condicion inicial u(x,0) for j=1:Nx u(j,1) = sin(2*pi*x(j)); end %Condiciones de contorno u(1,1:Mt)=u(1,1); u(Nx,1:Mt)=u(Nx,1); %Formula de recurrencia explicita for m=1:Mt-1 for j=2:Nx-1 u(j,m+1)=u(j,m) + s*(u(j+1,m) - 2*u(j,m) + u(j-1,m)) end end %Solucion exacta 33 34 Capítulo 3. Anexo for i=1:101 Xx(i)=(i-1)/100; Z(i)=exp(-4*pi*pi*0.1)*sin(2*pi*Xx(i)); end %Dibujamos ambas soluciones hold on grid on plot(x,u(:,Mt), ’m--’,’linewidth’, 2) plot(Xx, Z, ’linewidth’, 2) xlabel(’u(x,t)’) end %Final del programa Sección 3.1.1.b Se calcula también el error cometido con el método para los casos estudiados gráficamente, sin más que añadir al programa anterior el cálculo del error y definir la solución exacta en todos los nodos. Xe = zeros(1,Nx); for i=1:Nx Xe(i) = exp(-4*pi*pi*0.1)*sin(2*pi*x(i)); end error = max(abs(Xe’ - u(:,Mt))) Sección 3.1.3 Posteriormente, se modifica ligeramente el programa, para resolver ahora la ecuación de Fisher mediante un método nuevamente explícito. En este caso no se puede calcular la solución exacta, luego no podemos compararlas. Se trata de cambiar los valores de la difusión y la proliferación para ver como cambian distintos tumores a lo largo del tiempo en función de qué características tengan. function u = FKi(Nx, Mt, L, T, D, Rho, U0, sigma) hx = L/(Nx-1); %paso espacial en x ht = T/(Mt-1); %paso temporal s = D*ht/hx^2; %Inicializamos la matriz solucion y vectores posicion y tiempo u=zeros(Nx,Mt); x=zeros(1,Nx); t=zeros(1,Mt); for j=1:Nx x(j) = -(L/2) + (j-1)*hx; end for m=1:Mt t(m) = (m-1)*ht; Modelos matemáticos en oncología Simulación numérica 35 end %Imponemos la condicion inicial u(x,0) for j=1:Nx u(j,1) = U0*exp(-x(j)*x(j)/sigma); end %Formula de recurrencia explicita for m=1:Mt-1 for j=2:Nx-1 u(j,m+1)=u(j,m) + s*(u(j+1,m) - 2*u(j,m) + u(j-1,m))+ ht*Rho*(1-u(j,m))*u(j,m); end u(1,m+1) = u(1,m) + s*(2*u(2,m)-2*u(1,m)) + ht*Rho*(1-u(1,m))*u(1,m); u(Nx,m+1) = u(Nx,m) + s*(2*u(Nx-1,m)-2*u(Nx,m)) + ht*Rho*(1-u(Nx,m))*u(Nx,m); end %Dibujar surf(t,x,u) xlabel(’t’, ’fontname’, ’Times new Roman’, ’fontsize’, 20) %Eje t ylabel(’x’, ’fontname’, ’Times new Roman’, ’fontsize’, 20) %Eje x xlabel(’u(x,t)’, ’fontname’, ’Times new Roman’, ’fontsize’, 20) colormap hsv colorbar end %Final del programa Sección 3.2.1 Por otra parte, se ha programado un método IMEX en dos dimensiones para la resolución de nuevo de la ecuación de Fisher-Kolmogorov. El objetivo en este caso, es observar la evolución de un único tumor a lo largo de un año con todos los valores de los parámetros fijos. function u = Fk2d(Nx, Mt, L, T, D, Rho, U0, sigma) hx = L/(Nx-1); %paso espacial en x ht = T/(Mt-1); %paso temporal s = D*ht/hx^2; %Inicializamos la matriz solucion y vectores posicion y tiempo u=zeros(Nx*Nx,Mt); x=zeros(1,Nx*Nx); t=zeros(1,Mt); x1=-L/2:hx:L/2; y1=-L/2:hx:L/2; [X Y]=meshgrid(x1,y1); for j=1:Nx for i=1:Nx x(i,j)=-L/2+(i-1)*hx; Autora: Marta Gómez Gómez 36 Capítulo 3. Anexo y(i,j)=-L/2+(j-1)*hx; end end for m=1:Mt t(m) = (m-1)*ht; end %Imponemos la condicion inicial u(x,0) for j=1:Nx for i=1:Nx ind=(j-1)*Nx+i; u(ind,1) = U0*exp(-(x(i,j)*x(i,j)+y(i,j)*y(i,j))/sigma); end end %Definicion de los nodos en la malla (suponemos una malla cuadrada) A=sparse(Nx*Nx, Nx*Nx); i = 1; A(1,1) = 4*s + 1; A(1,2) = -2*s; A(1,1+Nx) = -2*s; for i=2:Nx-1 A(i,i) = 4*s + 1; A(i,i+1) = -s; A(i, i-1) = -s; A(i, i+Nx) = -2*s; end i = Nx; A(Nx,Nx) = 4*s + 1; A(Nx, Nx-1) = -2*s; A(Nx, 2*Nx) = -2*s; for j = 2:(Nx-1) i=i+1; A(i,i) = 4*s + 1; A(i,i+1) = -2*s; A(i,i+Nx) = -s; A(i,i-Nx) = -s; for k = 2:(Nx-1) i=i+1; A(i,i) = 4*s + 1; A(i,i+1) = -s; A(i,i-1) = -s; A(i,i+Nx) = -s; A(i,i-Nx) = -s; end Modelos matemáticos en oncología Simulación numérica 37 i=i+1; A(i,i)= 4*s + 1; A(i,i-1)= -2*s; A(i,i+Nx)= -s; A(i,i-Nx)= -s; end i=i+1; A(i,i) = 4*s + 1; A(i,i+1) = -2*s; A(i,i-Nx) = -2*s; for k = 2:Nx-1 i=i+1; A(i,i) = 4*s + 1; A(i,i+1) = -s; A(i, i-1) = -s; A(i, i-Nx) = -2*s; end i=i+1; A(i,i) = 4*s + 1; A(i,i-1) = -2*s; A(i,i-Nx) = -2*s; %Resolucion del sistema for m=1:Mt-1 for j=1:Nx*Nx b(j) = u(j,m)+ ht*Rho*u(j,m)*(1-u(j,m)); end u(:,m+1) = A\b’ end %Dibujar for i=1:Mt um=reshape(u(:,i),Nx,Nx) surf(X,Y,um) end end %Final del programa Radioterapia y casos prácticos Sección 3.3 Por último, se ha programado la ecuación que nos da el crecimiento del número de células cancerosas de un tejido en un período de tiempo cuando un paciente se somete a radioterapia y posteriormente. El objetivo es comparar en función de las características de un tumor, cual es la forma más óptima de distribuirle al paciente las dosis de radioterapia. Autora: Marta Gómez Gómez 38 Capítulo 3. Anexo %DATOS INICIALES NECESARIOS PARA EL PROGRAMA %En primer lugar Sf nos indica la fraccion de supervivencia de las celulas %tumorales tras cada sesion de radioterapia, pudiendo valer 0.48 o 0.85 en %funcion de como de sensible sea el tumor a las radiaciones. Sf = 0.85; %Sf = 0.48; %Por otro lado, definimos rho, el parametro de proliferacion de las celulas %tumorales, el cual hemos deducido que viene dado por rho = log(2) / Tdupl, %pudiendo ser el tiempo de duplicacion 20 o 80, en funcion de como de %agresivo es el tumor. %rho = log(2)/20; rho = log(2)/80; %En cuanto al numero de celulas tumores, conocemos el numero de celulas %tumorales en el instante inicial y el numero maximo de celulas tumorales %que soporta nuestro tejido. N(1) = 1e9; Nmax = 3e9; %Realizaremos la simulacion durante 1000 dias y con un reparto de las 30 dosis %de radioterapia continuo o semanalmente Tend = 1000; t = [7:37]; %Continuo %t = [0 7:7:210]; % Cada semana %PROGRAMA PRINCIPAL length(t) for j=1:length(t)-1 N(j+1) = Sf*N(j)*exp(rho*(t(j+1)-t(j)))/(1+(N(j)/Nmax)*(exp(rho*(t(j+1)-t(j)))-1)); Tn = t(j+1); Nn = N(j+1); end; time = [t(end):0.1:Tend]; Number = Nn*exp(rho*(time-Tn))./(1+(Nn/Nmax)*(exp(rho*(time-Tn))-1)); DL = 0.7*Nmax; size(N) size(t) N %REPRESENTACION GRAFICA DEL VALOR DE N A LO LARGO DE LOS 1000 DIAS plot(t,log(N),’o’,time,log(Number),’-’,[0 Tend],[log(DL) log(DL)],’--r’); N(end) Modelos matemáticos en oncología Simulación numérica Bibliografía [GG] RA. Gatenby,ET. Gawlinski,“A reaction-difussion model of cancer invasion”, Cancer Research 56 (1996), 5745-5753. [HV] W. Hundsdorfer,J.G. Verwer,Numerical Solution of Time-Dependent Advection-DiffusionReaction Equations, Springer, 2003. [JV] M. Joiner,A. Van der Kogel,Basic Clinical Radiobiology, Fourth Edition, Hodder Arnold, 2009. [L] R. Leveque,Finite difference methods for differential equations, SIAM, Philadephia, 2007. [M] J. Murray,Mathematical Biology I: An Introduction, Third Edition, Springer, 2002. [M] J. Murray,Mathematical Biology II: Spatial Models and Biomedical applications, Third Edition, Springer, 2003. [SRA] KR. Swanson,RC. Rostomily,EC. Alvord Jr,“A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle”, Brithish Journal of Cancer 98 (2008), 113-119 . [T] J.W. Thomas,Numerical Partial Differential Equations, Finite Difference Methods, Texts in Applied Mathematics 22, Springer, 1995. [TW] A. Tveito,R. Winther,Introduction to partial differential equations, A computational approach, Texts in Applied Mathematics 29, Springer, 1998. 39