Modelos da bioloxía matemática

Traballo Fin de Grao Modelos da bioloxía matemática Adrián García Muiño 2020/2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Modelos da bioloxía matemática Adrián García Muiño Xullo 2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Traballo proposto Área de Coñecemento: Matemática aplicada. Título: Modelos da bioloxía matemática. Breve descrición do contido Proporanse o estudo e a resolución numérica, e no seu caso, tamén teórica, de algúns modelos da bioloxía matemática baseados en ecuacións diferenciais. Recomendacións Dominio dunha linguaxe de programación como Matlab. Outras observacións iii Índice xeral Resumo vii Introdución ix 1. Dúas poboacións en competencia. Modelo do cancro 1 1.1. Modelo matemático do cancro . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2. Simulacións numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2. Interacción cancro-sistema inmunolóxico 13 2.1. Modelo matemático do cancro . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2. Simulacións numéricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3. Terapias do cancro 25 3.1. Viroterapia .................................... 25 3.2. Tratamento con GM-CSF ............................ 28 3.3. Simulación numérica do modelo de terapia do cancro . . . . . . . . . . . . . 29 Conclusións 37 Apéndice Ecuacións diferenciais 38 Bibliografía 51 v Resumo A bioloxía matemática é unha área de coñecemento que estuda os principios que rexen a estrutura, o desenvolvemento e o comportamento dos sistemas biolóxicos, e que ten como obxectivo a representación matemática, tratamento e modelado destes sistemas. Veremos como xurdiu, e cales foron os descubrimentos máis importantes nos primeiros pasos desta área. Ademais, tamén explicaremos o seu compoñente principal, o proceso de modelaxe. Nos seguintes capítulos centrarémonos en distintos modelos sobre o cancro, unha enfermidade caracterizada pola transformación das células, de modo que proliferan de maneira anormal e incontrolada, e as súas terapias, para reflectir como as matemáticas poden axudar no campo da oncoloxía. Mostraremos ademais os códigos Matlab empregados para resolver e interpretar os distintos modelos. Abstract Mathematical biology is an area of knowledge that studies the principles that govern the structure, development and behaviour of biological systems, and whose objective is the mathematical representation, treatment and modelling of these systems. We will see how it came about, and what were the most important discoveries in the first steps of this area. In addition, we will also explain its main component, the modelling process. In the following chapters we will focus on different models of cancer, a disease characterized by the transformation of cells so that they proliferate in an abnormal and uncontrolled way, and its therapies, to reflect how mathematics can help in the field of oncology. We will also show the Matlab codes used to solve and interpret the different models. vii xiv INTRODUCIÓN E, por último, a bioloxía organizacional, que trata sobre a xerarquía de estruturas e sistemas biolóxicos complexos que definen a vida mediante unha aproximación reduccionista, onde a natureza humana está determinada polos xenes, e as propiedades dos individuos e as súas accións son consecuencia inevitábel destes. Dentro de todas estas áres de investigación da bioloxía matemática, un esperanzador desafío está vinculado a entender as dinámicas do cancro dende as perspectivas morfolóxica, xenómica, proteómica e matemática, e trasladar os modelos e datos á práctica clínica [24]. Neste traballo centrarémonos nos distintos modelos sobre o cancro e as terapias para este, para reflectir como a matemática pode axudar no campo da bioloxía, e máis concretamente nunha enfermidade que actualmente é unha das principais causas de morte no mundo. O cancro, según o dicionario, é unha enfermidade que se caracteriza pola transformación anormal das células, de modo que proliferan de maneira anormal e incontrolada. Ao contrario que as células normais, as células canceríxenas non morren logo dun período de tempo programado, e divídense case sen límite. Esta multiplicación forma unhas masas, chamadas tumores, que poden destruír e substituír as células normais. Algúns cancros poden non formar tumores, e algúns tumores poden ser malignos, é dicir, producen metástase. E hai outros que crecen a un ritmo lento e que non se infiltran nos tecidos veciños, que son os chamados tumores benignos [25]. Aínda que xa existían algúns modelos matemáticos para o cancro anteriormente, como o publicado por Peter Armitage (1924) e Richard Doll (1912-2005) en 1954 [1], que trataba a primeira teoría multiestadía do cancro (según a cal se deberían superar entre seis e sete cambios xenéticos nas células antes de adquirir un fenotipo completamente canceríxeno), foi nos anos 70 cando se comezaron a desarrollar máis modelos sobre esta enfermidade e as súas terapias. Principalmente, os modelos de ecuacións diferenciais ordinarias e en derivadas parciais son as ferramentas máis empregadas no estudo do crecemento de tumores e na forma na que se difunden sobre os tecidos que os rodean. A pesar do gran número de traballos sobre este campo que se poden atopar actualmente, aínda non existe un modelo que proporcione unha predición e caracterización do comportamento para o crecemento de tumores canceríxenos nas súas múltiples formas e para calquera tipo de poboación, nin tampouco existen modelos para algún tipo de terapia que consiga eliminar o tumor en calquera tipo de situación [11]. Así, neste traballo trataremos de mostrar algúns dos modelos que existen, tratando tamén de ensinar as principais vías de estudo que se están a seguir actualmente neste campo. No primeiro capítulo trataremos de amosar un primeiro modelo xeral baseado na dinámica de dúas poboacións en competencia, que no noso caso serán as células sanas e as canceríxenas. Este tipo de modelo é un dos máis coñecidos, sobre todo no contexto da eco- INTRODUCIÓN xv loxía, e trata de ver se dúas poboacións serán capaces de sobrevivir xuntas, ou se algunha delas acabará coa outra [4]. No segundo capítulo, trataremos un modelo máis específico, onde amosaremos a interacción entre o cancro e o sistema inmunolóxico. Veremos os distintos tipos de células involucradas no proceso, e como actúa cada unha delas. Mostraremos, empregando códigos en Matlab, as gráficas das solucións dos modelos estudados. Por último, no terceiro capítulo, veremos dous modelos de terapias para o cancro, a viroterapia e o tratamento co fármaco GM-CSF. Á diferencia dos outros dous capítulos, nos que tratamos con ecuacións diferenciais, neste caso empregaremos ecuacións en derivadas parciais. Análogamente ao anterior capítulo, volveremos a mostrar as gráficas das solucións dos modelos empregando códigos en Matlab. Finalmente, para facilitar a lectura do traballo, incluimos un apéndice cun resumo da teoría clásica de sistemas de ecuacións diferenciais ordinarias e en derivadas parciais. Capítulo 1 Dúas poboacións en competencia. Modelo do cancro Nesta sección, seguiremos o capítulo 6 do libro Introduction to mathematical biology [4]. Comenzaremos estudando problemas dos chamados depredador-presa, onde dúas poboacións compiten por ter o máximo número de individuos. O fin é obter un modelo sobre o cancro pensando que as células normais e canceríxenas son dúas poboacións en competencia deste tipo. Neste sentido a competencia é unha interacción entre dous organismos que comparten os mesmos recursos para poder sobrevivir. Cando se coñecen suficientes datos sobre a competencia entre dous organismos, podemos axudarnos das matemáticas para predicir si tan só un ou ambos conseguirán sobrevivir a longo prazo. Para poder chegar a un modelo depredador-presa suficientemente completo, imos comezar por ver os primeiros modelos deste tipo, para a continuación, despois dalgunhas modificacións, chegar a modelos óptimos de dúas poboacións en competencia. Como xa mencionamos anteriormente, no século XVIII, Malthus publicou o seu ensaio sobre a dinámica de poboacións baseándose nun crecemento exponencial [13]. Alí, considerando xo tamaño da poboación no instante t, e kunha constante de proporcionalidade, él propuxo a seguinte ecuación: dx dt =k x, é dicir, supoñemos que a tasa de aumento da poboación é proporcional á poboación nese instante. De maneira elemental, podemos probar que a solución desta ecuación é x=Cek t, onde Cé unha constante arbitraria, é dicir, unha función exponencial, de aí que falemos do modelo exponencial de crecemento de Malthus. No ano 1838, Verhulst, depois de ler o ensaio de Malthus decatouse de que a miúdo 1 2CAPÍTULO 1. DÚAS POBOACIÓNS EN COMPETENCIA. MODELO DO CANCRO non é realista que a poboación creza indefinidamente, pois hai casos onde o número de individuos se estanca. Entón propuxo un novo modelo de crecemento loxístico [20], onde considerou ra tasa de crecemento da poboación e ka capacidade de carga do entorno. Podemos escribir a ecuación que goberna o modelo da forma seguinte: dx dt =rx(1 −x k), é dicir, o crecemento do modelo de Malthus amortíguase cun factor 1−x/k, que na práctica debe ser unha cantidade entre 0e1. Despois diso, como ben mencionamos tamén anteriormente, Lotka e Volterra desarrollaron de maneira independente un sistema de ecuacións que describen a relación entre dúas especies que comparten os mesmos recursos [7][21]. Partindo ambos da ecuación anterior, e baseándonse no que ocorre nunha reacción química, tanto Lotka como Volterra, axudado este último por un estudo do seu cuñado, o biólogo italiano Umberto D’Ancona (1896- 1964), tiveron a idea de aplicarlle a esta ecuación a lei química da acción de masas (que establece a relación entre as masas de reactivos e productos nun equilibrio químico a unha temperatura determinada). Este modelo, que surxe de xuntar ambos, e que actualmente se coñece co nome de Lotka-Volterra, é o primeiro de moitos modelos de interacción entre especies, e defínese mediante un sistema formado polas seguintes ecuacións diferenciais ordinarias: dx dt =r1x−a1xy, dy dt =a2xy −r2y, sendo a primeira a ecuación que modela a poboación das presas, e a segunda a de depredadores. Neste caso, xserá o tamaño da poboación de presas, yo de depredadores, r1er2as tasas de crecemento das presas e dos depredadores, respectivamente, a1o éxito na caza do depredador, que afecta á presa, e a2o éxito na caza que afecta ó depredador. Na primeira ecuación, as presas teñen un crecemento exponencial, amortiguado pola interacción da caza dos depredadores sobre estas. Na segunda, o crecemento do número de depredadores depende en primeiro lugar da interacción coas presas, e en segundo lugar dun decrecemento exponencial. Inspirándonos nos modelos vistos anteriormente, e noutros que viñeron despois, podemos centrarnos no modelo que nos ocupa nesta sección, que veremos primeiro aplicado a dúas poboacións xerais, e máis adiante adaptaremos ao caso do cancro, dx dt =r1x(1 −x k1 )−b1xy, (1.1) dy dt =r2y(1 −y k2 )−b2xy. (1.2) 3 Neste caso, en (1.1), r1é a tasa de crecemento da especie x,k1é a capacidade de carga que limita o seu crecemento, e b1é a tasa á que o organismo ymata ó organismo x. De maneira similar sucede na ecuación (1.2). Temos pois dúas poboacións que compiten polo mesmo espazo, e na que ambas poden ser presas ou depredadores. En ambas ecuacións, o primeiro termo da dereita reflicte un crecemento loxístico de cada especie, amortiguado por un segundo termo que reflexa a interacción entre ambas. Partindo do modelo anterior, imos ver algunhas das súas propiedades, para poder aplicalas despois ao modelo do cancro. A primeira que imos estudar será os puntos de equilibrio e o seu tipo. Para iso temos que ver en que puntos se anulan simultaneamente as ecuacións (1.1) e (1.2). Para poder estudar todos os puntos de equilibrio do sistema temos que resolver o seguinte sistema de ecuacións alxébricas: r1x(1 −x k1 )−b1xy = 0, r2y(1 −y k2 )−b2xy = 0. Sacando factor común xda primeira ecuación, e factor común yda segunda ecuación, cun razoamento elemental chegamos aos puntos de equilibrio (0,0),(k1,0) e(0, k2). O cuarto punto de equilibrio é a solución do seguinte sistema de ecuacións lineais: r1 k1 x+b1y=r1,(1.3) b2x+r2 k2 y=r2.(1.4) Empregamos a regra de Cramer para encontrar a súa solución de maneira directa. O determinante do sistema é: D= r1 k1b1 b2r2 k2 =r1r2 k1k2 −b1b2=r1r2 k1k2 (1 −β1β2) onde para abreviar poñemos β1=k1b1/r1,β2=k2b2/r2. Finalmente obtemos que a primeira incógnita vén dada por x=1 D r1b1 r2r2 k2 =1 D(r1r2 k2 −r2b1) = 1 1−β1β2 k1k2 r1r2 (r1r2 k2 −r2b1) = k1−k2β1 1−β1β2 , e análogamente para a segunda incógnita temos que: y=k2−k1β2 1−β1β2 . 4CAPÍTULO 1. DÚAS POBOACIÓNS EN COMPETENCIA. MODELO DO CANCRO Imos ver agora de que tipo é cada un dos puntos, sabendo que (k1,0) significa que a segunda poboación se extingue, e (0, k1)correspóndese con unha situación na que a primeira poboación queda extinguida. Empregando a matriz Xacobiana (ver apéndice), podemos observar que o punto (0,0) é inestable, (k1,0) é estable se k1> r2/b2, o que quere dicir que a capacidade de carga do ecosistema para o organismo xten que ser maior que a tasa que proporciona o cociente entre o crecemento da especie ye a tasa á que o organismo xmata ó organismo y, e (0, k2)é estable se k2> r1/b1, é dicir, a mesma situación que para o punto anterior, pero intercambiando as especies. Por último temos o punto estable no que ambas poboacións poden coexistir, que vén dado pola solución do sistema (1.3)-(1.4), é dicir (k1−k2β1 1−β1β2 ,k2−k1β2 1−β1β2 ),onde βi=kibi ri (i= 1,2) Este punto estable é de relevancia biolóxica tan só se as dúas compoñentes son positivas. Iso ocorre tan só se: k1−k2β1>0, k2−k1β2>0e1−β1β2>0,(1.5) ou, k1−k2β1<0, k2−k1β2<0e1−β1β2<0.(1.6) Agora ben, as primeiras condicións (1.5) cúmprense cando k1−k2β1>0⇐⇒ k1−k2 k1b1 r1 >0⇐⇒ k1(1 −k2b1 r1 )>0⇐⇒ 1−k2b1 r1 >0⇐⇒ k2>r1 b1 k2−k1β2>0⇐⇒ k2−k1 k2b2 r2 >0⇐⇒ k2(1 −k1b2 r2 )>0⇐⇒ 1−k1b2 r2 >0⇐⇒ k1>r2 b2 . Un cálculo elemental mostra que en (1.5) a terceira condición verifícase sen hipótesis adicionais. Traballando de maneira análoga en (1.6), obtemos desigualdades análogas pero co signo cambiado. En resumo, chegamos a que as dúas compoñentes serán positivas se k1>r2 b2 e k2>r1 b1 ou k1<r2 b2 e k2<r1 b1 . Así, como conclusión do modelo xeral, podemos dicir que ambas especies coexistirán sempre que a tasa de matanza bi, sexa maior que ri/kj, a tasa de crecemento dividida pola capacidade de carga, tanto para j= 1, i = 2 como para j= 2, i = 1; ou ben, que sexa menor en ambos casos. 1.1. MODELO MATEMÁTICO DO CANCRO 5 1.1. Modelo matemático do cancro Partindo do modelo anterior, imos a centrarnos agora en describir e buscar as propiedades dun modelo para o cancro baseándonos nos modelos de dúas poboacións en competencia. Lembremos que o modelo onde o crecemento loxístico xeneralizado dunha poboación cunha densidade x, cun termo de amortiguación é dx dt =rx(1 −x K)−µx onde ré a tasa de crecemento, µé a tasa de morte e Ké a capacidade de carga media determinada polos recursos dispoñibles para a poboación. O termo −µx non pode amortiguar totalmente o termo anterior, porque se fora µ>r, entón dx dt +µx −rx =−r Kx2≤0sempre. Calcularemos agora a solución desta ecuación seguindo os pasos para resolver unha ecuación diferencial de Bernoulli, por ser da forma dx dt +P(t)x=Q(t)xα. onde P(t) = µ−r,Q(t) = −r Keα= 2. Comezamos dividindo a ecuación diferencial entre x2e chegamos a 1 x2 dx dt + (µ−r)1 x=−r K. Facendo agora o cambio de variable z=x−1, entón dz dt =−1 x2 dx dt . Substituíndo: −dz dt + (µ−r)z=−r K, e cambiando de signo dz dt −(µ−r)z=r K, onde tanto µ−rcomo r/K son constantes. Esta é unha ecuación diferencial lineal, cuxa ecuación homoxénea correspondente é dz dt −(µ−r)z= 0, e ten por solución xeral z=Ce(µ−r)t, onde Cé unha constante arbitraria. Por outro lado, unha solución particular da ecuación inicial non homoxénea é: zpart =−r K(µ−r) 6CAPÍTULO 1. DÚAS POBOACIÓNS EN COMPETENCIA. MODELO DO CANCRO e entón, a solución xeral da ecuación lineal inicial é: z=Ce(µ−r)t−r K(µ−r). Entón, a solución da ecuación de Bernoulli ven dada por: x=1 Ce(µ−r)t−r K(µ−r) . Así, se µ>r, o denominador é crecente, o que implica que xdecrece co tempo. Polo tanto, para que xcreza co paso do tempo, temos que ter µ>r. O caso onde µ=rnon pode darse porque neste caso un denominador se anula. Se dúas poboacións xeycoexisten no mesmo medio e seguen un crecemento loxístico xeneralizado, entón o modelo podería ser dx dt =r1x(1 −x+y K)−µ1x, dy dt =r2y(1 −x+y K)−µ2y, onde r1er2son as tasas de crecemento das poboacións xeyrespectivamente, e µ1eµ2son as súas respectivas tasas de mortalidade. Supoñemos que as dúas poboacións comparten o mesmo medio, polo tanto, o termo (x+y)/K representa a carga total da poboación x+y no medio de capacidade K. Aplicaremos este modelo ao cancro no tecido humano, onde x representa a densidade de células sanas e yrepresenta a densidade de células canceríxenas no mesmo tecido, coas dúas poboacións competindo polo mesmo espazo. Dado que as células canceríxenas proliferan máis rápido que as normais, tomamos r2> r1. Por simplicidade, asumimos que µ1=µ2=µe tomamos r1> µ para que as dúas poboacións persistan. Asumimos tamén que as condicións iniciais x(0) ey(0) son non negativas. Entón o modelo queda: dx dt =x[r1(1 −x+y K)−µ],(1.7) dy dt =y[r2(1 −x+y K)−µ],(1.8) de onde observamos que non pode existir un punto de equilibrio (¯x, ¯y)con ¯x > 0,¯y > 0. En efecto, se tal punto existira, entón (¯x+ ¯y)/k sería igual a 1−µ/r1e1−µ/r2, o cal é imposible. Agora ben, de maneira inmediata vemos que o punto (0,0) é un punto de equilibrio con Xacobiano 1.1. MODELO MATEMÁTICO DO CANCRO 7 J(0,0) = r1−µ0 0r2−µ!. Dado que ambos valores propios son positivos, (0,0) é un nodo inestable (ver apéndice). Queda entón por considerar os puntos de equilibrio ((1 −µ r1 )K, 0) e(0,(1 −µ r2 )K), que son tamén solucións do sistema de ecuacións seguinte x[r1(1 −x+y K)−µ]=0, y[r2(1 −x+y K)−µ] = 0. Calculando o Xacobiano en cada un deles, atopámonos con que o primeiro é inestable, e o segundo estable. Iso significa que o estado no que o corpo está libre de células canceríxenas é inestable, mentras que o estado onde todas as células son canceríxenas é estable. A gráfica (1.1) mostra o retrato de fases [6] do modelo do cancro (1.7)-(1.8). Figura 1.1: Retrato de fases para o modelo do cancro (1.7)-(1.8) O lugar xeométrico Γ1dos puntos tales que o segundo membro da ecuación (1.7) é nulo ven dado por Γ1: 1 −x+y K=µ r1 14 CAPÍTULO 2. INTERACCIÓN CANCRO-SISTEMA INMUNOLÓXICO cancro maligno. E é este o estado de cancro que provoca a inmensa maioría de mortes. Algúns exemplos dos máis de cen tipos diferentes de cancro según o seu tecido de orixe son: os carcinomas, que comezan nas células epiteliais, as células que recubren as superficies do corpo; os sarcomas, que surxen dos tecidos conectivos (fascias), os músculos e a vasculatura (conxunto de vasos sanguíneos do corpo); leucemias, que son o cancro do sistema hematopoiético (células sanguíneas), linfomas, do sistema inmunolóxico, e os retinoblastomas, cancros dos ollos. Imos profundizar agora un pouco máis nos procesos que sofren os distintos tipos de células involucradas no proceso para poder crear un modelo adecuado. A transformación de células normais en células canceríxenas é debida a mutacións en xenes que regulan o crecemento e os cambios na estrutura ou na función que conducen á súa especialización. No contexto do cancro, estes xenes clasifícanse en oncóxenos, xenes que promoven o crecemento e a reprodución, e os xenes supresores de tumores, que son aqueles que inhiben a división celular e a súa supervivencia. A transformación maligna ocorre cando os oncóxenos se sobreexpresan anormalmente, ou cando os xenes supresores de tumores se volven infraexpresados ou discapacitados. Polo xeral, a transformación dunha célula normal nunha célula canceríxena ocorre despois de non unha, senón de varias mutacións xenéticas. Comúnmente, crese que a maioría destas mutacións que conducen ao cancro poden ser debidas a diversas condicións como fumar, factores diabéticos, contaminantes ambientais, exposición a radiación e certas infeccións. Tamén existen factores que son hereditarios. A malignidade do tumor induce, típicamente, unha resposta inmune celular moderada, pero as células canceríxenas tratan de evadir esta resposta inducindo cambios favorábeis no fenotipo das células inmunes. A interacción entre as células canceríxenas e o sistema inmunolóxico é complexa, e afecta á eficacia dos fármacos quimioterápicos. Para determinar a eficacia dos medicamentos contra o cancro, necesitamos desarrollar un modelo matemático de interacción entre as células canceríxenos e o sistema inmune e logo empregalo para evaluar dita eficacia, e será este o obxectivo deste capítulo. 2.1. Modelo matemático do cancro Imos comezar definindo os distintos tipos de células que interveñen na interacción entre células canceríxenas e non canceríxenas. O primeiro tipo do que falaremos son os macrófagos, onde distinguimos entre dous fenotipos: os macrófagos proinflamatorios M1, que producen unha citocina inflamatoria, chamada interleucina IL-12, e os macrófagos antiinflamatorios M2, que producen unha citocina antiinflamatoria, denominada interleucina IL-10. A citocina IL-12 activa as células T, encargadas de destruír as células canceríxenas, 2.1. MODELO MATEMÁTICO DO CANCRO 15 mentras que a citocina IL-10 ten o efecto contrario, inhibe a súa activación. Para evadir o sistema inmunolóxico, as células canceríxenas producen un factor de crecemento transformante β(chamado TGF-β, polas súas siglas en inglés, Transforming Growth Factor) que se adhiren á membrana dos macrófagos M1e comezan un proceso de transformación que cambia o seu fenotipo a macrófagos M2, o que reduce a destrucción de células canceríxenas por parte das células T. O seguinte diagrama resume este proceso, onde a punta de frecha quere dicir produción ou activación, e o acabado en raia significa inhibición ou morte. Figura 2.1: Interacción sistema inmunolóxico - cancro Denotamos por Cá cantidade de células caceríxenas, e por M1,M2,Tá cantidade de macrófagos M1,M2e células T, respectivamente. Baseándonos no proceso descrito anteriormente, podemos escribir as seguintes ecuacións diferenciais para as células dC dt =λCC(1 −C C0 )−µCTC, (2.1) dM1 dt =κ1−˜γTβ ˜ K1+Tβ M1−µM1,(2.2) dM2 dt = ˜γTβ ˜ K1+Tβ M1−µM2,(2.3) dT dt = ˜κT I12 ˜ K2+I10 −µTT, (2.4) 16 CAPÍTULO 2. INTERACCIÓN CANCRO-SISTEMA INMUNOLÓXICO e, para as citocinas dI12 dt =λ12M1−µ12I12,(2.5) dI10 dt =λ10M2−µ10I10,(2.6) dTβ dt =λβC−µβTβ.(2.7) con coeficientes constantes λC, µC, κ1,˜γ, µ, ˜ K1,˜κT,˜ K2,µT,λ12,µ12,λ10,µ10,λβeµβque asumimos positivos para manter o carácter crecente ou decrecente, según conveña. Estas ecuacións reflicten os procesos polos que pasan as células das que falamos anteriormente. Na ecuación (2.1) asumimos un crecemento loxístico xeneralizado das células canceríxenas, C, cun factor de crecemento λCe a eliminación das células do cancro polas células Ta unha tasa µC. Na ecuación (2.2) asumimos unha tasa de produción constante κ1e tasa de mortalidade µdo macrófago M1.Tβcambia o fenotipo de M1aM2, e isto explica o termo ˜γTβ ˜ K1+Tβ M1. O primeiro termo da ecuación (2.3) reflicte, como acabamos de explicar, o cambio de fenotipo de M1aM2, positivo neste caso por tratarse da ecuación de M2, e supoñemos que a tasa de morte µpara os macrófagos M2é a mesma que para os macrófagos M1. Na ecuación (2.4) o primeiro termo representa a activación das células Tpor I12, un proceso inhibido por I10 que aparece no factor 1 ˜ K2+I10 , e o segundo termo explica a morte das células Ta unha tasa µT. Nas ecuacións (2.5)-(2.7) o primeiro termo do lado dereito é un termo de produción das células correspondentes, e o segundo termo explica a degradación. Simplificamos o modelo (2.5)-(2.7) indicando que a dinámina das citocinas é moito máis rápida que a dinámica das células canceríxenas. Polo tanto, podemos supoñer estados estacionarios nas ecuacións (2.5)-(2.7), é dicir, dI12 dt =dI10 dt =dTβ dt = 0 e polo tanto, I12 =cte ×M1, I10 =cte ×M2e Tβ=cte ×C. Substituindo estas relacións nas ecuacións (2.2)-(2.4), o sistema (2.2)-(2.7) redúcese ó sistema de catro ecuacións dC dt =λCC(1 −C C0 )−µCTC, (2.8) dM1 dt =κ1−γC K1+CM1−µM1,(2.9) dM2 dt =γC K1+CM1−µM2,(2.10) dT dt =κT M1 K2+M2 −µTT, (2.11) con coeficientes constantes γ, κTeK1, K2, que de novo asumimos positivos pola mesma razón. 2.1. MODELO MATEMÁTICO DO CANCRO 17 Agora ben, un fármaco quimioterapéutico común é o inhibidor TGF-β. O efecto deste fármaco é aumentar µβ, que aparece na ecuación (2.7), e, polo tanto, diminuir γ. Para determinar si o inhibidor TGF-βpode erradicar o cancro no modelo, imos obter a continuación unhas estimacións sobre a súa solución. Imos introducir o seguinte teorema, cuxa demostración pódese ver en [4], necesario para demostrar as devanditas estimacións. Teorema 2.1. Si unha función x(t)satisfai a inecuación diferencial: dx dt +µx ≤bpara t > 0, ou, dx dt +µx ≥bpara t > 0, onde µebson constantes. Entón, existe un  > 0suficientemente pequeno tal que: x(t)<b µ+, ou, respectivamente, x(t)>b µ+. Do teorema anterior, deducimos as estimacións seguintes: Estimación 1. Para calquera  > 0suficientemente pequeno, cúmprese a seguinte desigualdade: M1(t)>κ1 µ+γ−sendo tsuficientemente grande, (2.12) poñamos, t>T1. Proba.- De feito, como C K1+C<1,(2.13) deducimos da ecuación (2.9) dM1 dt > κ1−(γ+µ)M1. De onde, aplicando o teorema (2.1) se sigue a desigualdade (2.12). Estimación 2. Para calquera  > 0suficientemente pequeno, M2(t)<γκ1 µ2+sendo tsuficientemente grande, (2.14) poñamos, t>T2. 18 CAPÍTULO 2. INTERACCIÓN CANCRO-SISTEMA INMUNOLÓXICO Proba.- De feito, aplicando (2.13) á ecuación (2.10) obtemos dM2 dt < γM1−µM2. Aplicando de novo (2.13), neste caso á ecuación (2.9) temos dM1 dt < κ1−µM1, o que leva, aplicando o teorema (2.1), a M1(t)<κ1 µ+1. Substituíndo isto na desigualdade diferencial de M2, obtemos dM2 dt <γκ1 µ+γ1−µM2, de onde, aplicando de novo o teorema (2.1) obtemos a desigualdade (2.14). Estimación 3. Para calquera  > 0suficientemente pequeno, cúmprese a seguinte desigualdade T(t)>γκ1 α(γ)−para todo tsuficientemente grande, (2.15) poñamos, t>T3, onde α(γ) = 1 κ1 µT(µ+γ)(K2+γκ1 µ2).(2.16) Proba.- Dividindo a desigualdade (2.12) entre a desigualdade (2.14) previa suma da constante K2aos dous membros, obtemos que M1 K2+M2 >(κ1 µ+γ−)1 K2+γκ1 µ2+(2.17) si t>T1+T2. Denotemos A=K2+γκ1 µ2, sexa f(x) = 1 x+A, entón f0(x) = −1 (x+A)2, f00(x) = 2 (x+A)3, evaluando en cero f(0) = 1 A, f0(0) = −1 A2, f00(xi) = 2 (xi+A)3, xientre 0ex. 2.1. MODELO MATEMÁTICO DO CANCRO 19 Pola fórmula de Taylor f(x) = f(0) + xf0(0) + x2 2f00(xi), xientre 0ex, é dicir, 1 x+A=1 A−x A2+x2 (xi+A)3, xientre 0ex, e facendo x= > 0, entón xi>0tamén. Así 2 (xi+A)3>0, e deducimos que 1 A+>1 A− A2>1 A− A− A2=1 A−(1 A+1 A2)=1 A−C1, onde C1=1 A+1 A2, tomando suficientemente pequeno como para que C1 < 1.Usando isto en (2.17), obtemos M1 K2+M2 >κ1 µ+γ−1 A−C1=κ1 µ+γ 1 A−κ1C1 µ+γ+1 A+C12> >κ1 µ+γ 1 A−κ1C1 µ+γ+1 A > κ1 µ+γ 1 A−C2 onde C2=κ1C1 µ+γ+1 A. Si substituímos isto na ecuación (2.11), obtemos a desigualdade dT dt >[κT κ1 µ+γ 1 A−κTC2]−µTT, para t > T1+T2, onde 1 µT κ1 µ+γ 1 A=1 α(γ), pola definición de α(γ)(2.16). Deducimos agora, aplicando o teorema (2.1) que T(t)>κT α(γ)−1 µT κTC2− si t−(T1+T2)é suficientemente grande, poñamos si t>T3. A continuación, observamos que dado que pode tomarse como calquera número positivo pequeno, tamén o poderíamos tomar de maneira que 1 µT κTC2+ 20 CAPÍTULO 2. INTERACCIÓN CANCRO-SISTEMA INMUNOLÓXICO sexa un número arbitrario pequeno. De aí obtemos que a desigualdade (2.15) verificase con un pequeno arbitrario sempre que tsexa suficientemente grande. Recordemos que o efecto do fármaco anticanceríxeno inhibidor de TGF-βé diminuír o parámetro γque aparece nas ecuacións (2.10) e (2.16). A máxima eficacia do medicamento (ignorando os efectos secundarios negativos) é cando γ= 0, de (2.16) obtemos α(0) = 1 κ1 µTµK2.(2.18) O seguinte teorema da unha condición suficiente baixo a cal o tratamento co inhibidor TGF-βcurará o cancro. Teorema 2.2. Si κTµC> λCα(γ)entón C(t)→0cando t→ ∞. Proba.- Dando conta de que 1−C C0 ≤1, e substituíndo (2.15) na ecuación (2.8), obtemos dC dt ≤λCC−µCTC < (λC−µC κT α(γ)−)C, si t > T3e, por hipótese, o segundo membro é máis pequeno que −βC para algún β > 0, sempre que se elixa o suficientemente pequeno. Resulta que C(t)≤C(T3)e−β(t−T3)t>T3, polo que C(t)→0cando t→ ∞ . O parámetro κTen (2.9) depende da forza do sistema inmunolóxico: vendo a ecuación (2.11), chegamos á conclusión de que o incremento deste parámetro está ligado á movilización das células Tpolos macrófagos para matar as células canceríxenas. O parámetro µCdepende da eficacia das células Ten recoñecer e eliminar as células canceríxenas. Así, en conxunto, o producto κTµCmide a forza do sistema inmunolóxico na loita contra as células do cancro. A función α(γ)depende da eficacia do inhibidor de TGF-β, sendo a máxima eficacia α(0). Polo tanto, o teorema (2.2) dinos que o fármaco por si só non pode garantizar a erradicación do cancro: o sistema inmunolóxico debe ser o suficientemente forte como para que κTµC>λC κ1 µTµK2,(2.19) de maneira que a hipótese do teorema (2.2) podería satisfacerse baixo tratamento con γo suficientemente pequeno como para asegurar que o cancro está erradicado. Sen embargo, 2.2. SIMULACIÓNS NUMÉRICAS 21 observamos que esta hipótese nos da unha condición bastante estricta para a erradicación do cancro; numéricamente pódese derivar unha condición suficiente máis axustada á realidade. Por outro lado, si a desigualdade (2.19) se invirte, entón o teorema (2.2) non pode aplicarse, é dicir, o cancro C(t)pode non desaparecer cando t→ ∞ incluso si γ= 0 (é dicir, baixo calquera tratamento cun fármaco inhibidor de TGF-β). De feito, si asumimos que κTµC<λC κ1 µTµK2.(2.20) e que γ= 0, de tal maneira que M2= 0, entón a ecuación (2.10) deixase fora do modelo (non hai conversión de macrófagos M1a macrófagos M2polo que estes últimos desaparecen), e a ecuación (2.9) convírtese en dM1 dt =k1−µM1. Ademais diso, si (2.20) se verifica e γnon é igual a cero pero é suficientemente pequeno, podemos comprobar facendo un estudo máis amplo [4], que existe unha solución única de estado estacionario da forma: (¯ C+O(γ),¯ M1+O(γ), M2=O(γ),¯ T+O(γ)) onde as funcións O(γ)dependen de γe satisfan |O(γ)| ≤ cte ×γ. Ademais, escribindo explícitamente a matriz Xacobiana neste punto, pódese comprobar que todos os valores propios do polinomio característico son negativos, o que garantiza a estabilidade do punto de equilibrio [4]. 2.2. Simulacións numéricas Imos usar entón o modelo (2.8)-(2.11) para simular canto de efectivo é o medicamento. Empregando os datos de [4]: λC= 10−2/día, µC= 10−5/célula/día, C0= 106célula/cm3, µ= 0.3/día, k1= 3000/cm3/día, γ = 200/día, kT= 3300/célula/día, K1= 0.05C0, K2= 105célula/cm3, µT= 0.2/día. E as seguintes condicións iniciais: C(0)=102célula/cm3, M1(0) = 5 ×104células/cm3, M2(0) = 0, T(0) = 0, 22 CAPÍTULO 2. INTERACCIÓN CANCRO-SISTEMA INMUNOLÓXICO para 0≤t≤400 días. Empregaremos o código Matlab que sigue [4]: principal_inmunoloxico_cancro c l o s e a l l , c l e a r a l l , global lambda_c C_0 mu_c k_1 gamma K_1 K_2 k_T mu mu_T %% parametros lambda_c = 10^(−2) ; % / dia mu_c = 10^(−5) ; % / c e l u l a / dia C_0 = 10^6; % / c e l u l a /cm^3 mu = 0 . 3 ; % / dia k_1 = 3000; % c e l u l a /cm^3/ dia gamma = 200; % / dia mu_T = 0 . 2 ; % / dia k_T = 3300; % / c e l u l a / dia K_1 =0.05∗C_0 ;% c e l u l a /cm^3 K_2 = 10^5; % c e l u l a /cm^3 %% c o nd i ci o ns i n i c i a i s C= 10^2; %c e l u l a /cm^3 M_1 = 5∗10^4; %c e l u l a /cm^3 M_2 = 0 ; T= 0; z_ini = [ C M_1 M_2 T ] ; %co nd ici ons i n i c i a i s tspan = [ 0 , 6 0 ] ; %% solucionador ODE [t,z] = ode15s('funcion_inmunoloxico_cancro ',tspan ,z_ini ) ; %% debuxo tvec = { 'C','M_1','M_2','T'} ; f o r i= 1 : 4 subplot(2 ,2 , i) plot(t,z( : , i) ) , hold on xlabel('t') , ylabel(tvec(i) ) end funcion_inmunoloxico_cancro fun c ti o n dz =funcion_inmunoloxico_cancro(t,z) global lambda_c C_0 mu_c k_1 gamma K_1 K_2 k_T mu mu_T 2.2. SIMULACIÓNS NUMÉRICAS 23 dz=ze ro s (4 ,1 ) ; C=z(1) ; M_1 =z(2) ; M_2 =z(3) ; T=z(4) ; dz (1) = lambda_c ∗C∗(1−C/C_=) −mu_c∗T∗C; dz (2) = k_1 −gamma∗M_1 ∗C/ ( K_1+C)−mu∗M_1 ; dz (3) = gamma∗M_1 ∗C/ ( K_1+C)−mu∗M_2 ; dz (4) = k_T ∗M_1 / ( K_2+M_2 )−mu_T∗T; Figura 2.2: Progreso das células do cancro e das que participan na súa eliminación. Nas gráficas podemos comprobar que se verifican os resultados teóricos anteriormente demostrados. Primeiramente, imos comprobar que os nosos datos cumpren a ecuación (2.19). En efecto kTµC= 3300 ×10−5= 0.033 λC k1 µTµK2=10−2 3000 ×0.2×0.3×105= 0.02 o cal podemos ver tamén gráficamente, porque na primeira gráfica, despois dun tempo t= 400, as células do cancro desaparecen. Isto é debido a que, ademáis de que se cumpre (2.19), γé o suficientemente pequeno como para que se erradique o cancro. Na segunda 30 CAPÍTULO 3. TERAPIAS DO CANCRO Polo tanto ∂ ∂t =∂ ∂T −r˙ R(t) R(t)2 ∂ ∂ρ =∂ ∂T −ρ˙ R(t) R(t) ∂ ∂ρ,∂ ∂r =1 R(t) ∂ ∂ρ. onde se emprega un punto (˙ R)para referirse ás derivadas que involucran o tempo. Isto leva a que 1 r2 ∂ ∂r(r2uf) = u∂f ∂r +f r2 ∂ ∂r(r2u) = u R(t) ∂f ∂ρ +λx −µ(θ−x−y) θf, 1 r2 ∂ ∂r(r2∂v ∂r ) = 1 ρ2R(t)2 1 R ∂ ∂ρ(ρ2R(t)21 R(t) ∂v ∂ρ) = 1 R2 ∂2v ∂ρ2+2 ρR(t)2 ∂v ∂ρ. Definindo agora F=λx −µ(θ−x−y) θ, o sistema (3.1)-(3.6) convírtese en ∂x ∂T +u−ρ˙ R(t) R(t) ∂ ∂ρ(x) = λx −βxv −Fx, (3.19) ∂y ∂T +u−˙ R(t) R(t) ∂ ∂ρ(y) = βxv −δy −Fy, (3.20) n=θ−x−y, (3.21) ∂v ∂T −(ρ˙ R(t) R(t)+2D ρR(t)2)∂v ∂ρ −D R2 ∂2v ∂ρ2=bδy −γv, (3.22) 1 ρ2R ∂ ∂ρ(ρ2u) = F. (3.23) As condicións de contorno fixas e móviles son ˙ R(t) = u(1, t), u(0, t)=0,∂v ∂r (0, t) = 0, e ∂v ∂r (1, t) = 0. Imos entón facer a análise e as simulacións numéricas cun exemplo. Pretendemos calcular R(t)para 0< t < T, T = 15 días para o seguinte conxunto de parámetros e datos iniciais, λ= 2 ·10−2/h, β =7 10 ·10−9mm3/h por virus, δ=1 18/h;µ=1 48/h;D= 3.6·10−2mm2/h; b= 50 virus/célula;γ= 2.5·10−2/h;θ= 106células/mm3 R(0) = 2 mm;x(r, 0) = 0.84 ·106células/mm3;y(r, 0) = 0.1·106células/mm3; v(r, 0) = Ae−r2/4células/mm3, A = 5 ·108. 3.3. SIMULACIÓN NUMÉRICA DO MODELO DE TERAPIA DO CANCRO 31 Para resolver o problema numéricamente é preciso ter en conta que os parámetros teñen un rango moi amplo de valores que poden causar problemas. Polo tanto, eliminamos as dimensións do problema seguindo a seguinte escala: ˜x=x θ,˜y=y θ,˜v=v 102θ,˜ T=T 102h,˜ρ=ρ 1mm,˜u= 102u. Obtemos así as seguintes ecuacións: ∂˜x ∂˜ T+˜u−ρ˙ R(˜ T) R(˜ T) ∂ ∂˜ρ˜x=˜ λ˜x−˜ β˜x˜v−˜ F˜x, (3.24) ∂˜y ∂˜ T+˜u−ρ˙ R(t) R(t) ∂ ∂˜ρ(˜y) = ˜ b˜x˜v−˜ δ˜y−˜ F˜y, (3.25) ∂˜v ∂˜ T−(˜ρ˙ R(t) R(t)+2˜ D ˜ρR(t)2)∂˜v ∂˜ρ−˜ D R2 ∂2˜v ∂˜ρ2=b˜ δ˜y−˜γ˜v, (3.26) 1 ˜ρ2R ∂ ∂˜ρ(˜ρ2˜u) = ˜v, (3.27) onde F=˜ λ˜x−˜µ(1 −˜x−˜y). As condicións de contorno en ˜ρ= 0 e˜ρ= 1 son ˙ ˜ R(˜ t) = u(1,˜ t), u(0,˜ t)=0,∂v ∂˜ρ(0,˜ t) = 0,e∂v ∂˜ρ(1,˜ t)=0, e os parámetros que empregaremos para probar o código serán, seguindo [10]: ˜ λ= 2,˜ β= 0.07; ˜ δ=100 18 ; ˜µ=100 48 ;˜ D= 3.6; b= 50 virus/célula; ˜γ= 2.5; ˜ R(0) = 2; x(r, 0) = 0.84; y(r, 0) = 0.1; v(r, 0) = Ae−r2/4, A = 5. Denotemos a solución numérica no n-ésimo paso por (Xn, Y n, V n, Un, Rn). Empregaremos o método Adams-Bashforth [3] para avanzar no tempo. Para calcular a solución no (n+1)-ésimo paso, o esquema emprega a solución tanto no n-ésimo paso coma no (n−1)-ésimo paso. Aplicando entón o método Adams-Bashforth a ˙ R(˜ t) = u(1,˜ t) obtemos Rn+1 =Rn+∇t 2(3Un−Un−1). 32 CAPÍTULO 3. TERAPIAS DO CANCRO A continuación, a ecuación de tipo hiperbólico, por exemplo, (3.19), resólvese con un esquema dos denominados en inglés “leap-frog integration” [14]: Xn+1 j−Xn−1 j 2∇t+An j Xn j+1 −Xn j−1 2∇ρ=λXn j−βXn jVn j−Fn jXn j. Nos dous extremos, a ecuación convírtese na seguinte ecuación diferencial ordinaria Xn+1 j−Xn−1 j 2∇t=λXn j−βXn jVn j−Fn jXn j. Un problema común asociado a este método é que aplicado a ecuacións non lineais, vólvese inestable debido a que os puntos de malla están completamente desacoplados [10]. Por sorte, esta inestabilidade pódese arreglar acoplando dúas mallas a través do promedio simple no tempo Xn j=1 2(Xn+1 j+Xn−1 j), que mantén a precisión de segunda orde no método. O mesmo esquema emprégase para a ecuación (3.20). Para calcular Un+1, usamos a regra do trapecio ρ2 j+1Un+1 j+1 −ρ2 jUn+1 j=Rn+1 2∇ρ[ρ2 j+1Fn+1 j+1 +ρ2 jFn+1 j]. Finalmente, tratamos a ecuación parabólica (3.22). Escribimos a ecuación na forma ∂v ∂T +A1 ∂v ∂ρ +A2 ∂2v ∂ρ2=RHS. e empregamos para a súa resolución o esquema 3Vn+1 j−4Vn j+Vn−1 j 2∇t+(A1)n+1 j Vn+1 j−1−Vn+1 j−1 2∇ρ+(A2)n+1 j Vn+1 j+1 −2Vn+1 j+Vn+1 j−1 (∇ρ)2= (RHS)n+1 j. Esta ecuación discretizada pódese escribir na forma bjVn+1 j−1+djVn+1 j+ajVn+1 j+1 =Sj, j = 1,2, ..., J, que é un sistema de ecuacións lineais con matriz tridiagonal que debemos resolver. Imos agora entón a mostrar os códigos e as figuras para resolver o problema en Matlab seguindo [10]. O primeiro algoritmo que imos amosar trátase dunha función necesaria para poder implementar o resto dos códigos. Esta función resolve a ecuación matricial Ax =b mediante eliminación gaussiana, onde Aé a matriz tridiagonal, con da diagonal, la diagonal inferior e ua superior; así, dterá nelementos, e leu,n−1elementos cada unha. 3.3. SIMULACIÓN NUMÉRICA DO MODELO DE TERAPIA DO CANCRO 33 Función tridiagSolve fun c ti o n x=tridiagSolve(l,d,u,b) n=length(d) ; % ce r os debaixo da diagon al f o r i=2:n ratio=l(i−1)/d(i−1) ; d(i)=d(i)−ratio ∗u(i−1) ; b (i)=b(i)−ratio ∗b(i−1) ; end x=ze ro s (n, 1 ) ; x(n)=b(n) /d(n) ; % s u bs ti tu t o in ve rs o para r e s o l v e r f o r i=n−1:−1:1 x(i)=(b(i)−u(i)∗x(i+1) ) /d(i) ; end x=x; end Algoritmo Viroterapia N= 48+1; % numero de puntos en x dx = 1/(N−1) ; xx = ( 0 : ( N−1)) '∗ dx ; dt = 0 . 5 ∗(dx ) ^2; Tmax = 15∗24/100; max_iter =ceil(Tmax/dt ) +1; dt =Tmax/( max_iter−1) ; t= ( 0 : ( max_iter −1)) ∗dt ; R=ze r o s (1 , max_iter ) ; U=ze r o s (N, 1 ) ; X=z er os (N, 1 ) ; Y=z e r o s (N, 1 ) ; V=ze r o s (N, 1 ) ; U1 =z e ro s (N, 1 ) ; X1 =zer os (N, 1 ) ; Y1 =z e ro s (N, 1 ) ; V1 =ze r os (N, 1 ) ; U2 =z e ro s (N, 1 ) ; X2 =zer os (N, 1 ) ; Y2 =z e ro s (N, 1 ) ; V2 =ze r os (N, 1 ) ; lambda = 2 . 0 ; plot_iter = 0; col ='bgrcmybgrcmy '; f o r burst = [ 25 50 100 150 2 0 0 ] ; plot_iter =plot_iter+1; beta = 0.07∗burst ;D= 3 . 6 ; delta = 100/18; gamma = 2 . 5 ; mu = 100/48; theta = 1; % i n i c i a l i z a c i o n R(1) = 2 ; %mm en t=0 % p rim eir a i t e r a c i o n iter = 1; X( : ) = 0 . 8 4 ; Y( : ) = 0 . 1 0 ; a= 5 ; V( : ) = ( a∗exp(−xx .^2/4) ) ; U(1) = 0 ; 34 CAPÍTULO 3. TERAPIAS DO CANCRO f o r i= 1 : N−1 Fp = ( lambda∗X(i+1)−mu ∗(theta−X(i+1)−Y(i+1)) ) / theta ; Fm=(lambda∗X(i)−mu ∗(theta−X(i)−Y(i) ) ) /theta ; U(i+1) = (xx (i) ^2∗U(i)+R(iter) /2∗dx ∗(xx (i+1)^2∗Fp+xx (i) ^2∗Fm ) ) /( xx (i+1)^2) ; end R(iter+1) = R(iter)+dt ∗(U(N) ) ; f o r i= 1 : N A= ( U(i)−xx (i)∗U(N) ) /R(iter+1) ; F= ( lambda∗X(i)−mu ∗(theta−X(i)−Y(i) ) ) /theta ; i f i==1 | i==N X1 (i) = X(i)+ dt ∗(lambda∗X(i)−beta ∗V(i)∗X(i)−F∗X(i) ) ; Y1 (i) = Y(i)+ dt ∗(beta ∗V(i)∗X(i)−delta ∗Y(i)−F∗Y(i) ) ; e l s e X1 (i) = X(i)−dt ∗A∗(X(i)−X(i−1) ) /dx+dt ∗(lambda∗X(i)−beta ∗V(i)∗X(i)−F∗X(i) ) ; Y1 (i) = Y(i)−dt ∗A∗(Y(i)−Y(i−1) ) /dx+dt ∗(beta ∗V(i)∗X(i)−delta ∗Y(i)−F∗Y(i) ) ; end end U1 (1) = 0; f o r i= 1 : N−1 Fp = ( lambda∗X1 (i+1)−mu ∗(theta−X1 (i+1)−Y1 (i+1)) ) ; Fm = ( lambda∗X1 (i)−mu ∗(theta−X1 (i)−Y1 (i) ) ) ; U1 (i+1) = ( xx (i) ^2∗U1 (i)+R(iter+1)/2∗dx ∗(xx (i+1)^2∗Fp+xx (i) ^2∗Fm ) ) /( xx (i+1)^2) ; end A1 =−(xx .∗U1 (N) /R(iter+1)+2∗D./(R(iter+1)^2) . / xx ) ; A1 (1) = 0; A2 =−D./ R(iter+1)^2; S= 2∗V+2∗delta ∗dt ∗Y1 ; L1 = 2∗dt/dx /dx ∗(A2 )−dt /dx∗A1 ; D1 = (2−4∗dt /dx ^2∗A2+2∗gamma∗dt )∗ones (N, 1 ) ; UU1 = 2∗dt /dx ^2∗A2+dt /dx ∗A1 ; L1 (N) = 4∗dt/dx /dx ∗(A2 )−dt /dx∗A1 (N) ; UU1 (1) = 4∗dt/dx ^2∗A2+dt /dx ∗A1 (1) ; V1 =tridiagSolve(L1 (2:N) , D1 (1:N) , UU1 (1:N−1) , S) ; f o r iter = 2 : max_iter−1 i f mod (iter , 2 0 )==0 figure(1) ; plot(xx ,X1 ,xx ,Y1 ,xx ,U1 ,xx ,V1 ) legend('X','Y','U','V') t i t l e ( [ 'T = 'num2str ( ( iter−1)∗dt )'R = 'num2str(R(iter ) ) ] ) ; drawnow end R(iter+1) = R(iter) +0.5∗dt ∗(3∗U1 (N)−U(N) ) ; f o r i= 1 : N A= ( U1 (i)−xx (i)∗U1 (N) ) /R(iter+1) ; F= ( lambda∗X1 (i)−mu ∗(theta−X1 (i)−Y1 (i) ) ) ; i f (i==1) | ( i==N) X2 (i) = X(i)+ 2∗dt ∗(lambda∗X1 (i)−beta ∗V1 (i)∗X1 (i)−F∗X1 (i) ) ; Y2 (i) = Y(i)+ 2∗dt ∗(beta ∗V1 (i)∗X1 (i)−delta ∗Y1 (i)−F∗Y1 (i) ) ; V2 (i) = (4∗V1 (i)−V(i)+ 2∗dt ∗delta ∗Y2 (i) ) /(3+2∗gamma∗dt ) ; 3.3. SIMULACIÓN NUMÉRICA DO MODELO DE TERAPIA DO CANCRO 35 e l s e X2 (i) = X(i)−2∗dt ∗A∗(X1 (i+1)−X1 (i−1) ) /2/dx+2∗dt ∗(lambda∗X1 (i)−beta ∗V1 (i)∗X1 (i)−←- F∗X1 (i) ) ; Y2 (i) = Y(i)−2∗dt ∗A∗(Y1 (i+1)−Y1 (i−1) ) /2/dx+2∗dt ∗(beta ∗V1 (i)∗X1 (i)−delta ∗Y1 (i)−F←- ∗Y1 (i) ) ; end end U2 (1) = 0; f o r i= 1 : N−1 Fp = ( lambda∗X2 (i+1)−mu ∗(theta−X2 (i+1)−Y2 (i+1)) ) / theta ; Fm = ( lambda∗X2 (i)−mu ∗(theta−X2 (i)−Y2 (i) ) ) /theta ; U2 (i+1) = (xx (i) ^2∗U2 (i)+R(iter+1)/2∗dx ∗(xx (i+1)^2∗Fp+xx (i) ^2∗Fm ) ) /( xx (i+1)^2) ; end A1 =−(xx .∗U2 (N) /R(iter+1)+2∗D./(R(iter+1)^2) . / xx ) ; A1 (1) = 0 ; A2 = −D./ R(iter+1)^2; S= 4∗V1−V+2∗delta ∗dt ∗Y2 ;L1 = 2∗dt /dx/dx ∗A2−dt/dx ∗A1 ; D1 = (3−4∗dt /dx ^2∗A2+2∗gamma∗dt )∗ones (N, 1 ) ; UU1 = 2∗dt /dx ^2∗A2+dt /dx ∗A1 ;L1 (N) = 4∗dt/dx /dx ∗(A2 )−dt /dx∗A1 (N) ; UU1 (1) = 4∗dt/dx ^2∗A2+dt /dx ∗A1 (1) ; V2 =tridiagSolve(L1 (2:N) , D1 (1:N) , UU1 (1:N−1) , S) ; X=X1 ;Y=Y1 ;U=U1 ;V=V1 ;X1 =X2 ;Y1 =Y2 ;U1 =U2 ;V1 =V2 ; i f mod (iter , 5 0 ) == 0 X1 = 0 . 5 ∗(X+X2 ) ; Y1 = 0. 5∗(Y+Y2 ) ; U1 = 0. 5 ∗(U+U2 ) ; V1 = 0 .5 ∗(V+V2 ) ; end end figure(2) ; hold on ;plot(t∗100/24 , R,col (plot_iter) ) ; t i t l e ('R') save ( [ 'burst_ 'num2str(burst )'. mat '] ) end Na primeira figura mostrase Xa densidade de células canceríxenas, Ya densidade de células canceríxenas infectadas polo virus, Ua velocidade radial, e Va densidade de virus libres. Executando o código en Matlab, podemos ver que conforme vai aumentando o tempo T, o radio do tumor Rvai diminuindo considerablemente, asi como o número de virus e de células canceríxenas infectadas polo virus ou non. Na segunda podemos ver a diminución do radio do tumor en 15 días. 36 CAPÍTULO 3. TERAPIAS DO CANCRO Figura 3.1: Diminución do radio do cancro en 15 días. Conclusións No primeiro capítulo viamos un modelo do cancro basado no modelo depredador-presa. Para iso, consideramos as células normais e canceríxenas como dúas poboacións en competencia deste tipo. Como conclusión observamos que, sen tratamento, as células canceríxenas acabarán por encher todo o tecido. No seguinte capítulo estudamos a interacción entre as células do cancro e o sistema inmunolóxico, o encargado da defensa natural do corpo contra os axentes infecciosos. Obtendo as condicións que debe exercer o sistema inmunolóxico sobre o cancro para que así este poida ser erradicado baixo tratamento co fármaco anticanceríxeno inhibidor de TGF- β. De non cumprirse esas condicións, o cancro pode non desaparecer baixo tratamento cun fármaco deste tipo. Por último, modelamos dous tipos de terapias do cancro, a viroterapia e o tratamento co fármaco GM-CSF. E fixemos simulacións numéricas do caso da viroterapia para poder ver como afectaría este tratamento. 37 Apéndice Ecuacións diferenciais Lembraremos aquí a teoría das ecuacións diferencias necesaria para o tratamento dos modelos presentados ao longo do traballo. Para iso, seguiremos a estrutura do libro Introduction to mathematical biology [4] para as tres primeiras seccións; e a obra Introducción a las Ecuaciones en Derivadas Parciales (EDP’s) [18] para proporcionar información sobre as ecuacións en derivadas parciais na última sección. Ecuación diferencial de primeira orde As ecuacións diferenciais de primeira orde (EDO’s) teñen a forma dx dt =f(x, t)(1) onde f(x, t)é unha función dada. De feito, existen numerosas solucións deste tipo, pero a solución será única si preescribimos unha condición inicial x(t0) = x0(2) para os valores x0ey0. O conxunto das ecuacións (1)-(2) denomínase problema de valor inicial. Existen diferentes tipos de ecuacións diferenciais que se poden resolver explícitamente Ecuacións lineais dx dt +p(t)x=g(t) onde p(t)eg(t)son funcións de t. Separación de variabeis dx dt =g(x)h(t). 39 46 Apéndice: Ecuacións diferenciais Figura 2: Retrato de fases do sistema (12). Para resolver unha ecuación lineal non homoxénea sd2x dt2+bdx dt +cx =f(t) con unha función dada f(t), primeiro precisamos encontrar unha solución especial ˜xe, logo, a solución xeral é a suma de ˜xe a solución xeral da ecuación homoxénea. O mesmo procedemento se aplica ós sistemas lineais non homoxéneos. Sistemas de dúas ecuacións diferenciais Imos estudar sistemas xerais de dúas ecuación diferencias de primeira orde, dx1 dt =f1(x1, x2),dx2 dt =f2(x1, x2)(21) onde f1(x1, x2), f2(x1, x2)son funcións dadas, non necesariamente lineais. Un punto (a, b) tal que f1(a, b) = 0, f2(a, b) = 0. chámase punto de equilibrio, punto estacionario ou punto estable do sistema (21). A línea nula x1é a curva que consta dos puntos que satisfan a ecuación f1(x1, x2) = 0. Apéndice: Ecuacións diferenciais 47 De maneira similar,a línea nula x2é a curva que consta dos puntos que satisfan a ecuación f2(x1, x2) = 0. Os puntos de equilibrio do sistema (21) son os puntos onde as dúas lineas nulas se intersecan. Para ter unha idea de cómo se comportan as traxectorias cerca dun punto estacionario (a, b), precisamos linealizar o sistema. Establecemos X1=x1−a, X2=x2−b. Entón, pola fórmula de Taylor, para i= 1,2 fi(x1, x2) = fi(a+X1, b +X2) = fi(a, b) + ∂fi ∂x1 X1+∂fi ∂x2 X2+termos de orde superior, onde abreviamos ∂fi ∂x1 =∂fi ∂x1 (a, b),∂fi ∂x2 =∂fi ∂x2 (a, b). Si definimos aij =∂fi ∂xj (a, b) entón o sistema (21) ten a forma dXi dt =ai1X1+ai2+termos de orde superior (i= 1,2) onde X1, X2están próximos a 0. Polo tanto, as traxectorias de (a, b)espérase que se comporten dunha maneira similar á das traxectorias de dXi dt =ai1X1+ai2X2, i = 1,2.(22) En consecuencia disto, o punto de equilibrio (a, b)de (21) dise que é asintóticamente estable (ou, brevemente, estable) se o punto de equilibrio x= 0 de (22) é asintóticamente estable, é dicir, se as partes reais dos valores propios da matriz A=aij son negativas. Un punto de equilibrio que non é estable dise que é inestable. Concluímos que o punto de equilibrio (a, b)do sistema (21) é estable si e só si as seguintes desigualdades se manteñen en (a, b): ∂f1 ∂x1 +∂f2 ∂x2 <0,∂f1 ∂x1 ∂f2 ∂x2 −∂f1 ∂x2 ∂f2 ∂x1 >0, é dicir, a traza de (∂fi ∂xj)<0e o determinante de (∂fi ∂xj)>0. A matriz (∂fi ∂xj(a, b)) chámase matriz Xacobiana no punto de equilibrio (a, b). 48 Apéndice: Ecuacións diferenciais Ecuaciones en derivadas parciais Chamamos ecuación diferencial en derivadas parciais (EDP) á ecuación da forma F(x1, x2, ..., xn, u, ∂u ∂x1 , ..., ∂u ∂xn , ..., ∂mu ∂k1x1∂k2x2...∂knxn ) = 0 (23) que permite conectar as variabeis independientes xi,∀i= 1,2, ..., n, a función que se busca e as súas derivadas parciais. Cúmplese que: ki,∀i= 1,2, ..., n son enteiros non negativos tales que: k1+k2+... +kn=m. A función Fé a función prefixada dos seus argumentos. Chámase orde dunha EDP á orde superior das derivadas parciais que figuran na ecuación. Empréganse as notacións ux≡∂u ∂x, uy≡∂u ∂y , uxx ≡∂2u ∂x2, ... Dada a EDP definida en (23) de orden m, chámase solución de dita EDP en certa rexión Dde variación das xi,∀i= 1,2, ..., n a unha función calquera u=u(x1, x2, ..., xn)∈ Cm(D) (conxunto das funcións continuas na rexión Dxunto con todas as derivadas de ata orde m incluidas), tal que ó substituir u, e as súas derivadas en (23), a última se convirte na identidade respecto a xi,∀i= 1,2, ..., n na rexión D. Dada unha EDP de orde n, unha solución que conteña nfuncións arbitrarias chámase solución xeral, e calquera solución obtida desta solución xeral por seleccións particulares das funcións arbitrarias chámase solución particular. En moitas ocasións precisamos determinar solucións de EDP’s que satisfagan condicións dadas. Notación Nas ecuacións diferenciais en derivadas parciais é moi común denotar as derivadas parciais empregando subíndices: u0 x=∂u ∂x, u00 xy =∂2u ∂y∂x =∂ ∂y ∂u ∂x Especialmente na física matemática, soese preferir o operador nabla (que en coordenadas cartesianas se escribe como ∇= (∂x, ∂y, ∂z)para as derivadas parcias e un punto ( ˙u)para as derivadas que involucran o tempo. Ecuacións en derivadas parciais lineais A ecuación en derivadas parciais chámase lineal, si esta é lineal respecto á función buscada e todas as súas derivadas que forman parte da ecuación. En caso contrario, chámase non lineal. Apéndice: Ecuacións diferenciais 49 A EDP lineal de segunda orde para a función de dúas variabeis independentes xeyno caso xeral ten a forma A(x, y)∂2u(x, y) ∂x2+ 2B(x, y)∂2u(x, y) ∂x∂y +C(x, y)∂2u(x, y) ∂y2+(24) a(x, y)∂u(x, y) ∂x +b(x, y)∂u(x, y) ∂y +c(x, y)u(x, y) = f(x, y), sendo A(x, y), B(x, y), C(x, y), a(x, y), b(x, y), c(x, y)funcións das variabeis xeynunha rexión D⊂R2, e a función incógnita u=u(x, y). Se f(x, y) = 0 en D⊂R2, a ecuación (24) chámase homoxénea. Dada a EDP de segunda orde (24) nunha certa rexión contida Den R2(plano OXY ), dise que é unha función Hiperbólica en , si δ=B2−AC > 0en D, Parabólica en , si δ=B2−AC = 0 en D, Elíptica en , si δ=B2−AC < 0en D. Bibliografía [1] Armitage, Peter and Doll, Richard The Age Distribution of cancer and a multi-stage theor of carcinogenesis, 1954. [2] Bernoulli, Daniel Essai d’une nouvelle analyse de la mortalité causée par la petite vérole, et des avantages de l’inoculation pour la prévenir. Histoire de l’Acad., Roy. Sci.(Paris) avec Mem, Des Math. And Phis., Mem, 1-45, 1760. [3] Butcher, John Charles and Goodwin, Nicolette Numerical methods for ordinary differential equations Wilwy Online Library, 2, 2008. [4] Chou, Ching-Shan and Friedman, Avner, Introduction to Mathematical Biology. Springer, 2008. [5] Díaz, José and Álvarez, Elena, Breve historia de las biomatemáticas en los siglos XX y XXI. Inventio, la génesis de la cultura universitaria en Morelos 4(7): 61–68. 2008. [6] Di Prima, Boyce and Boyce, W Ecuaciones diferenciales y problemas con valores en la frontera Limusa, Grupo Noriega Editores, 2000. [7] Dublin, Loots I and Lotka, Alfred J On the true rate of natural increase: As exemplified by the population of the United States, 1920 Journal of the American statistical association, Taylor & Francis Group, 20(151): 305-339, 1925. [8] Euler, Leonhard Découverte d’un nouveau principe de mécanique. Mémoires de l’académie des sciences de Berlin, 185-217, 1752. [9] Fibonacci, Leonardo, Liber abacci. 1202. [10] Friedman, Avner and Kao, Chiu-Yen Mathematical modeling of biological processes, Springer, 2014. [11] Gómez Gómez, Marta Modelos matemáticos en oncología. Simulación numérica (traballo fin de grao). Universidad Zaragoza, 2014 51 52 BIBLIOGRAFÍA [12] Lindstrom, Mary J and Bates, Douglas M Nonlinear mixed effects models for repeated measures data, Biometrics, JSTOR, 673-687, 1990. [13] Malthus, Thomas R An essay on the theory of population, Oxford: Oxford University Press, 1798. [14] Morton, K.W., Mayers, D. F. Numerical solution of Partial Differential Equations Cambridge University Press, 2, 2005. [15] Ozores, Antón Lombardero, Un vistazo a la Biomatemática. Números: Revista de didáctica de las matemáticas, Sociedad Canaria Isaac Newton de Profesores de Matemáticas, 86: 29–38. 2014. [16] Rashevsky, Nicolas and others Mathematical biophysics, University of Chicago Press, 1938. [17] Ribba, Benjamin and Kaloshi, Gentian and Peyre, Mathieu and Ricard, Damien and Calvez, Vincent and Tod, Michel and Čajavec-Bernard, Branka and Idbaih, Ahmed and Psimaras, Dimitri and Dainese, Linda and others A tumor growth inhibition model for low-grade glioma treated with chemotherapy or radiotherapy, Clinical Cancer Research, AACR, 18(18): 5071-5080, 2012. [18] Romero, Sixto e J.Moreno, Francisco e M.Rodríguez, Isabel Introducción a las Ecuaciones en Derivadas Parciales (EDP’s) Universidad de Huelva, 2001. [19] Turing, Alan M Philosophical Transactions of the Royal Society of London, Series B, Biological Sciences, 237(641): 37-72, 1952. [20] Verhulst, Pierre-François Notice sur la loi que la population suit dans son accroissement Corresp. Math. Phys., 10: 113-126, 1838. [21] Volterra, Vito Variazioni e fluttuazioni del numero d’individui in specie animali conviventi Societá anonima tipografica “Leonardo da Vinci”, 1926. [22] https://www.cancer.gov/ [23] http://www.ecologia.unam.mx/web/index.php/investigacion/ ecologia-evolutiva [24] https://www.icmat.es/newsletter/2015/eighth_es.pdf [25] https://seom.org/informacion-sobre-el-cancer/que-es-el-cancer-y-como-se-desarrolla BIBLIOGRAFÍA 53 [26] https://es.wikipedia.org/wiki/Biologia_matematica [27] https://es.wikipedia.org/wiki/Biologia_de_sistemas [28] https://fr.wikipedia.org/wiki/Biomathematique [29] https://en.wikipedia.org/wiki/Mathematical_and_theoretical_biology [30] https://es.wikipedia.org/wiki/Neurociencia_computacional [31] https://es.wikipedia.org/wiki/Simetria_radial_(Biologia)