UNIVERSIDAD COMPLUTENSE DE MADRID

GRADO DE MATEMÁTICAS

MODELOS DE INTERACCIÓN ENTRE CÉLULAS
TUMORALES Y EL SISTEMA INMUNE

TRABAJO DE FIN DE GRADO

Realizado por: Carlos Gallardo Arroyo

Dirigido por: Ana Carpio Rodŕıguez

Departamento de Matemática Aplicada

Curso 2016-2017



A mi familia y amigos,
que me enseñaron a no rendirme.

ii



Índice general

1. Interacciones entre el Sistema Inmune y Cáncer. Modelos matemáticos sin
dependencia espacial. 1
1.1. Introducción [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2. Teoŕıa del estudio cualitativo de los sistemas dinámicos. 3
2.1. Diagrama de fases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2. Puntos cŕıticos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3. Bifurcaciones [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.3.1. Bifurcaciones tipo nodo-silla. . . . . . . . . . . . . . . . . . . . . . . . 6
2.3.2. Bifurcaciones tipo transcŕıtico. . . . . . . . . . . . . . . . . . . . . . . 6
2.3.3. Bifurcación tipo Pitchfork. . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3.4. Bifurcaciones tipo Hopf. . . . . . . . . . . . . . . . . . . . . . . . . . 8

3. Modelo de una ecuación: Crecimiento tumoral. 10
3.1. Modelo de Gompertz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.2. Descripción del modelo de una ecuación diferencial. . . . . . . . . . . . . . . 11

3.2.1. Ejemplo del modelo loǵıstico. . . . . . . . . . . . . . . . . . . . . . . 12

4. Modelos de ecuaciones: Interacciones entre células tumorales y sistema
inmune. 15
4.1. Desarrollo del modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.2. Modelo teórico de Kuznestov (1994). . . . . . . . . . . . . . . . . . . . . . . 16

4.2.1. Estados estacionarios. . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.2.2. Bifurcaciones en el sistema. . . . . . . . . . . . . . . . . . . . . . . . 19

4.3. Simulación numérica del modelo de Kurznestov. . . . . . . . . . . . . . . . . 21
4.3.1. Eliminación del tumor sin recidivas. . . . . . . . . . . . . . . . . . . . 21
4.3.2. Crecimiento del tumor. . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.3.3. Crecimiento y decrecimiento tumoral periódico. . . . . . . . . . . . . 24

4.4. Modelo teórico de Sotolongo-Costa (2003). . . . . . . . . . . . . . . . . . . . 25
4.4.1. Estados estacionarios. . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.4.2. Dinámica y análisis de la estabilidad [14]. . . . . . . . . . . . . . . . . 25
4.4.3. Interpretación biológica [14]. . . . . . . . . . . . . . . . . . . . . . . . 27

4.5. Simulación numérica del modelo de Sotolongo-Costa. . . . . . . . . . . . . . 29
4.5.1. Descontrol del sistema inmune. . . . . . . . . . . . . . . . . . . . . . 29
4.5.2. Control de las células malignas. . . . . . . . . . . . . . . . . . . . . . 31
4.5.3. Crecimiento y decrecimiento peródico significativo. . . . . . . . . . . 32

5. Modelo de tres ecuaciones. Interacciones entre células canceŕıgenas y dos
componentes del ambiente tumoral (S.I) 33
5.1. Crecimiento tumoral controlado por dos tipos de células efectoras. . . . . . . 34

5.1.1. Simulación numérica del modelo Moore-Li. . . . . . . . . . . . . . . . 34

iii



ÍNDICE GENERAL ÍNDICE GENERAL

5.2. Crecimiento tumoral controlado por células efectoras y células normales. . . 37
5.3. Crecimiento tumoral controlado por células efectoras y citoquinas. . . . . . . 37

5.3.1. Simulación numérica del modelo Krischner - Panneta. . . . . . . . . . 39

6. Discusión sobre los modelos de ODEs estudiados y posibles enfoques fu-
turos. 41
6.1. Incremento de la complejidad de los sistemas. . . . . . . . . . . . . . . . . . 41
6.2. Posibles elementos a incorporar en los modelos. . . . . . . . . . . . . . . . . 42

6.2.1. Discretización y modelos estocásticos. . . . . . . . . . . . . . . . . . . 42
6.2.2. Agentes individuales y microambiente tumoral. . . . . . . . . . . . . 43
6.2.3. Estructura espacial. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

7. Anexo I

iv



Abstract

Resumen

Este trabajo de fin de grado es un acercamiento a modelos dinámicos basados en ecua-
ciones diferenciales con los cuales poder modelizar la interacción entre las células tumorales
y el sistema inmune. Haremos un recorrido por distintos modelos planteados a lo largo de los
últimos años de una, dos y tres ecuaciones estudiando algunas de sus caracteŕısticas, tanto
de forma cualitativa como anaĺıtica. También haremos algunas simulaciones numéricas que
nos ayudarán a entender cómo funcionan estos modelos.

Palabras clave:

Cáncer, sistema inmune, citoquinas, ecuaciones diferenciales, bifurcaciones, simulación,
diagrama de fases.

Abstract

This end-of-degree paper is an approach to dynamic models based on differential equa-
tions with which to model the interaction between tumor cells and the immune system. We
will take a tour of some different models proposed over the last years of one, two and three
equations studying some of their characteristics, both qualitatively and analytically. We will
also do some numerical simulations that will help us understand how these models work.

Key words:

Cancer, differential equations, inmune system, cytokinas, bifurcations, simulation, phase
portrait.

v



vi Abstract

vi



CAPÍTULO 1.

Interacciones entre el Sistema Inmune y Cáncer.

Modelos matemáticos sin dependencia espacial.

1.1. Introducción [1].

En los últimos años se ha observado que el sistema inmune es capaz de reconocer y eliminar
células tumorales.

Gran parte de las investigaciones se han centrado en cómo mejorar la respuesta anti-tumoral
estimulando el sistema inmune a través de vacunas o directamente inyectando células T o
citoquinas. Sin embargo, para crear inmunoterapias efectivas tenemos que entender primero
los mecanismos que gobiernan el crecimiento tumoral.

Los primeros pasos que se dieron para entender el crecimiento tumoral fueron en aras de
entender cómo las células sanas mutaban a células canceŕıgenas. Sin embargo este concepto
de células sanas no es adecuado ya que esto no depende únicamente de la célula en cuestión,
sino de la interacción entre múltiples factores, lo que se conoce como microsistema tumoral:
células inmunológicas, fibroblastos, tejido conectivo, celulas endoteliales, citoquinas, hormo-
nas... Las interacciones entre las células cancerosas y otros componentes están continuamente
cambiando ya que la fuerza de interacción depende de la densidad (o concentración).

Una forma de aproximarnos al entendimiento de cómo se desarrollan estas interacciones es a
través de modelos matemáticos incluso si algunos de los componentes de los modelos no son
del todo conocidos. Gracias a estos modelos podemos intuir los mecanismos más plausibles
que tienen estos procesos a través de las observaciones.

La mayoŕıa de los modelos están centrados en modelos espaciales descritos por ecuaciones
en derivadas parciales (PDEs) o autómatas celulares [2]. Muchos de estos modelos (basados
en la teoŕıa cinética) también utiliza ecuaciones integrales y diferenciales [3] [5].

Sin embargo, podemos describir estos procesos centrándonos en una perspectiva no espacial,
a través de ecuaciones diferenciales (ODEs). Esto nos permite simplificar los modelos para
poder centrarnos en las interacciones entre las células cancerosas y los diferentes tipos de
células inmunes. Para ello podemos utilizar una, dos o más ecuaciones según los elementos
que pretendamos incorporar al modelo.

Con todo esto veremos las ventajas e inconvenientes de utilizar las ecuaciones diferenciales
para el estudio de este tipo de modelos.

Las inmunoterapias se han centrado principalmente en la actividad anti-tumoral de los glóbu-

1



2 1.1. Introducción [1].

los blancos, especialmente las células T (CD8+), los ”natural killer”(NK) y los macrófagos.
Estos tipos de células inmunológicas son capaces de provocar la lisis en las células diana
(es decir, provocar la rotura de la membrana celular). Además se pueden incorporar células
tisulares o citoquinas aśı como terapias externas.

Este trabajo comienza con un resumen de ciertos conceptos teóricos sobre teoŕıa cualitativa
de ecuaciones diferenciales aśı como de teoŕıa de bifurcaciones los cuales nos ayudarán a
entender los modelos propuestos por los autores.

En cuanto al estudio de los modelos, empezamos planteando una ecuación referida al creci-
miento de la población de células tumorales. Esto nos va a permitir a parte de entender cómo
evoluciona un sistema de este tipo, qué otros factores necesitamos introducir en el modelo
para estudiar este tipo de sistemas cáncer-inmune.

Cuando hayamos estudiado las posibilidades que tiene los sistemas de una ecuación, introdu-
cimos otra a través de la cual describimos el funcionamiento de un tipo de células del sistema
inmune genéricas. De esta forma, podemos entender cómo interactúan estas células con las
células que conforman los tumores. Esta sección es la más amplia y la que se estudiará con
detalle.

Seguido a esto, podemos incorporar otra ecuación más (añadiendo otra variable nueva) para
modelizar la interacción, no solo entre los dos tipos de células anteriores, sino con otros tres
posibles tipos de células o agentes qúımicos. En este caso únicamente veremos cómo plan-
tear estos modelos y una simulación numérica con parámetros experimentales para ver cómo
funcionaŕıan en el caso de tratar con casos prácticos reales.

Por último, se proponen lineas de continuación de este trabajo utilizando distintas técni-
cas (algunas fuera del análisis de los sistemas de ecuaciones diferenciales) pero que seŕıan
igualmente claves para ser capaces de tener una visión global de la enfermedad que estamos
tratando.

2



CAPÍTULO 2.

Teoŕıa del estudio cualitativo de los sistemas

dinámicos.

Se considera un sistema de n ecuaciones diferenciales que se puede escribir de la forma:

dy

dt
= f(y, t)

de forma que y ∈ Rn , t ∈ R y f : Rn × R → Rn. Entonces vamos a ver una serie de
definiciones y propiedades que nos servirán para los modelos.

Definición 2.1. Se dice que un sistema de ecuaciones diferenciales es un sistema autónomo
si la función f no depende expĺıcitamente de la variable independiente t.

Propiedades 2.1. Todo sistema de ecuaciones se puede transformar en un sistema autóno-
mo añandiendo un cambio de variable.

Definición 2.2. Se denomina espacio de estados del sistema al subconjunto de Rn (Ω) en
el que existen soluciones del sistema.

2.1. Diagrama de fases.

Queremos crear un modelo geométrico del conjunto de todos los posibles estados del sistema
en el espacio de estados. En este modelo geométrico la dinámica determina una estructu-
ra celular de cuencas limitadas por separatrices. En cada cuenca existe un núcleo llamado
atractor.

Entonces el objetivo del marco teórico es la descripción del diagrama de fases, estudian-
do el comportamiento asintótico de las trayectorias, determinando trayectorias especiales y
atractores.

Definición 2.3. Se denomina diagrama de fase del sistema de ecuaciones diferenciales y′ =
f(y) al conjunto de todas las trayectorias de las soluciones del sistema.

Definición 2.4. Un atractor de un sistema de ecuaciones diferenciales es un conjunto cerra-
do y acotado hacia el cual se aproxima, cuando el tiempo t→∞ la órbita de las soluciones.

El diagrama de fases consta de trayectorias que no se cortan por lo que pueden ser de tres
tipos:

Curvas abiertas y simples.

Curvas cerradas y simples (soluciones peródicas).

Puntos de equilibro (soluciones constantes).

3



4 2.2. Puntos cŕıticos.

2.2. Puntos cŕıticos.

Se estudian los puntos cŕıticos o estacionarios del sistema autónomo que son las soluciones
de f(y) = 0. Cada punto cŕıtico es una órbita del sistema formada sólo por dicho punto.

Definición 2.5. Una solución estacionaria de un sistema de ecuaciones diferenciales es una
solución constante, es decir, y(t) = y(t0).

Definición 2.6. Un punto y0 ∈ Rn es un punto cŕıtico del sistema autónomo y′ = f(y) si
f(y0) = 0. En caso contrario se dice que es un punto regular.

Ahora vamos a ver la clasificación de los puntos cŕıticos:

Definición 2.7.

v es un nodo estable si para cada ε > 0 existe un δ > 0 tal que si y es solución tal que
‖y(t0)− v‖ < δ entonces ‖y(t)− v‖ < ε ∀t > t0

v es un nodo asintóticamente estable si v es estable y existe un δ > 0 tal que para cada
solución y, si ‖y(t0)− v‖ < δ, entonces ĺımt→∞ ‖y(t)‖ = v. ES UN ATRACTOR.

v es un nodo inestable si no es estable. ES UN REPULSOR.

Por otra parte, algunos sistemas ecuaciones dependen de un parámetro (o más):

dy

dt
= f(y, λ)

donde los parámetros no dependen de la variable independiente pero se modifica el com-
portamiento de las soluciones según los valores que tome el parámetro. En general, una
modificación pequeña del valor del parámetro supone una variación también pequeña de
las soluciones. Sin embargo, existen ocasiones en las que una pequeña variación produce
un cambio drástico en la dinámica del sistema. En estas ocasiones se dice que la ecuación
diferencial tiene una bifurcación.

Definición 2.8. Una ecuación diferencial que depende de un parámetro tiene una bifurca-
ción cuando existe una variación cualitativa en el comporamiento de las soluciones para un
determinado valor del parámetro.

Dos teoremas importantes que utilizaremos son:

Teorema 2.2. (Green): Sea Ω ∈ C(R,R) cerrado y sea D el interior de C (acotado). Si
existen dos funciones f, g sobre un abierto que contenga a D y tienen parciales continuas
entonces: ∫

Ω

(f(x, y)dx+ g(x, y)dy) =

∫ ∫
D

(
∂f

∂x
− ∂g

∂y

)
dxdy

Teorema 2.3. (Dulac-Bendixon): Sea una función ϕ(x, y) ∈ C1 de forma que satisface
la siguiente condición:

∂(ϕf)

∂x
+
∂(ϕg)

∂y
tiene el mismo signo a.e.(x,y) de un conexo de R2

4



2. Teoŕıa del estudio cualitativo de los sistemas dinámicos. 5

entonces el sistema autónomo dado por:{
dx
dt

= f(x, y)
dy
dt

= g(x, y)

no tiene soluciones periódicas.

Demostración.
Sin pérdida de generalidad vamos a suponer que existe una función ϕ(x, y) se verifica:

∂(ϕf)

∂x
+
∂(ϕg)

∂y
> 0

en una región conexa del plano denotada por Ω.

Ahora sea C una trayectoria cerrada del plano autónomo en Ω (es decir, estamos suponiendo
que existe una solución orbital del sistema). Y ahora sea D el interior del conjunto C, es
decir:

D = int(C)

Entonces aplicamos el teorema de Green:∫ ∫
D

(
∂(ϕf)

∂x
+
∂(ϕg)

∂y

)
dxdy = ...

... =

∫
C

(−ϕgdx+ϕfdy) =

∫
C

(−ϕgdx+ϕfdy) =

∫
C

ϕ

(
−dy
dt
dx+

dx

dt
dy

)
=

∫
C

ϕ(−ẏdx+ẋdy)

pero como en Ω tenemos que:

dx = f(x, y)dt = ẋdt

dy = g(x, y)dt = ẏdt

luego sustituimos y tenemos que el valor de la integral es:∫
C

ϕ(−ẏdx+ ẋdy) =

∫
C

ϕ(−ẏẋ+ ẋẏ)dt = 0

pero la integral no puede anularse porque estamos integrando en una órbita cerrada (debeŕıa
ser positiva). Luego no pueden existir órbitas.

2.3. Bifurcaciones [4].

El comportamiento de las soluciones del modelo puede variar con los parámetros de los que
depende el mismo, es decir, para determinados valores cŕıticos puede variar completamente
la dinámica. Dichos valores cŕıticos de los parámetros se denominan bifurcaciones.
Hay varios tipos de bifurcaciones que aparecen con frecuencia. Estudiaremos algunos casos.

5



6 2.3. Bifurcaciones [4].

2.3.1. Bifurcaciones tipo nodo-silla.

Este tipo de bifurcación proporciona el mecanismo por el cual los puntos de equilibrio se
crean o se destruyen.

El prototipo de bifurcación nodo-silla lo define la ecuación:

x′ = r + x2

donde r ∈ R es el parámetro. Entonces, dependiendo del valor de r tenemos:

Si r > 0 no hay equilibrios.

Si r = 0 hay un único equilibrio (donde se encuentra la bifurcación).

Si r < 0 hay dos equilibrios (
√
−r y −

√
−r).

La ecuación diferencial pasa de tener un equilibrio que es atractor y otro que es repulsor a
una situación sin equilibrios donde las trayectorias se van a infinito.

Figura 2.1: Ejemplo de diagrama de bifurcación tipo: nodo-silla [4]

2.3.2. Bifurcaciones tipo transcŕıtico.

La bifurcación transcritica proporciona un mecanismo para el cambio de estabilidad de un
punto de equilibrio que persiste al variar los parámetros. Sin embargo śı puede cambiar de
estabilidad.

6



2. Teoŕıa del estudio cualitativo de los sistemas dinámicos. 7

La forma normal para una bifurcación transcŕıtica es:

x′ = rx− x2

donde x, r pueden ser positivos o negativos. Vemos la naturaleza de los equilibrios en función
del parámetro:

Si r < 0 hay un equilibrio estable en x = 0 y otro inestable en x = r

Si r = 0 los dos equilibrios colisionan en x = 0 que es marginalmente estable, es decir,
es estable o inestable dependiendo de según por donde nos acerquemos.

Si r > 0 hay un equilibrio inestable en x = 0 y otro estable en x = r

Figura 2.2: Ejemplo de diagrama de bifurcación tipo: transcŕıtico [4]

2.3.3. Bifurcación tipo Pitchfork.

Existen dos tipos de bifurcaciones Pitchfork: supercŕıtica y subcŕıtica.

La bifurcación supercŕıtica tiene la forma normal:

x′ = rx− x3

que en función del parámetro tenemos:

Si r < 0 hay un equilibrio exponencialmente estable en x = 0

Si r = 0 hay un equilibrio estable en x = 0

7



8 2.3. Bifurcaciones [4].

Figura 2.3: Ejemplo de diagrama de bifurcación tipo: pitchfork supercŕıtico [4]

Si r > 0 entonces x = 0 se hace inestable y aparecen dos equilibrios estables simétricos
en x = ±

√
r

La bifurcación subcŕıtica tiene la forma normal:

x′ = rx+ x3

que en función del parámetro tenemos:

Si r > 0 hay un equilibrio inestable en x = 0

Si r = 0 hay un equilibrio estable o inestable en x = 0 dependiendo de por donde nos
aproximemos.

Si r < 0 entonces x = 0 se hace inestable y aparecen dos equilibrios inestables simétricos
en x = ±

√
−r

2.3.4. Bifurcaciones tipo Hopf.

Las bifurcaciones en torno a puntos de equilibrio estudiadas hasta ahora persisten en siste-
mas dinámicos de dimensión mayor. Las variaciones de interés se concentran en un espacio
unidimensional, mientras que el resto de las componentes contribuyen a la atracción o re-
pulsión a ese subespacio.

En sistemas dinámicos bidimensionales a parte de los puntos de equilibrio podemos tener
trayectorias cerradas:

8



2. Teoŕıa del estudio cualitativo de los sistemas dinámicos. 9

Figura 2.4: Ejemplo de diagrama de bifurcación tipo: pitchfork subcŕıtico [4]

Estables: atraen las trayectorias situadas fuera y dentro.

Inestables: repelen las trayectorias situadas fuera y dentro.

Marginalmente estables: atraen las trayectorias de un lado pero no del otro.

Las bifurcaciones tipo Hopf proporcionan el mecanismo por el que se crea un ciclo ĺımite a
partir de un estado estacionario. Éstas son un fenómeno no lineal.

Figura 2.5: Ejemplo de diagrama de bifurcación tipo: hopf [4]

9



CAPÍTULO 3.

Modelo de una ecuación: Crecimiento tumoral.

El primer paso para comprender el crecimiento del tumor es describir los patrones de cre-
cimiento que presentan este tipo de sistemas. Por ello se utiliza gran cantidad de datos
procedentes de pruebas de rayos X y mamograf́ıas (entre otros) para ver como son estos
patrones.

Lo que los datos muestran es que el crecimiento del tumor (antes de que sea posible el
diagnóstico por imagen) es bastante más rápido que la proliferación de las células malignas
cuando el tumor alcanza cierto tamaño y ya es subceptible de ser detectado a través de
pruebas rutinarias. En este aspecto los datos muestran que hasta el 77 % de los pacientes
con cáncer de mama experimentan la mayor parte del crecimiento entre dos pruebas anuales.

3.1. Modelo de Gompertz.

El primero en explicar este tipo de crecimiento fue Gompertz (1825) a través de modelos
matemáticos de replicación y muerte celular (y nada más).

Este modelo proporciona una curva sigmoidal de crecimiento de la población que muestra
una aceleración en el crecimiento en los primeros instantes y una desaceleración en el creci-
miento para grandes poblaciones celulares.

En el caso del desarrollo tumoral este tipo de dinámica puede ser explicada por la cantidad
de recursos de los que disponen las células, siendo mayores cuanto menor es la población (y
de ah́ı su rápida proliferación) y agotándose cuanto más grande es el tumor y mas recursos
necesita para expandirse.

Descripción del modelo:

El modelo teórico es un modelo matemático que en su principio y su final sufre un crecimien-
to lento pero tiene un desarrollo rápido. En referencia al crecimiento canceŕıgeno la primera
fase de crecimiento lento se obviaŕıa y se estudiaŕıa el sistema a partir de que el crecimiento
comience a ser significativo.

La ecuación del modelo se puede expresar de la siguiente forma:

N(t) = kC−ae
−bt

t ≥ 0

tal que:

N(t): Tamaño de la población para cierto instante de tiempo t.

10



3. Modelo de una ecuación: Crecimiento tumoral. 11

k: Constante positiva que indica el tamaño máximo que puede alcanzar la población.

b: Constante relacionada con la capacidad de proliferación de la población.

Propiedades del modelo de Gompertz:

Posee un valor inicial cuando t = −∞ y el valor k cuando t = +∞

A diferencia del modelo loǵıstico esta curva no es simétrica.

El punto de inflexión de la curva se encuentra en N = an0

Ejemplo 3.1. Una curva de crecimiento podŕıa estar descrita de la siguiente manera:

N(t) = kC−ae
−bt

de forma que a = Ln

(
k

N0

)
condición inicial del problema.

Figura 3.1: Crecimiento de la población para distintos parámetros

3.2. Descripción del modelo de una ecuación diferencial.

Desde el trabajo de Gompertz muchos modelos matemáticos han sido derivados del descrito
anteriormente para explicar los datos del crecimiento tumoral, predecir la supervivencia del
paciente y ser capaces de sugerir el mejor tratamiento adaptado al tipo de tumor del paciente.

Algunos de estos modelos asumen que el crecimiento es exponencial pero la mayoŕıa conside-
ran distintos tipos de crecimiento que deceleran con el tiempo (cuando t→∞). La ecuación
general que describe la dinámica de crecimiento tumoral puede escribirse como [1]:

dx

dt
= xf(x)

11



12 3.2. Descripción del modelo de una ecuación diferencial.

Figura 3.2: Esquema de interacción de las células canceŕıgenas x

donde x es el tamaño de la población en un tiempo t, x(0) > 0 y el factor f(x) es relativo a
la dependencia de la densidad en la proliferación y muerte de las células tumorales.

La dependencia de la densidad la podemos escribir como:

f(x) = p(x)− d(x)

donde p(x) describe la proliferación celular y d(x) describe la muerte celular.

Los modelos más simples pueden expresarse a través de potencias o funciones potenciales:

p(x) = axα

d(x) = bxβ

Si tomamos:

α = 0 y β = 1 (Modelo loǵıstico) [7]

α = −1/3 y β = 0 (Modelo de Bertanffly [6])

p(x) = a y d(x) = bLn(x) (Modelo de Gompertz) [8]

3.2.1. Ejemplo del modelo loǵıstico.

Para el modelo loǵıstico reescalamos la función multiplicado por 1/b, b 6= 0 entonces tenemos:

F (x) = rx− x2

donde r = a
b

En este caso tenemos que:

Si r < 0⇔ a
b
< 0:

x = 0 ESTABLE

x = a/b NO ESTABLE

12



3. Modelo de una ecuación: Crecimiento tumoral. 13

Si r > 0⇔ a
b
> 0

x = 0 NO ESTABLE

x = a/b ESTABLE

Luego en 0 tenemos una bifurcación transcŕıtica.

Figura 3.3: Caso k < 0

Figura 3.4: Caso k > 0

Estos modelos como hemos dicho se pueden utilizar para el crecimiento tumoral lo que nos
permite ver cómo afecta el tratamiento, si hay recidivas en el tiempo...

Por ejemplo el modelo de Gompertz tiene un crecimiento no acotado cuando la densidad
tiende a cero:

ĺım
x→0+

f(x) = +∞

Entonces para grandes poblaciones de células está eventualmente acotado pero este modelo
no sirve para describir tumores pequeños (esto es una limitación bastante grande).

13



14 3.2. Descripción del modelo de una ecuación diferencial.

Conocemos soluciones anaĺıticas para la ecuación antes descrita tomando cualesquiera paráme-
tros α y β lo que hace que sea sencillo predecir la dinámica de los tumores dada una medida
del tamaño del tumor para cierto tiempo (y estimando los parámetros del modelo).

De todas maneras, estos modelos han dado buenos resultados en la explicación de los patro-
nes de crecimiento [10]. Además esto los hace idóneos para clasificar los tumores según su
agresividad y cuantificar la relación entre el crecimiento y la edad del paciente.

Con esta ecuación podemos modelar los siguientes casos:

Crecimiento del tumor p(x) > d(x)

Decrecimiento del tumor p(x) < d(x)

Latencia del tumor p(x) = d(x)

Sin embargo una sola ecuación no puede modelar una situación en las que estén combinados
los tres anteriores ya que el signo de xf(x) siempre es el mismo. No puede haber crecimiento
en un tiempo t1 y decrecimiento en otro t2.

En cambio, esto se puede solventar incluyendo en el modelo un término de tratamiento
dependiente del tiempo φ(t):

dx

dt
= x[p(x)− d(x)]− aφ(t)x

donde a > 0 representa la fuerza del agente administrado y φ(t) representa la concentración
del agente durante el tratamiento. También puede ser interpretado como respuesta inmune.

De todas maneras este modelo aún no puede explicar los mecanismos que lideran la inter-
acción que pretendemos estudiar. Para ello, lo que vamos a estudiar en el siguiente caṕıtulo
es la interacción con otras células, las células inmunes, las cuales son capaces de regular
y controlar (a veces) el crecimiento del tumor. Esto lo hacemos introduciendo una nueva
variable y una nueva ecuación al sistema.

14



CAPÍTULO 4.

Modelos de ecuaciones: Interacciones entre células

tumorales y sistema inmune.

4.1. Desarrollo del modelo.

La hipótesis de supervivencia formulada en los años 50 dice que el sistema inmune es capaz
de inhibir el crecimiento de tumores pequeños y eliminarlos antes de que sean cĺınicamente
detectables. Esto nos permite establecer el siguiente modelo matemático de dos ecuaciones
diferenciales de primer orden. Para ello, vamos a ver cómo interactuan estas células:

Descripción del modelo [11] [12]:

Las células tumorales provocan un efecto positivo en las células inmunes, es decir,
provocan un est́ımulo en su crecimiento descrito con el término py(x, y).

Las células tumorales también provocan un efecto negativo en las células inmunes, es
decir, provocan un descenso de su población. Esto viene dado por el término dy(x, y).

Las células inmunes provocan un efecto negativo en las células canceŕıgenas, es decir,
provocan el descenso de su población descrito con el término dx(x, y).

Además, las células canceŕıgenas provocan sobre ellas mismas un efecto tanto positivo
como negativo lo que viene descrito por el término xf(x).

A parte de esto, hay que tener en cuenta las posibles interacciones externas (al ciclo cáncer-
sistema inmune) provocadas por la actuación por ejemplo de quimioterapia u otras con-
diciones descritas con un término φy(t) las cuales dependen del tiempo de actuación o de
exposición.

Lo último que vamos a tener en cuenta por ahora es la muerte de células inmunes descrita
con el término ay(y).

Considerando todo esto, lo que vamos a hacer es aplicar el modelo presa-depredador de
ecuaciones diferenciales para estudiar la relación entre las células canceŕıgenas (x) y las
células inmunes (y): {

x′ = xf(x)− dx(x, y)

y′ = py(x, y)− dy(x, y)− ay(y) + φy(t)

donde las células tumorales actúan como presa y las células inmunes actúan como depreda-
dor. En este sistema x representa la densidad o la cantidad de células tumorales e y representa

15



16 4.2. Modelo teórico de Kuznestov (1994).

Figura 4.1: Esquema de las interacciones entre poblaciones x e y

la densidad o cantidad de células inmunes.

Entonces la población x depende de la autorregulación impuesta por xf(x) y por la actua-
ción del sistema inmune sobre estas células.

La población y, sin embargo, depende por un lado de la producción de células inmunes pro-
vocadas por la aparición de células tumorales (a modo de est́ımulo) py(x, y) y del aporte
externo de distintas terapias φ(t) pero, por otro lado, se ve reducida por la actuación de las
células canceŕıgenas en dos términos (directo e indirecto) dy(x, y) y ay(y).

Vamos a empezar a estudiar una serie de modelos concretos para las ecuaciones planteadas.

4.2. Modelo teórico de Kuznestov (1994).

Modelo concreto propuesto por Kuznetsov [13] de funciones para nuestro sistema de ecua-
ciones diferenciales:

f(x) = a(1− βx)
dx(x, y) = nxy
py(x, y) = ρxy

g+x

dy(x, y) = mxy
ay(y) = dy
φ(t) = s

En este sistema propuesto vamos a utilizar el criterio de Dulac-Bendixon para estudiar las
posibles órbitas del sistema. Para ello necesitamos previamente el Teorema de Green.

16



4. Modelos de ecuaciones: Interacciones entre células tumorales y sistema
inmune. 17

Entonces aplicando el teorema de Dulac-Bendixon el resultado al modelo tenemos que, si
consideramos el intput externo nulo φ(t) = 0. Entonces si tomamos la función ϕ(x) = 1

xy
por

el teorema de Dulac-Bendixon tenemos que:

∂

∂x

(
ϕ
dx

dt

)
+

∂

∂y

(
ϕ
dy

dt

)
= −

(
s

xy2
+
aβ

y

)
Entonces si asumimos que los parámetros son positivos es siempre negativo luego no hay
ciclos ĺımite o conexiones homocĺınicas en el sistema. En esta misma linea tenemos también
que carece de bifurcaciones de Hopf.

4.2.1. Estados estacionarios.

Nos interesan únicamente las soluciones de las poblaciones x e y no negativas. De la misma
forma también vamos a asumir que los parámetros son no negativos. Entonces estudiamos
las curvas que verifican:

dx

dt
= 0 y

dy

dt
= 0

Entonces tenemos estados estacionarios en la intersección de ambas.
Para dx

dt
= 0 tenemos:

x = 0

y = a(1− βx) ≡ g(x)

Para dy
dt

= 0 tenemos:

y =
s

d+mx− ρx
g+x

≡ f(x)

Entonces para la intersección de x = 0 y f(x) obtenemos el punto:

(x, y) =
(

0,
s

d

)
Ahora para la intersección de las funciones f(x) y g(x) vemos que obtendŕıamos un polinomio
de grado 3:

C3x
3 + C2x

2 + C1x+ C0 = 0

donde:
C0 = g

(s
a
− d
)

C1 =
s

a
+ ρ−mg − d+ dmβ

C2 = −m+ (mg + d− ρ)β C3 = mβ

Por lo que obtendŕıamos mas estados estacionarios (de 1 a 3 dependiendo del valor de los
parámetros). Veamos algunas caracteŕısticas de las funciones f y g para determinar esto:

La función g es siempre una recta con pendiente −aβ

Dos aspectos de la función f :

Propiedades 4.1. La función f tiene aśıntotas horizontales cuando el parámetro ρ verifica
alguna de las condiciones siguientes:

ρ >
(√

(gm) +
√
d
)2

ó ρ <
(√

(gm)−
√
d
)2

17



18 4.2. Modelo teórico de Kuznestov (1994).

Demostración.
Tenemos aśıntotas horizontales cuando el denominador se anula, entonces:

d+mx− ρx

g + x
= 0

xasintota =
(ρ− d− gm±

√
(d+ gm− ρ)2 − 4gdm)

2m

Entonces tenemos dos soluciones cuando (d + mg − ρ)2 > 4gdm. Por lo tanto nos resultan
las condiciones del principio.

Propiedades 4.2. La función f tiene un máximo en x > 0 cuando ρ > gm

Demostración.
Tomando derivada de f tenemos que:

m =
ρg

(g + x)2

entonces...

xmáx =

√
ρg

m
− g

Figura 4.2: Distintas posibilidades de la función f en función de lo estudiado [13]

Sin embargo, para los parámetros estimados por los autores estamos en el caso en el que
existen 4 equilibrios (dos estables y dos inestables). Si los calculamos numéricamente con el

18



4. Modelos de ecuaciones: Interacciones entre células tumorales y sistema
inmune. 19

polinomio antes descrito (para los valores concretos de los parámetros) aśı como ls autovalores
asociados (con el jacobiano) tenemos:
Equilibrio inestable (silla):

A = (0. 32, 0) λ1 = −0. 374 λ2 = 1. 321

Equilibrio estable (espiral):

B = (1. 6, 8. 19) λ1 = −0. 501 + 0. 573i λ2 = −0. 501− 0. 573i

Equilibrio inestable (silla):

C = (0. 76, 267. 8) λ1 = −1. 356 λ2 = 0. 325

Equilibrio estable (nodo):

D = (0. 17, 447. 13) λ1 = −1. 693 λ2 = −0. 453

Figura 4.3: Diagrama de fases para el conjuto de parámetros dado por los autores.

4.2.2. Bifurcaciones en el sistema.

Vamos ahora a estudiar las posibles bifurcaciones respecto al parámetro s. Escogemos es-
te parámetro ya que se refiere a la tasa del tratamiento, por lo cual, en la práctica, es el
parámetro que vamos a poder modificar y ajustar su valor para obtener los mejores resulta-
dos posibles en cuanto a la reducción del tumor en el paciente se refiere [9].

Existen dos bifurcaciones que las vemos en el siguiente diagrama (figura 4.4).

La bifurcación (1) es de tipo nodo-silla mientras que la bifurcación (2) es de tipo transcŕıtico.

Entonces podemos observar los diagramas de fases dependiendo del valor de s para observar
los cambios en la dinámica [9] [13]. Estos valores del parámetro los tomaremos en las tres
posibles regiones:

19



20 4.2. Modelo teórico de Kuznestov (1994).

Figura 4.4: Parámetro bifuraciones (s)

1. Antes de la bifurcación (1) En este caso es en el que se encuentran los parámetros
estimados.

2. Entre la bifurcación (1) y (2)

3. Después de la bifurcación (2)

Figura 4.5: Caso 1, s = 0,118

20



4. Modelos de ecuaciones: Interacciones entre células tumorales y sistema
inmune. 21

Figura 4.6: Caso 2, s = 0,3

Figura 4.7: Caso 3, s = 1

4.3. Simulación numérica del modelo de Kurznestov.

En esta sección vamos a ver cómo se desarrollaŕıa el tumor en distintos pacientes los cua-
les presentan caracteŕısticas distintas. Esto se ve reflejado en los parámetros del modelo
obtenidos de la experimentación (Kuznestov [13]):

4.3.1. Eliminación del tumor sin recidivas.

Los parámetros que tomamos son:

21



22 4.3. Simulación numérica del modelo de Kurznestov.

d = 0,1908 β = 2x10−3 ρ = 1,131 g = 20,19
m = 0,00311 a = 1,636 n = 1 s = 0,318

Figura 4.8:

Figura 4.9:

Las células tumorales decaen a 0 rápidamente por lo que no hay evidencias de cáncer. Las
células inmunes se mantienen constantes pero con densidad baja ya que no hay agentes
externos (como nuevas células canceŕıgenas) que lo activen.

22



4. Modelos de ecuaciones: Interacciones entre células tumorales y sistema
inmune. 23

4.3.2. Crecimiento del tumor.

Los parámetros que tomamos son:

d = 2,0 β = 4x10−3 ρ = 1,131 g = 20,19
m = 0,00311 a = 1,636 n = 1 s = 0,318

Figura 4.10:

Figura 4.11:

En este caso las células tumorales crecen rápidamente mientras que las células inmunes son
incapaces de eliminarlo y caen rápidamente.

23



24 4.3. Simulación numérica del modelo de Kurznestov.

4.3.3. Crecimiento y decrecimiento tumoral periódico.

Los parámetros que tomamos son:

d = 0,3743 β = 2x10−3 ρ = 1,131 g = 20,19
m = 0,00311 a = 1,636 n = 1 s = 0,318

Figura 4.12:

Figura 4.13:

En este caso podemos ver como la densidad e células canceŕıgenas aumenta para tiempo
t = 15 y va disminuyendo aunque con oscilaciones. No llega a desaparecer el tumor si no que
se queda en torno al valor x = 10.

24



4. Modelos de ecuaciones: Interacciones entre células tumorales y sistema
inmune. 25

4.4. Modelo teórico de Sotolongo-Costa (2003).

El modelo propuesto por los autores es el siguiente [14]:{
x′ = αx− xy
y′ = xy − 1

α
y − kx+ σ + V cos2(βt)

En este caso, el factor tratamiento se incorpora como una función sinusoidal, lo que nos
permitiŕıa modelizar la incorporación de diversas terapias que aumentan y disminuyen en
el tiempo (como por ejemplo ocurre con las quimioterapias aplicadas cada cierto periodo de
tiempo cuyo efecto está en función del lapso de tiempo que haya pasado desde la aplicación).

4.4.1. Estados estacionarios.

Igual que en el modelo anterior vamos a ver qué puntos verifican:

dx

dt
= 0 y

dy

dt
= 0

Para dx
dt

= 0:
x = 0

y = α

Para dy
dt

= 0:

y =
kx− σ
x− 1

α

Por lo tanto tenemos los puntos:
(x, y) = (0, ασ)

y

(x, y) =

(
1− σ − V cos2(βt)

α− k
, α

)

4.4.2. Dinámica y análisis de la estabilidad [14].

Para estudiar la dinámica de este sistema vamos a tomar V = 0, es decir, un sistema en el
que no incluimos tratamiento, solo consideramos las interacciones entre células.

dx

dt
= αx− xy

dy

dt
= xy − 1

α
− kxσ

con x(0) = x0 e y(0) = y0 condiciones iniciales.

Sustituyendo la primera ecuación en la segunda:

d2x

dt2
+

(
1

α
− x− 1

x

dx

dt

)
= (k − α)x2 + αx

25



26 4.4. Modelo teórico de Sotolongo-Costa (2003).

donde x(0) = x0 y x′(0) = v0 como condiciones iniciales.
El potencial es:

U(x) = −1

3
(k − α)x3 − 1

2
σx2

El potencial tiene dos extremos:
x1 = 0

x2 =
σ − 1

k − α
Si solo tomamos los casos en que x1 > 0 y x2 > 0 (para que tenga sentido) tenemos los
extremos:

x1 = 0 minimo:

σ > 1⇒ k

α
> 1 máximo

x1 = 0 máximo:

σ < 1⇒ k

α
> 1 mı́nimo

El análisis de los puntos fijos en el plano de fases de las ecuaciones muestran dos estados
estacionarios. Vemos los puntos fijos y los autovalores asociados [14]:

P1 = (0, ασ)

λ± =
α2(1− σ)− 1

2α
± |α

2(1− σ) + 1

2α
|

Si σ < 1 es un punto de silla.

Si σ > 1 es un nodo estable.

P2 =

(
1− σ
α− k

, α

)

λ± =
k − ασ

2α(α− k)
±

√(
k − ασ

2α(α− k)

)2

− (1− σ)

El análisis de los autovalores nos da la siguiente tabla de casos:

26



4. Modelos de ecuaciones: Interacciones entre células tumorales y sistema
inmune. 27

Figura 4.14:

4.4.3. Interpretación biológica [14].

Vamos a empezar por el caso en que k
α
< 1:

Para σ < k
α

se produce un descontrol del crecimiento tumoral. Si k
σ
< 1 el sistema evoluciona

con una masa de células tumorales de forma oscilatoria amortiguada. Esto también se puede
interpretar como un estado latente del tumor.

Figura 4.15: (a): Descontrol del crecimiento tumoral; (b) Crecimiento amortiguado

En ambos casos casos existe una población de células tumorales. Sin embargo, el primer caso
se ha visto que es posible para cualquier tipo de condición inicial mientras que el segundo
sólo se da cuando hay una respuesta inmune débil.

27



28 4.4. Modelo teórico de Sotolongo-Costa (2003).

Figura 4.16: Plano de fases: células malignas frente linfocitos σ = 0,09

Figura 4.17: Plano de fases: células malignas frente linfocitos σ = 0,25

Ahora vemos el caso contrario, cuando k
α
> 1:

En este caso hay dos posibilidades para los valores de σ. Puede ser 1 < k
α
< σ y 1 < σ < k

α
.

En ambos casos existen puntos de silla lo que significa que si la población de células can-
ceŕıgenas está por debajo de la linea que separa las dos dinámicas (separatriz) la situación
es irreversible (color rosa) y el cáncer se desarrolla hasta que la respuesta inmune se puede
considerar nula. Esto puede que no cause la muerte del paciente pero si que el sistema in-
mune sea incapaz de responder a otras enfermedades.

Sin embargo, si está por encima de la separatriz (color amarillo), los experimentos cĺınicos
han demostrado que puede haber una regresión del tumor.

28



4. Modelos de ecuaciones: Interacciones entre células tumorales y sistema
inmune. 29

Figura 4.18: k
α
> 1

La simulación numérica nos permite ver la existencia de una ”barrera inmunológica”, es
decir, la barrera hasta la cual el sistema inmune es capaz de eliminar el cáncer, pero que
si es sobrepasada se pueden generar tumores grandes. También, en función de otro tipo de
modelos y parámetros que no hemos estudiado puede aparecer cierto tipo de ”barrera in-
munológica”bajo la cual el sistema inmune es capaz de controlar la densidad de población y
reducirla al mı́nimo pero sin ser capaz de eliminarla.

El análisis de estos modelos de dos ecuaciones también revela que el modelo de crecimiento
de Gompertz no es compatible con la hipótesis de inmuno-supervivencia puesto que el siste-
ma inmune no puede erradicar completamente las células canceŕıgenas que crecen según el
modelo Gompertziano, aunque como dijimos, pueda ajustarse bien a los datos.

4.5. Simulación numérica del modelo de Sotolongo-Costa.

Vemos cómo pueden ser las trayectiorias de las distintas componentes (población células
malignas y población de células inmunes) para el modelo estudiado previamente propuesto
por Sotolongo-Costa:

4.5.1. Descontrol del sistema inmune.

α = 2,0 k = 0,2 σ = 0,05
x0 = 2,1 y0 = 2,7 V = 0 β = 0,34

29



30 4.5. Simulación numérica del modelo de Sotolongo-Costa.

Figura 4.19: Caso 1 (Sotolongo-Costa)

Figura 4.20: Caso 1 (Sotolongo-Costa)

En este caso vemos como el sistema inmune se activa de manera exagerada para controlar
las células malignas. Se produce un estado periódico en el que el ataque del sistema inmune
puede dañar al paciente (enfermedad autoinmune).

30



4. Modelos de ecuaciones: Interacciones entre células tumorales y sistema
inmune. 31

4.5.2. Control de las células malignas.

α = 2,0 k = 0,2 σ = 0,05
x0 = 5,3 y0 = 6,7 V = 0,25 β = 0,34

Figura 4.21: Caso 2 (Sotolongo-Costa)

Figura 4.22: Caso 2 (Sotolongo-Costa)

En este caso, el sistema inmune es capaz de adaptarse a la actividad tumoral y erradicarla
cuando se produce un pico de existencia de células canceŕıgenas. El resultado es una estabili-
zación de la masa tumoral. Además el término de tratamiento no nulo también hace posible
esta estabilización.

31



32 4.5. Simulación numérica del modelo de Sotolongo-Costa.

4.5.3. Crecimiento y decrecimiento peródico significativo.

α = 2,0 k = 0,2 σ = 0,05
x0 = 5,3 y0 = 6,7 V = 0,25 β = 0,32

Figura 4.23: Caso 3 (Sotolongo-Costa)

Figura 4.24: Caso 3 (Sotolongo-Costa)

Al incluir el término de tratamiento al sistema tenemos estos crecimientos y decrecimientos
peródicos. En el momento en el que la densidad de población de células malignas crece el
sistema inmune responde y actúa el fármaco o terapia aplicada al paciente.

32



CAPÍTULO 5.

Modelo de tres ecuaciones. Interacciones entre células

canceŕıgenas y dos componentes del ambiente tumoral

(S.I)

Hemos visto en el caṕıtulo anterior las interacciones entre dos elementos del ambiente tu-
moral: las células canceŕıgenas y células generéricas del sistema inmune. Sin embargo este
tipo de modelo, a pesar de ser bastante eficiente a la hora de clasificar tipos de tumores,
ver su desarrollo, etc, podemos complicarlo para que tenga otro tipo de interacciones que si
que se sabe que se dan, por ello será nuestro punto de partida para añadir una nueva ecuación.

El componente que vamos a incorporar podŕıa ser de naturaleza distinta: o bien otro tipo
diferente de célula inmune o procedente de tejido sano o bien citoquinas presentes en el
microambiente celular.

A pesar de que estos tres tipos de interacción tienen el mismo resultado (la lisis de la mem-
brana celular de las células cancerosas) el modo de modelizar estas interacciones es distinto
debido a que la forma en que se produce esta lisis y, en definitiva, la forma de interactuar
con las células es distinta.

De forma resumida, en este caṕıtulo vamos a ver los siguientes tres modelos:

1. Interacciones entre células canceŕıgenas y dos tipos de células efectoras (como pueden
ser las células T CD8+ y los nátural Killers NK) [16].

2. Interacciones entre células canceŕıgenas, células efectoras y células tisulares normales
[18].

3. Interacciones entre células cancerosas, células efectoras y citoquinas (como pueden ser
TGF-β, IFN-γ) [17].

Sin embargo hay otros modelos de tres ecuaciones en los que se investigan las interacciones
entre células cancerosas, células efectoras y células ”naive”. También los hay que estudian
celulas cancerosas y distintos tipos de anticuerpos. Sin embargo es imposible abarcar todos
los posibles modelos.

La mayoŕıa de estos están basados en ODEs como los que estudiamos en este trabajo. Sin
embargo seŕıa interesante incorporar efectos estocásticos en las interacciones pues nos daŕıan
una visión más realista de lo que puede ocurrir.

También incluiremos en estos modelos distintos tratamientos como la inyección continua de
citoquinas, de células efectoras, o de distintos tipos de medicamentos.

33



34 5.1. Crecimiento tumoral controlado por dos tipos de células efectoras.

5.1. Crecimiento tumoral controlado por dos tipos de células efec-
toras.

El modelo general que simula las interacciones entre las células canceŕıgenas y dos poblaciones
de células distintas seŕıa:

x′ = xf(x)− dx(x, y, z)

y′ = py(x, y, z)− dy(x, y)− ay(y) + φy(t)

z′ = pz(x, y, z)− dy(x, z)− az(z) + φz(t)

Como en los modelos anteriores, f(x) describe el crecimiento de células tumorales sometido
a cierta autoregulación. El factor dx(x, y, z) proporciona el ratio de muerte de células tumo-
rales debidas a las interacciones entre los otros dos tipos de células, y y z.

En la segunda y tercera ecuación, los términos py(x, y, z) y pz(x, y, z) nos proporcionan el
crecimiento de las células y y z respectivamente en presencia de las células canceŕıgenas x.
Los términos ay(y) y az(z) describen la muerte celular (apoptosis), mientras que los térmi-
nos dy(x, y) y dz(x, z) describen la capacidad de las células tumorales de inactivar las células
efectoras (término de competición entre estas células: x e y aśı como z e y).

Los dos factores que restan, φy(t) y φz(t) son los términos referidos a los tratamientos que
aplicamos al paciente, los cuales son dependientes del tiempo (como ocurre en una terapia
real) a pesar de que puedan ser distintos en la forma de actuar.

Como hemos visto en la introducción de este caṕıtulo, este sistema de ecuaciones modela
las interacciones entre células cancerosas y dos tipos de elementos presentes en el ambiente,
escogidas en este caso las células CD8+ T (y) y células T näıve (z).

Modelo de Moore - Li (2003) [15]
x′ = rcxLn

(
Cmax

x

)
− dcx− γcxy

y′ = αnknz
(

x
x+η

)
+ αey

(
x

x+η

)
− dey − γext

z′ = sn − dnz − knz
(

x
x+η

)

5.1.1. Simulación numérica del modelo Moore-Li.

Tomamos el siguiente conjunto de valores para los parámetros del modelo:

sn = 0,071 dn = 0,05 de = 0,12 dc = 0,68
kn = 0,063 η = 43; αn = 0,56 αe = 0,93

Cmax = 19 · 104 γe = 1,9 · 10−3 γc = 0,048

Obteniendo las siguientes gráficas:

34



5. Modelo de tres ecuaciones. Interacciones entre células canceŕıgenas y dos
componentes del ambiente tumoral (S.I) 35

Figura 5.1: rc = 0,23

Figura 5.2: x(0) = 104 y(0) = 50 z(0) = 2500

Figura 5.3: rc = 0,23

35



36 5.1. Crecimiento tumoral controlado por dos tipos de células efectoras.

Figura 5.4: x(0) = 104 y(0) = 30 z(0) = 2500

Figura 5.5: rc = 0,43

Figura 5.6: x(0) = 104 y(0) = 50 z(0) = 2500

36



5. Modelo de tres ecuaciones. Interacciones entre células canceŕıgenas y dos
componentes del ambiente tumoral (S.I) 37

5.2. Crecimiento tumoral controlado por células efectoras y células
normales.

En muchos casos es importante, a parte de estudiar las células inmunes que producen un efec-
to claro en los tumores, estudiar las células sanas de los distintos tejidos que se encuentran
en torno al tumor. Esto se debe a que la competición por los recursos y el tratamiento que se
le pueda aplicar a los pacientes puede afectar gravemente al tejido sano por lo que, a parte
de querer maximizar el número de células malignas que eliminamos queremos, casi con el
mismo nivel de importancia, minimizar la destrucción de células sanas o tejidos involucrados.

La competición entre estos tres tipos de células la podemos resumir en los siguientes aspectos:

1. Competición por los recursos entre células normales y células cancerosas (e incluso
algunas células inmunológicas).

2. Competición de tipo predador-presa entre cáncer y sistema inmune.

Este tipo de interacciones pueden ser descritas de forma parecida al modelo de la sección
anterior donde, ahora, la población z representa la densidad de células normales sanas/ te-
jido. En particular podemos presentar el siguiente modelo:

Modelo de Pillis y Radunskaya (2001 - 2003)[16]

f(x) = r1(1− b1x)
ay(y) = d1y
az(z) = 0

dy(x, y) = c1xy
dz(x, z) = c4xz
φy(t) = s
φz(t) = 0

py(x, y, z) = ρxy
α+x

pz(x, y, z) = r2z(1− b2z)

5.3. Crecimiento tumoral controlado por células efectoras y citoqui-
nas.

Uno de los primeros modelos en investigar el rol de actuación de las cytokinas en la regresión
tumoral fue descrito por Krischner y Panetta, quienes investigaron el efecto de IL-2 y las
células T citotóxicas en la dinámica de interacción del sistema inmune con el tumor.

Las ecuaciones que describen la dinámica entre las células cancerosas (x), células inmunes
(y) y las citoquinas (z) puede ser descrita de forma general como:

37



38 5.3. Crecimiento tumoral controlado por células efectoras y citoquinas.


x′ = xf(x)− dx(x, y, z)

y′ = py(x, y, z)− ay(y) + φy(t)

z′ = pz(x, y, z)− az(z) + φz(t)

Particularizamos con el siguiente modelo:

Modelo Krischner - Panetta (1998)[17]

f(x) = r2(1− bx)
ay(y) = µ2y
az(z) = µ3z
dy(x, y) = 0
dz(x, z) = 0

dx(x, y) = axy
g2+x

φy(t) = s1

φz(t) = s2

py(x, y, z) = cx+ p1yz
g1+z

pz(x, y, z) = p2xy
g3+x

Este modelo pretende investigar los resultados dados dos tipos de tratamientos continuos:

1. La inyección de citoquinas (s1 > 0).

2. La transferencia de células inmunes ya sea de un paciente a otro o la re-inyección de
sus propias células (s2 > 0).

Los resultados anaĺıticos y numéricos muestran que en ausencia del tratamiento, es decir,
cuando s1 = s2 = 0 el modelo puede mostrar tres tipos de comportamientos:

Persistencia de tumores macroscópicos.

Oscilación entre tumores macroscópicos y microscópicos.

Persistencia de tumores latentes.

Entonces, ahora si nos focalizamos en la actuación de los tratamientos tenemos que:

Si s1 > 0 y s2 = 0 entonces puede mantenerse un estado de ausencia tumoral estable.
El valor cŕıtico del tratamiento se alcanza en s1 > s1

crit = r2g2µ2
a

.

Si s1 = 0 y s2 > 0 donde el valor cŕıtico del tratamiento está en s2
crit = µ2µ3g1

p1−µ2 . En estos
casos también puede darse la posibilidad de que el sistema inmune tenga una respuesta
exagerada (y →∞) lo que causaŕıa daño al paciente.

Si aplicamos los dos tratamientos al mismo tiempo, s1 > 0 y s2 > 0 entonces mejoran
sustancialmente los resultados aunque el mejor tratamiento en este caso se obtiene
cuando s2 está por debajo de su valor cŕıtico. En dicho caso podemos mantener un
estado libre de actividad tumoral.

38



5. Modelo de tres ecuaciones. Interacciones entre células canceŕıgenas y dos
componentes del ambiente tumoral (S.I) 39

5.3.1. Simulación numérica del modelo Krischner - Panneta.

Los parámetros tomados en este modelo serán:

p1 = 0,1245 p2 = 5 g1 = 2 · 107 g2 = 105 g3 = 103

a = 1 b = 10−9 r2 = 0,18 µ2 = 0,03 µ3 = 10

Obteniendo los siguientes patrones de crecimiento para cada valor de c:

Figura 5.7: c = 0,05

Figura 5.8: x0 = 10 y0 = 100000 z0 = 0

39



40 5.3. Crecimiento tumoral controlado por células efectoras y citoquinas.

Figura 5.9: c = 0,02

Figura 5.10: x0 = 10 y0 = 100000 z0 = 0

40



CAPÍTULO 6.

Discusión sobre los modelos de ODEs estudiados y

posibles enfoques futuros.

6.1. Incremento de la complejidad de los sistemas.

Hemos empezado el trabajo considerando el modelo más simple de ecuaciones (el dado por
una única ecuación diferencial ordinaria). Este modelo sólo es capaz de describir un tipo de
población celular (en este caso se ha escogido las células canceŕıgenas) en donde el proceso
esta sujeto a la autoregulación y a la dependencia de la densidad de población.

Estos modelos son capaces, a pesar de su simplicidad, de ajustar gran cantidad de los datos
recogidos en distintos experimentos aśı como también son útiles a la hora de predecir el
comportamiento del tumor en ausencia de tratamiento.

Sin embargo, como hemos visto, con estos modelos no podemos estudiar los comportamien-
tos que tienen ciertos tumores de crecimiento y de recidiva del tumor, aśı como, obviamente,
no nos permite estudiar las distintas interacciones con agentes presentes en el medio celular
que los experimentos dicen que se producen.

Por lo tanto, este tipo de modelos nos sirve como acercamiento a la idea de ser capaces
de modelizar procesos que ni si quiera hoy son comprendidos en su totalidad por la ciencia
médica.

En los siguientes caṕıtulos, a través del conocimiento biológico que tenemos de estos sistemas
hemos sido capaces de incrementar el grado de detalle.

Considerando los distintos agentes que interactúan con las células malignas como pueden
ser células con capacidad inmune, células sanas procedentes de tejidos colindantes, agentes
como citoquinas, natural killers, etc podemos crear modelos que, como antes, aunque no estén
totalmente descritos por las ciencias experimentales, los modelos matemáticos nos ayudan a
entender todas las posibles interacciones que se pueden producir entre todos ellos. Esto es
importante por dos razones:

Ser capaces de entender el desarrollo de la enfermedad en un grupo de pacientes para
aśı establecer las ventajas de distintos tratamientos frente a otros aśı como la predicción
del crecimiento o receso del tumor en función del tiempo.

Focalizar las nuevas investigaciones médicas en aspectos que se puedan predecir con
este tipo de sistemas con el fin de encontrar distintos modos de terapias que sean más
eficaces a la hora de tratar estas enfermedades (u otras que sean susceptibles a la mo-

41



42 6.2. Posibles elementos a incorporar en los modelos.

delización matemática), ahorrando recursos y principalmente ahorrando tiempo.

El desarrollo de este trabajo ha sido incorporar una ecuación nueva en cada caṕıtulo, dando
a entender que estos modelos se pueden detallar hasta niveles insospechados con la inclusión
de nuevas ecuaciones y nuevos términos de interacción relacionados con los distintos elemen-
tos que se incorporen.

Obivamente, cuanto más detalle se implemente en los sistemas, mas recursos de análisis
numérico necesitaremos para estudiar estos modelos en función de los parámetros de estudio,
por lo tanto, un desarrollo en este tipo de modelos tiene que ir de la mano del desarrollo de
nuevas técnicas de simulación que sean más eficientes a la hora de estudiarlos.

6.2. Posibles elementos a incorporar en los modelos.

En este trabajo nos hemos restringido al estudio de modelos de ecuaciones diferenciales or-
dinarias con una población celular que hemos asumido que está homogéneamente dispuesta.
Estas aproximaciones parecen apropiadas para el estudio de las interacciones de los sistemas
cáncer-inmune en el contexto en que lo tratamos como una enfermedad sistémica que no está
aislada en un único tumor, aunque como hemos visto no estamos interesados de momento
en la distribución espacial de este tipo de células.

Sin embargo, podŕıamos introducir en estos modelos (de hecho son más que adecuados para
ello) factores que nos permitan modelizar tumores aislados. En ellos habŕıa que introducir
una estructura espacial del tumor lo cual, muchos autores han hecho de múltiples formas
distintas.

6.2.1. Discretización y modelos estocásticos.

Cualquier modelo de ODEs revisado en los anteriores caṕıtulos puede ser re-descrito como
un modelo de Monte Carlo, que trabaja con números enteros de células. Aqúı entra en juego
el estudio de los procesos estocásticos (por ejemplo, muy usado el Algoritmo de Gillespie en
estos modelos).

Sin embargo estos modelos son sobre todo utilizados para grandes poblaciones de células.
Esto se debe a que el comportamiento en estos casos converge a los sistemas de ecuaciónes
homólogos. No son tan utilizados en microtumores o densidades de poblaciones bajas ya que
este tipo de procesos tiene un comportamiento distinto en estos casos.

Cabe destacar la importancia de estos procesos a la hora de introducir factores que modelicen
las posibles mutaciones celulares debido al componente de aleatoriedad que tienen estos
procesos. Este también es un aspecto importante en el terreno de la oncoloǵıa ya que un
mismo cáncer puede mutar a otro tipo y su comportamiento tendŕıa que ser modelizado de
otra forma.

42



6. Discusión sobre los modelos de ODEs estudiados y posibles enfoques futuros.43

6.2.2. Agentes individuales y microambiente tumoral.

Los sistemas inmunes reales son extremadamente complejos en donde se incluyen gran canti-
dad de tipos celulares distintos aśı como diferentes tipos de moléculas de señalización. Estos
modelos en los que incorporamos gran cantidad de detalle solo son abarcables a través de
técnicas de simulación muy avanzadas que se escapan a este trabajo.

En la mayoŕıa de estos modelos lo que se hace es recrear el microambiente que rodea a
las células malignas añadiendo todo tipo de células y moléculas aśı como los nutrientes los
cuales juegan un papel imprescindible a la hora del crecimiento o recesión de las células
canceŕıgenas (es lo que proporciona enerǵıa para su desarrollo). En estos casos se realizan
gran cantidad de simulaciones para los parámetros que se desean controlar y con ello, ser
capaz de probar ciertas hipótesis concretas.

6.2.3. Estructura espacial.

Como hemos visto en estos modelos no hemos incorporado una estructura espacial para
estudiar los tumores. Una forma de solucionar este problema y seguir con este trabajo seŕıa
la incorporación de ciertas ecuaciones en derivadas parciales relativas a los problemas con
las cuales, a través de diversas técnicas anaĺıticas, podŕıamos estudiar la distribución celular.
Las técnicas que podŕıan usarse en dichos casos podŕıan ir desde el estudio de los posibles
tipos de solución que pueden presentar estas ecuaciones hasta la estimación del crecimiento
tumoral a través de ecuaciones de ondas y transporte.

43



CAPÍTULO 7.

Anexo

Algunos códigos escritos en MATLAB para generar algunas de las gráficas presentes en el
trabajo son:
Krischner - Panetta (datos)

1 % Datos t r e s e cuac ione s
2

3 g l o b a l p1 p2 g1 g2 g3 a b c mu2 mu3 r2 s1 s2
4 f = @tumor growth KP ; i n t e r v a l o = [ 0 3 0 0 ] ;
5 x0 = [10 100000 0 ] ;
6

7 %D i s t i n t o s v a l o r e s para l o s p a r m e t r o s
8

9 p1 = 0 . 1 2 4 5 ; p2 = 5 ; g1 = 2∗(10ˆ7) ; g2 = 10ˆ5 ; g3 = 10ˆ3 ; a = 1 ; b
= (10ˆ−9) ;

10 r2 = 0 . 1 8 ; mu2 = 0 . 0 3 ; mu3 = 10 ;
11

12 %c = 0 . 0 0 2 ;
13 c = 0 . 0 5 ;
14 %c = 0 . 0 2 ;
15

16 %s1 = r2∗g2∗mu2/a ;
17 s1 = 0 ;
18 s2 = 0 ;
19 %s2 = mu2∗mu3∗g1 /(p1−mu2) ;

Krischner - Panetta (solución)

1

2

3 datos KP ;
4 [ t , x]=ode45 ( @tumor growth KP , i n t e r v a l o , x0 ) ;
5

6 f i g u r e (1 )
7 p lo t ( t , x ( : , 1 ) , ’ r ’ , t , x ( : , 2 ) , ’ b ’ , t , x ( : , 3 ) , ’ g ’ )
8 l egend ( ’ C l u l a s tumorales ’ , ’ C l u l a s inmunes ’ , ’ Cytokinas ’ )
9 f i g u r e (2 )

10 subplot ( 3 , 1 , 1 )
11 p lo t ( t , x ( : , 1 ) , ’ r ’ )
12 t i t l e ( ’ C LULAS TUMORALES’ )
13 subplot ( 3 , 1 , 2 )

i



ii

14 p lo t ( t , x ( : , 2 ) , ’b ’ )
15 t i t l e ( ’ C LULAS INMUNES ’ )
16 subplot ( 3 , 1 , 3 )
17 p lo t ( t , x ( : , 3 ) , ’ g ’ )
18 t i t l e ( ’CYTOKINAS ’ )

Krischner - Panetta (ecuaciones)

1

2 f unc t i on f = tumor growth KP ( t , x )
3

4 \% Modelo de t r e s ecuac ione s
5

6 g l o b a l p1 p2 g1 g2 g3 a b c mu2 mu3 r2 s1 s2
7

8 f 1 = x (1) ∗ r2∗(1−b∗x (1 ) ) − ( ( a∗x (1 ) .∗ x (2 ) ) /( g2 + x (1) ) ) ;
9 f 2 = s1 + c∗x (1 ) + ( ( p1∗x (2 ) .∗ x (3 ) ) /( g1 + x (3) ) ) − mu2∗x (2 ) ;

10 f 3 = s2 + ( ( p2∗x (1 ) .∗ x (2 ) ) /( g3 + x (1) ) ) − mu3∗x (3 ) ;
11

12 f = [ f 1 ; f 2 ; f 3 ] ;

Campo vectorial

1

2 f unc t i on v e c t o r f i e l d ( deqns , xval , yval , t )
3 i f narg in ==3;
4 t =0;
5 end
6 m = length ( xval ) ;
7 n = length ( yval ) ;
8 x1 = ze ro s (n ,m) ;
9 y1 = ze ro s (n ,m) ;

10 f o r a=1:n
11 f o r b = 1 : n
12 pts = f e v a l ( deqns , t , [ xval ( a ) ; yval (b) ] ) ;
13 x1 (b , a ) = pts (1 ) ;
14 y1 (b , a ) = pts (2 ) ;
15 end
16 end
17

18 arrow=s q r t ( x1 .ˆ2+y1 . ˆ 2 ) ;
19 qu iver ( xval , yval , x1 . / arrow , y1 . / arrow , 0 . 5 , ’ r ’ ) ;
20 a x i s t i g h t ;
21 end

Diagrama de fases

1

2 \%Diagrama de f a s e s

ii



7. Anexo iii

3

4 \%d = 2 . 0 ; beta = 4∗(10ˆ−3) ; p = 1 . 1 3 1 ; g = 2 0 . 1 9 ; m = 0 .00311 ; a
= 1 . 6 3 6 ; n = 1 ; s = 0 . 3 1 8 ;

5 c l e a r
6 sys = i n l i n e ( ’ [ 1 . 6 3 6∗ x (1 ) ∗(1−4∗0.001∗x (1 ) ) −x (1 ) .∗ x (2 ) ;

1 .131 ./ (20 .19+ x (1) ) −0.00311∗x (1 ) .∗ x (2 ) − 2∗x (2 ) + 0.318 ] ’ , ’ t
’ , ’ x ’ )

7 v e c t o r f i e l d ( sys , 0 : 1 : 2 0 , 0 : 1 : 2 0 )
8 hold on
9 f o r x0 = 4 5 : 1 : 5 5

10 f o r y0 = 0 : 1 : 1 0
11 [ ts , xs ] = ode45 ( sys , [ 0 1 0 ] , [ x0 y0 ] ) ;
12 p lo t ( xs ( : , 1 ) , xs ( : , 2 ) )
13 end
14 end
15

16 f o r x0 = 4 5 : 1 : 5 5
17 f o r y0 = 0 : 1 : 1 0
18 [ ts , xs ] = ode45 ( sys , [ 0 1 0 ] , [ x0 y0 ] ) ;
19 p lo t ( xs ( : , 1 ) , xs ( : , 2 ) )
20 end
21 end
22

23 hold o f f
24 a x i s ( [ 0 6 0 6 ] )
25 f s i z e =15;
26 s e t ( gca , ’ x t i c k ’ , [ 0 : 1 : 6 ] , ’ FontSize ’ , f s i z e )
27 s e t ( gca , ’ y t i c k ’ , [ 0 : 1 : 6 ] , ’ FontSize ’ , f s i z e )
28 x l a b e l ( ’ x ( t ) ’ , ’ FontSize ’ , f s i z e )
29 y l a b e l ( ’ y ( t ) ’ , ’ FontSize ’ , f s i z e )
30 hold o f f

No se incluyen todos los códigos pues muchos son similares. Simplemente es un breve resumen
de la linea que se ha seguido en el trabajo.

iii



Bibliograf́ıa

[1] Interactions Between the Immune System and Cancer: A Brief Review of Non-spatial
Mathematical Models. Eftimie, R., Bramson, J.L. and Earn, D.J.D. Bull. Math. Biol.
(2011) 73: 2

[2] Araujo, R., McElwain, D., 2004. A history of the study of solid tumor growth: the
contribution of mathematical modeling. Bull. Math. Biol. 66, 1039–1091.

[3] Well posedness of an integrodifferential kinetic model of Fokker-Planck type for angio-
genesis, A. Carpio, G. Duro, Nonlinear Analysis-Real World Applications 30, 184-212,
2016

[4] Modelizacion y simulacion en sistemas dinámicos, A. Carpio, 2017, eprints.ucm.es

[5] Bellomo, N., Delitala, M., 2008. From the mathematical kinetic, and stochastic game
theory to modeling mutations, onset, progression and immune competition of cancer
cells. Phys. Life Rev. 5, 183–206.

[6] von Bertalanffy, L., 1957. Quantitative laws in metabolism and growth. Q. Rev. Biol.
32, 217–231.

[7] Hart, D., Shochat, E., Agur, Z., 1998. The growth law of primary breast cancer as
inferred from mammography screening trials data. Br. J. Cancer 78, 382–387.

[8] Laird, A., 1964. Dynamics of tumor growth. Br. J. Cancer 18, 490–502.

[9] Bifurcation Analysis of Tumor-Immune ODE System L.G. de Pillis and A.E. Raduns-
kaya; August 15, 2002

[10] Norton, L., 1988. A Gompertzian model of human breast cancer growth. Cancer Res.
48, 7067–7071.

[11] d’Onofrio, A., 2008. Metamodeling tumor-immune system interaction, tumor evasion
and immunotherapy. Math. Comput. Model. 47, 614–637.

[12] d’Onofrio, A., 2005. A general framework for modeling tumor-immune system compe-
tition and immunotherapy: mathematical analysis and biomedical inferences. Physica D
208, 220–235.

[13] Kuznetsov, V., Makalkin, I., Taylor, M., Perelson, A., 1994. Nonlinear dynamics of
immunogenic tumors: parameter estimation and global bifurcation analysis. Bull. Math.
Biol. 2(56), 295–321.

[14] Sotolongo-Costa, O., Molina, L.M., Perez, D.R., Antoraz, J., Reyes, M.C., 2003. Beha-
vior of tumors under nonstationary therapy. Physica D 178, 242–253.

44



BIBLIOGRAFÍA 45

[15] Moore, H., Li, N., 2004. A mathematical model for chronic myelogenous leukemia
(CML) and T cell interaction. J. Theor. Biol. 227, 513–523.

[16] de Pillis, L., Radunskaya, A., 2003a. The dynamics of an optimally controlled tumor
model: a case study. Math. Comput. Model. 37, 1221–1244.

[17] Kirschner, D., Panetta, J., 1998. Modeling immunotherapy of the tumor-immune in-
teraction. J. Math. Biol. 37, 235–252.

[18] Owen, M., Sherratt, J., 1998. Modeling the macrophage invasion of tumors: effects on
growth and composition. Math. Med. Biol. 15, 165–185.

45


	Interacciones entre el Sistema Inmune y Cáncer. Modelos matemáticos sin dependencia espacial.
	Introducción prin.

	Teoría del estudio cualitativo de los sistemas dinámicos.
	Diagrama de fases.
	Puntos críticos.
	Bifurcaciones ana.
	Bifurcaciones tipo nodo-silla.
	Bifurcaciones tipo transcrítico.
	Bifurcación tipo Pitchfork.
	Bifurcaciones tipo Hopf.


	Modelo de una ecuación: Crecimiento tumoral.
	Modelo de Gompertz.
	Descripción del modelo de una ecuación diferencial.
	Ejemplo del modelo logístico.


	Modelos de ecuaciones: Interacciones entre células tumorales y sistema inmune.
	Desarrollo del modelo.
	Modelo teórico de Kuznestov (1994).
	Estados estacionarios.
	Bifurcaciones en el sistema.

	Simulación numérica del modelo de Kurznestov.
	Eliminación del tumor sin recidivas.
	Crecimiento del tumor.
	Crecimiento y decrecimiento tumoral periódico.

	Modelo teórico de Sotolongo-Costa (2003).
	Estados estacionarios.
	Dinámica y análisis de la estabilidad sot.
	Interpretación biológica sot.

	Simulación numérica del modelo de Sotolongo-Costa.
	Descontrol del sistema inmune.
	Control de las células malignas.
	Crecimiento y decrecimiento peródico significativo.


	Modelo de tres ecuaciones. Interacciones entre células cancerígenas y dos componentes del ambiente tumoral (S.I)
	Crecimiento tumoral controlado por dos tipos de células efectoras.
	Simulación numérica del modelo Moore-Li.

	Crecimiento tumoral controlado por células efectoras y células normales.
	Crecimiento tumoral controlado por células efectoras y citoquinas.
	Simulación numérica del modelo Krischner - Panneta.


	Discusión sobre los modelos de ODEs estudiados y posibles enfoques futuros.
	Incremento de la complejidad de los sistemas.
	Posibles elementos a incorporar en los modelos.
	Discretización y modelos estocásticos.
	Agentes individuales y microambiente tumoral.
	Estructura espacial.


	Anexo