Facultad de Ciencias

Trabajo Fin de Grado

Grado en Matemáticas

Las matemáticas de la fisiología del corazón

Autora: Elena de la Vega de la Peña

Tutora: Begoña Cano Urdiales

Índice general

Índice general I

Índice de figuras II

Introducción 1

1. Presentación del modelo a imitar 5

1.1. Fisiologı́a del corazón . . . . . . . . . . . . . . . . . . . . . . . 5
1.2. Cualidades a imitar en el modelo matemático y variables clave
del mismo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2. Sucesivas aproximaciones al modelo de Zeeman 10

Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Ejemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Ejemplo 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Ejemplo 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Ejemplo 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Nuestra propuesta como modificación del modelo de Zeeman . 28

3. Teoremas de existencia o no de soluciones periódicas en el

plano 32
3.1. Aplicación de los teoremas a los ejemplos del capı́tulo 2 . . . . 40
3.2. Aplicación de los teoremas a la ecuación de Van der Pol . . . . 42

4. Simulación de un electrocardiograma 44

A. Programas de MATLAB 51

Bibliografı́a 56

I
Índice de figuras

1.1. Fisiologı́a del corazón . . . . . . . . . . . . . . . . . . . . . . . 6

1.2. Sistema circulatorio . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1. Ejemplo 1 para  = 1/10 . . . . . . . . . . . . . . . . . . . . . 11

2.2. Ejemplo 1 para  = 1/100 . . . . . . . . . . . . . . . . . . . . 12
2.3. Ejemplo 2 para  = 1/10 . . . . . . . . . . . . . . . . . . . . . 13
2.4. Ejemplo 2 para  = 1/100 . . . . . . . . . . . . . . . . . . . . 14
2.5. Ejemplo 3 para  = 1/10 . . . . . . . . . . . . . . . . . . . . . 15
2.6. Ejemplo 3 para  = 1/100 . . . . . . . . . . . . . . . . . . . . 16
2.7. Ejemplo 4 para  = 1/10 . . . . . . . . . . . . . . . . . . . . . 17
2.8. Ejemplo 4 para  = 1/100 . . . . . . . . . . . . . . . . . . . . 18
2.9. Ejemplo 5 para  = 1/5 . . . . . . . . . . . . . . . . . . . . . . 19
2.10. Ejemplo 5 para  = 1/100 . . . . . . . . . . . . . . . . . . . . 20
2.11. Ejemplo 6 para  = 1/5 . . . . . . . . . . . . . . . . . . . . . . 22
2.12. Ejemplo del sistema (2.1) con b0 = 0.1 . . . . . . . . . . . . . 23
2.13. Ejemplo del sistema (2.1) con b0 = 1 . . . . . . . . . . . . . . 24
2.14. Ejemplo del sistema (2.1) con b0 = −1 . . . . . . . . . . . . . 25
2.15. Ejemplo del sistema (2.2) con a = 0 . . . . . . . . . . . . . . 26
2.16. Ejemplo del sistema (2.2) con a = 0.2 . . . . . . . . . . . . . 27
2.17. Ejemplo del sistema (2.2) con a = 2 . . . . . . . . . . . . . . 28
2.18. Ejemplo del sistema (2.3) con x0 = 0.5 y a = 1 . . . . . . . . . 29
2.19. Ejemplo del sistema (2.3) con x0 = 0.5 y a = 2 . . . . . . . . 30
2.20. Ejemplo del sistema (2.3) con x0 = 0.5 y a = 0.5 . . . . . . . 31

4.1. Componentes de un electrocardiograma normal . . . . . . . . 44

4.2. Esquema de los tres osciladores en el caso normal . . . . . . . 45
4.3. Electrocardiograma generado a partir del sistema (2.3) . . . . 46
4.4. Electrocardiograma normal . . . . . . . . . . . . . . . . . . . . 48
4.5. Electrocardiograma aleteo ventricular . . . . . . . . . . . . . . 49

II
Introducción

En la actualidad, los modelos y simulaciones matemáticas se están convir-

tiendo en una herramienta cada vez más importante para el estudio de pro-
blemas de la vida real. Concretamente, si hablamos del campo de la medicina,
el desarrollo de estas técnicas nos permite comprender mejor los mecanismos
fisiológicos e incluso los comportamientos patológicos, permitiéndonos actuar
sobre ellos y predecir las posibles respuestas.

El objetivo principal de este Trabajo de Fin de Grado es proponer y

justificar matemáticamente un modelado de los movimientos del corazón.
Para ello, estudiamos la parte dinámica del comportamiento del corazón y
consideramos un sistema de ecuaciones diferenciales ordinarias en R2 para
representar dicho comportamiento. También se hace referencia a la relación
existente entre el modelo considerado y la ecuación de Van der Pol, puesto
que ambos presentan ciertas caracterı́sticas comunes, como la existencia de
comportamientos periódicos. En dicho modelado se ven reflejados los equi-
librios de sı́stole y diástole, ası́ como un determinado umbral que provoca
que se desencadene la onda electroquı́mica que da lugar a la contracción del
corazón. Para ello, nuestro objetivo inicial fue estudiar el modelo propuesto
por E. C. Zeeman [7], matemático británico nacido en 1925 y conocido por
sus trabajos sobre topologı́a geométrica, teorı́a del caos y teorı́a de singu-
laridades. Sin embargo, tras estudiar su propuesta final de modelo en dos
dimensiones, consideramos también una pequeña variación que nos parece
más acertada, y que describimos en el trabajo. En el mismo, se presentan
también ciertos teoremas y gráficas de mapas de fase que justifican e ilustran
cada modelo estudiado.

A lo largo del siguiente trabajo se usan resultados y conceptos relaciona-

dos con los sistemas de ecuaciones diferenciales ordinarias y con el modelado
matemático adquiridos principalmente en las asignaturas de “Ecuaciones Di-
ferenciales” y “Matemáticas Aplicadas a las Ciencias Naturales y Sociales”.
También conceptos básicos de álgebra matricial estudiados en la asignatura

1
“Álgebra y geometrı́a lineales I”.

En el primer capı́tulo, se exponen conceptos básicos de la fisiologı́a del

corazón y del papel que desarrolla el mismo en el sistema circulatorio. Estas
nociones van encaminadas a conseguir más adelante que el modelo matemáti-
co propuesto se ajuste todo lo posible a la realidad. En la segunda sección de
este capı́tulo se presentan las cualidades y variables que Zeeman considera
necesarias para su modelo. Además, se introducen ciertos fenómenos que son
provocados por posibles alteraciones en el funcionamiento del corazón, como
la ley de Starling o el experimento de Rybak. En el capı́tulo siguiente se
prueba cómo el modelo final es capaz de describir estos fenómenos.

En el segundo capı́tulo, se realiza un estudio de los distintos modelos en

términos de dos variables que aparecen en el artı́culo de Zeeman [7], hasta
llegar al que él considera más adecuado. Para ello estudiamos los correspon-
dientes sistemas de ecuaciones diferenciales en R2 y el comportamiento de sus
respectivos puntos de equilibrio utilizando el teorema de Hartman-Grobman.
Además, se acompaña cada modelo con su correspondiente mapa de fase que
ha sido generado con MATLAB y que refleja dicho comportamiento. Al final
de este capı́tulo se expone y justifica la modificación del modelo de Zeeman
que nosotras consideramos más adecuado para representar el movimiento del
corazón. Concluimos el capı́tulo realizando un estudio de cómo cambia el
comportamiento del modelo si hacemos variar el valor de la tensión del co-
razón. Esto justifica el hecho de que dicho modelo satisfaga los fenómenos
presentados en el capı́tulo anterior.

En el tercer capı́tulo, presentamos una serie de teoremas que justifican o

niegan la existencia de órbitas cerradas en un sistema de ecuaciones diferen-
ciales. Entre ellos exponemos el teorema generalizado de Poincaré-Bendixson,
un ejemplo que ilustre dicho teorema y el teorema de Liénard junto con su
demostración. Para ello seguimos la demostración que aparece en [3] y de-
tallamos algunos pasos. En la primera sección de este capı́tulo aplicamos los
teoremas presentados anteriormente a nuestro modelo, con el fin de demostrar
la existencia de una solución periódica y con ello probar que la modificación
del modelo de Zeeman considerada es capaz de representar la periodicidad
del latido del corazón. En la segunda sección introducimos la ecuación de
Van der Pol y su relación con el modelo que estamos considerando.

Aunque nuestro objetivo inicial fue únicamente estudiar el conocido mo-

delo de Zeeman, encontramos en la literatura posteriormente el modelo de
Gois [1] que imita mejor el comportamiento real del electrocardiograma. Por

2
ello, hemos añadido el capı́tulo 4. Éste comienza con una breve descripción de
lo que es un electrocardiograma y cómo se transmite el impulso eléctrico que
genera la contracción y la relajación del corazón a través del mismo. Además,
se presenta el nuevo modelo propuesto por Gois. Finalmente, comparamos
el electrocardiograma generado por el modelo modificado, modelo con el que
hemos trabajado en capı́tulos anteriores, con el electrocardiograma que se
genera con el modelo que propone Gois.

Por último, me gustarı́a agradecer a mi tutora, Begoña Cano, su tiempo

y su dedicación a lo largo de la elaboración de este trabajo.

3
4
Capı́tulo 1

Presentación del modelo a

imitar

1.1. Fisiologı́a del corazón

El corazón es un órgano muscular encargado de impulsar regularmente
la sangre por todo el cuerpo hasta el resto de órganos y tejidos. Junto con
los vasos sanguı́neos y la propia sangre constituyen el sistema circulatorio,
responsable del aporte de nutrientes y oxı́geno a los tejidos y de la recogida
de dióxido de carbono y sustancias de deshecho del metabolismo.
El corazón está formado por cuatro cavidades huecas (dos aurı́culas en la
parte superior y dos ventrı́culos en la parte inferior) y cuatro válvulas que
impiden que la sangre fluya en sentido contrario, como se aprecia en la figura
1.1. Cada aurı́cula está conectada con el ventrı́culo del mismo lado a través de
un orificio. La separación en dos cavidades en cada lado es necesaria, puesto
que el corazón sólo tiene capacidad para expulsar la sangre, no para hacerla
llegar a su interior. Esto se debe a que el corazón está hecho de un tejido
no-rı́gido. Para poder realizar una buena contracción, el corazón debe tener
el ventrı́culo totalmente lleno. Es por esto que existe la aurı́cula, la cual se
contrae ligeramente para llenar el ventrı́culo [2].

Además, el interior del corazón está divido en dos mitades separadas por
un tabique, diferenciando ası́ la parte derecha que recibe la sangre procedente
del cuerpo y la envı́a a la circulación pulmonar y la parte izquierda que recibe
la circulación de los pulmones y la distribuye a todo el organismo.
Por tanto, gracias al trabajo del corazón y a un gradiente de presiones,
la sangre sigue dos circuitos sanguı́neos. Uno de ellos es de baja presión (pa-
ra no dañar las delicadas membranas de los pulmones) e impulsa la sangre

5
 Figura 1.1: Fisiologı́a del corazón

pobre en oxı́geno desde el ventrı́culo derecho a través de la arteria pulmonar

a los pulmones, donde posteriormente se realiza el intercambio de gases y
se oxigena la sangre. El otro circuito es de mayor presión, puesto que debe
hacer llegar la sangre a todo el cuerpo y retirar a su vez las sustancias de
deshecho. Es el encargado de transportar la sangre desoxigenada desde los
distintos órganos del cuerpo al corazón. Se pueden ver estos dos circuitos en
la figura 1.2.

Durante el ciclo cardı́aco existen dos posibles estados de equilibrio. Uno

corresponde al momento de mayor contracción del corazón, llamado sı́stole,
mientras que el otro se corresponde con el estado de relajación o diástole.
Esta contracción se debe a una onda electroquı́mica que, en condiciones
normales, se desencadena de manera espontánea en el nodo sinoauricular
[6] (formado por células musculares cardı́acas especializadas que inician y
propagan el impulso cardı́aco). Ésta se expande primero por las aurı́culas,
provocando que las fibras musculares se contraigan y empujen la sangre a

6
 Figura 1.2: Sistema circulatorio

los ventrı́culos; y posteriormente se expande rápidamente por los ventrı́culos

provocando la contracción total e impulsando ası́ la sangre fuera del corazón.
Inmediatamente después las fibras musculares se relajan, devolviendo el co-
razón al estado de diástole.

1.2. Cualidades a imitar en el modelo ma-

temático y variables clave del mismo
Describiremos en este trabajo la parte básica del modelo que Zeeman
presenta en [7]. La principal caracterı́stica de este modelo es que se basa
en simular la parte dinámica del corazón (relativamente simple), más que
la bioquı́mica (la cual es más complicada). Debido a la fisiologı́a cardı́aca
descrita anteriormente, Zeeman busca que el modelo matemático cumpla las
siguientes cualidades dinámicas:

(1) Tener un equilibrio estable.

7
(2) Tener un umbral, que desencadene una acción, es decir, la contracción
del corazón.

(3) Rápida vuelta al equilibrio.

Zeeman considera el estado de relajación del corazón, diástole, como el

equilibrio estable de la cualidad 1. Cuando la onda electroquı́mica que pro-
voca la contracción del corazón en sı́stole alcanza cada fibra muscular, se
corresponde con la acción requerida en la cualidad 2. Cada fibra se man-
tiene contraı́da 1/5 segundos aproximadamente y después vuelve a relajarse
rápidamente, es decir, cumple la cualidad 3. Por tanto, estas tres cualidades
describen el comportamiento local de una fibra muscular del corazón.

De esta manera, Zeeman toma dos variables para representar su modelo

matemático. Considera la variable x como la longitud de una fibra muscular
cardı́aca y la variable b que representa una cierta variable de control electro-
quı́mico, correspondiente a cierta concentración o flujo de una determinada
sustancia quı́mica. Puesto que la tensión en la fibra generada por la presión
sanguı́nea también juega un papel importante, se considera el parámetro a
correspondiente a dicha tensión.

El modelo matemático que se obtiene al final es capaz de describir tam-

bién ciertos fenómenos que se producen debido a posibles cambios, como por
ejemplo, la capacidad de adaptarse a volúmenes cambiantes de sangre. El
ejemplo anterior es el principio o ley de Frank-Starling [7]: Cuanto mayor
es el volumen de llenado al final de la diástole, mayor es la fuerza de con-
tracción y por tanto, mayor es el volumen de sangre impulsado en la sı́stole.
Cuando aumenta el volumen diastólico final, las fibras cardı́acas en reposo se
alargan, lo que provoca que con la llegada del impulso cardı́aco la velocidad
de contracción aumente y con ella la fuerza de contracción ejercida. Por el
contrario, si el volumen diastólico es menor, la fuerza de contracción también
será menor.
Este sistema responde excelentemente ante situaciones extremas o de
emergencia. Supongamos que por una determinada causa sentimos miedo.
Esto genera un aumento de adrenalina en el torrente sanguı́neo. Esta adre-
nalina provoca una contracción de las arterias y con ello un incremento del
ritmo cardı́aco. Todo esto da lugar a un aumento de la presión sanguı́nea,
por lo que las aurı́culas impulsarán más sangre en los ventrı́culos. La ley de
Starling describe cómo esta extensión de los músculos del corazón da lugar
a una contracción más potente, lo que contrarresta el aumento de presión y
hace que la sangre circule más rápidamente.

8
 Pero este mecanismo tiene un lı́mite, puesto que a partir de cierto punto
de distensión de los ventrı́culos (como puede ocurrir por ejemplo cuando una
persona con presión sanguı́nea alta recibe una fuerte impresión), la fuerza
de contracción ya no se ve afectada, pudiendo dar lugar incluso a un fallo
cardı́aco.

Otro ejemplo de estas posibles alteraciones es el experimento que realiza

Boris Rybak en 1957 con una rana. Si el corazón es extraı́do del cuerpo, deja
de latir. Sin embargo, si se le aplica una leve tensión, vuelve a latir y continua
ası́ durante algunas horas. Si se deja de aplicar dicha tensión, el corazón se
para. Esto quiere decir que, sin tensión, el umbral que desencadena la acción
desaparecerı́a y por tanto el latido cesarı́a.
Algo parecido ocurre cuando el corazón solo hace de ’bypass’ en una ope-
ración, de forma que desaparece la presión sanguı́nea en el mismo, el latido
se vuelve muy lento y el corazón no cambia prácticamente de tamaño.

En el próximo capı́tulo se puede ver cómo los modelos presentados satis-

facen estos fenómenos.

9
Capı́tulo 2

Sucesivas aproximaciones al
modelo de Zeeman

En este capı́tulo haremos un estudio de los ejemplos que aparecen en el

artı́culo de Zeeman [7] hasta llegar a un modelo con las cualidades reque-
ridas que han sido descritas en el capı́tulo anterior. Representaremos estas
cualidades usando un sistema de ecuaciones diferenciales en R2 , donde las
coordenadas serán x la longitud de una fibra muscular del corazón y b una
variable de control electroquı́mico.

Debido a la cualidad 1, el sistema ha de tener un punto de equilibrio que

inicialmente consideraremos (0, 0) ∈ R2 .

Ejemplo 1
Zeeman propone empezar con un caso lineal, donde las ecuaciones son
(
ḃ = −b,
ẋ = −x

con  una constante pequeña positiva, para que ası́ el módulo de uno de
los autovalores sea más grande comparado con el del otro autovalor y den
lugar a que las órbitas no demasiado cercanas a x = 0 tengan un flujo rápido.

La matriz asociada al sistema anterior es

−1 0
!
A= 1 .
0 −


10
 1
Los autovalores correspondientes son λ1 = − y λ2 = −1. Como ambos

son reales y negativos, el único punto de equilibrio (0, 0) es un nodo estable.
 
0
Asociados a estos autovalores tenemos los siguientes autovectores v1 =
  1
1
y v2 = respectivamente. Puesto que λ1 < λ2 , el vector dominante
0
será v2 cuanto t crece.

Figura 2.1: Ejemplo 1 para  = 1/10

La solución del sistema es

b(t) = c2 e−t ,


x(t) = c1 e−t/ .
Las trayectorias están dibujadas en las figuras 2.1 y 2.2 utilizando un
programa en MATLAB que dibuja el mapa de fases correspondiente a un
sistema diferencial. Se observa que dichas trayectorias son tangentes a la rec-
ta correspondiente al autovector lento v2 (asociado al autovalor mayor, que
llamaremos autovalor lento) cuando se acercan al origen; y son casi paralelas

11
 Figura 2.2: Ejemplo 1 para  = 1/100

a la recta generada por el autovector rápido v1 (asociado al autovalor menor,

o autovalor rápido) cuando nos alejamos del origen. Es decir, cuanto más
alejados estemos de la recta x = 0, más verticales van a ser las trayectorias,
ya que la variable x decrece mucho más rápido que b y por tanto es la compo-
nente dominante. Puesto que x decrece más rápido denotaremos por ecuación
rápida a ẋ = −x y ecuación lenta a ḃ = −b. En la figura 2.2, se tiene que 
es más pequeño y por tanto, el decaimiento de x es más pronunciado.

Ejemplo 2
En segundo lugar Zeeman considera el siguiente sistema
(
ḃ = −b,
ẋ = −(x + b)
cuya matriz asociada es

12
 −1 0
!
A= 1 1 .
− −
 
El único equilibrio sigue siendo el punto (0, 0). En este caso los autovalores
1
siguen siendo λ1 = − y λ2 = −1 pero sus respectivos autovectores son
    
0 −1 + 
v1 = y v2 = .
1 1

Figura 2.3: Ejemplo 2 para  = 1/10

Puesto que de nuevo λ1 < λ2 < 0 el equilibrio (0, 0) sigue siendo un nodo
estable y el vector dominante es v2 cuando t → ∞, como podemos compro-
bar en la figura 2.3. A lo largo de la recta x + b = 0 las trayectorias pasan
con pendiente horizontal porque ẋ = 0. Donde esto ocurre, lo llamaremos
variedad lenta, puesto que ocurre que la ecuación rápida da variación cero en
esos puntos. La dibujaremos en rojo en todas las figuras, al igual que hicimos
en las figuras 2.1 y 2.2. La variedad invariante determinada por λ2 = −1 es
(1 − )x + b = 0, para la cual la variedad lenta x + b = 0 es una aproximación,

13
como puede observarse claramente en la figura 2.3.

Figura 2.4: Ejemplo 2 para  = 1/100

La solución del sistema es

b(t) = c2 (−1 + )e−t ,


x(t) = c1 e−t/ + c2 e−t .

Como en el caso anterior las trayectorias son tangentes a la recta generada
por el autovector lento v2 cuando t → ∞ (es decir, al acercarse al (0, 0)) y
aún más paralelas a v1 cuando t → −∞ (es decir, al alejarse del (0, 0)).
Además, cuanto más pequeño es , más cerca están la variedad lenta y la
variedad invariante asociada a λ2 = −1 y más verticales son las trayectorias
antes de acercarse a la variedad lenta, por donde prácticamente caen al origen
(0, 0). (Ver figuras 2.3 y 2.4).

Ejemplo 3
A continuación Zeeman introduce esta modificación del ejemplo 2, donde
cambia la ecuación lenta

14
 (
ḃ = x,
ẋ = −(x + b).
La matriz asociada al sistema es
!
0 1
A= 1 1 .
− −
 

Figura 2.5: Ejemplo 3 para  = 1/10

En este caso la variedad lenta es la misma que en el ejemplo 2, ya que la

segunda ecuación es la misma. Además, sobre ella tenemos que −b = x y por
tanto la primera ecuación es la misma y el comportamiento será el mismo en
la variedad lenta.
√ √
−1 + 1 − 4 −1 − 1 − 4
Los autovalores son λ1 = y λ2 = , ambos
2 2

15
reales negativos y los autovectores asociados son
 √   √ 
−1 − 1 − 4 −1 + 1 − 4
v1 =  2  , v2 =  2 
1 1
respectivamente.

Figura 2.6: Ejemplo 3 para  = 1/100

De esta manera se tiene que el mapa de fases de la figura 2.5 es muy

parecido al del ejemplo 2 (figura 2.3), solo con pequeñas diferencias a la vista
respecto de la variación de la variable b. De hecho, cuando  es muy pequeño,
las diferencias entre ambos mapas de fase son imperceptibles (comparar fi-
guras 2.4 y 2.6). Más adelante veremos por qué era necesario hacer el cambio
en la primera ecuación para conseguir la cualidad 3.

Ejemplo 4
Zeeman considera entonces un ejemplo no lineal

16
 (
ḃ = x,
ẋ = −(x3 + x + b).
Ocurre que la matriz jacobiana de este sistema es
 
0 1
J(b, x) =  1 3x2 + 1  .
− −
 

Figura 2.7: Ejemplo 4 para  = 1/10

Evaluando J(b, x) en el único punto crı́tico (0, 0) del sistema se obtiene

!
0 1
J(0, 0) = 1 1 ,
− −
 
que es la misma matriz del ejemplo 3. Entonces la aproximación lineal de este
sistema en el origen se corresponde con el ejemplo 3, por lo que aplicando
el teorema de Hartman-Grobman, el origen será un nodo estable y como se
ve en las figuras 2.7 y 2.8 el comportamiento local será igual. Sin embargo,

17
la variedad lenta x3 − x + b = 0 es ahora una curva. Las trayectorias lejos
del origen se comportan como en el ejemplo 1, son casi verticales. Es más,
podemos decir que los ejemplos 1, 2, 3 y 4 son difeomorfos.

Figura 2.8: Ejemplo 4 para  = 1/100

Ejemplo 5
A continuación Zeeman considera el sistema
(
ḃ = x,
ẋ = −(x3 − x + b)
donde se ha realizado un cambio de signo en uno de los términos de la segunda
ecuación. De esta manera la matriz jacobiana es de la forma
 
0 1
J(b, x) =  1 3x2 − 1  .
− −
 

18
 Figura 2.9: Ejemplo 5 para  = 1/5

Como en los casos anteriores, el único punto de equilibrio es el (0, 0), pero
evaluando la matriz jabociana en ese punto nos queda
!
0 1
J(0, 0) = 1 1 ,
−
 
cuyos autovalores son √
1± 1 − 4
λ= ,
2
√
y como 1 − 4 < 1 para  pequeño, entonces 1 ± 1 − 4 > 0 y tenemos
que los dos autovalores son positivos, √
dando lugar a un nodo inestable
√ en el
origen (0, 0). Llamaremos λ1 = (1 + 1 − 4)/2 y λ2 = (1 − 1 − 4)/2.
Debido al cambio de signo con respecto al ejemplo 4 tenemos que la variedad
lenta es ahora una curva cúbica con forma de S. Además, hemos perdido la
cualidad 1, puesto que ahora el origen en una fuente. Más concretamente, la
parte central de la variedad lenta ha pasado a ser repulsora. Los puntos de la

19
variedad lenta donde b es un máximo o mı́nimo frente a x vienen dados por
3x2 = 1. También ha aparecido una órbita cerrada a la que tienden el resto
de trayectorias, como se ve en las figuras 2.9 y 2.10. En el próximo capı́tulo
demostraremos la existencia de la misma aplicando el teorema de Liénard.
Los autovectores asociados a λ1 y λ2 son

Figura 2.10: Ejemplo 5 para  = 1/100

 √   √ 
1− 1 − 4 1 + 1 − 4
v1 =  2  , v2 =  2 .
1 1

En este caso, como 0 < λ2 < λ1 , el vector dominante es v2 cerca del

(0, 0) (cuando t → −∞). La pendiente de la recta correspondiente a v2 se
comporta como 1 − , ya que haciendo el desarrollo de Taylor de (1 − 4)1/2
nos queda
√
1 + 1 − 4 1 1 1 1
= + (1 − 4)1/2 ' + (1 − 2) = 1 − .
2 2 2 2 2

20
En cambio la pendiente de la curva x3 − x + b = 0 en (0, 0) es 1.

Como se puede ver en la Figura 2.9, en el mapa de fases las trayectorias

tienen pendiente horizontal a lo largo de la curva x3 − x + b = 0 y vertical en
x = 0. Por encima de la curva x3 −x+b = 0 las trayectorias van hacia abajo y
por debajo de la curva van hacia arriba. Por otro lado, en el semiplano x > 0,
apuntan hacia la derecha, mientras que en el semiplano x < 0 apuntan hacia
la izquierda.

Ejemplo 6
Por último Zeeman presenta el siguiente ejemplo
(
ḃ = x − x0 ,
ẋ = −(x3 − x + b),
√
considerando x0 > 1/ 3. La principal diferencia con el ejemplo 5 es que
ahora el punto de equilibrio está en la parte atractiva de la variedad lenta y
por tanto es atractor en lugar de repulsor. Dicho punto de equilibrio ya no
es el origen, es el punto (−x30 + x0 , x0 ). Denotaremos b0 = −x30 + x0 .

La matriz jacobiana correspondiente al sistema de ecuaciones anterior es

 
0 1
J(b, x) =  1 3x2 − 1  ,
− −
 
y evaluada en el punto de equilibrio nos queda
 
0 1
J(b0 , x0 ) =  1 3x2 − 1  .
− − 0
 
Los autovalores asociados a la matriz anterior son
p
−(3x20 − 1) ± (3x20 − 1)2 − 4
λ= .
2
√
Puesto que x0 > 1/ 3, se tiene que 3x20 − 1 > 0. Esto implica que
λ2 < λ1 < 0 y entonces el equilibrio (x0 , b0 ) es un nodo estable.
En este caso, los autovectores correspondientes son
 p   p 
−(3x20 − 1) − (3x20 − 1)2 − 4 −(3x20 − 1) + (3x20 − 1)2 − 4
v1 =  2
 , v2 = 
2
,
1 1

21
 Figura 2.11: Ejemplo 6 para  = 1/5

con v1 vector dominante cuando t → ∞. La inversa de su pendiente se com-

porta como −(3x20 − 1) + /(3x20 − 1), frente a la inversa de la pendiente de
la curva x3 − x + b = 0 en (b0 , x0 ) que es −(3x20 − 1).

Estudiando el mapa de fases las trayectorias tienen pendiente horizontal

a lo largo de la curva x3 − x + b = 0 y vertical en x = x0 . Por encima de
la curva x3 − x + b = 0 las trayectorias van hacia abajo y por debajo de la
curva van hacia arriba. Por otro lado, en el semiplano x > x0 , apuntan hacia
la derecha, mientras que en el semiplano x < x0 apuntan hacia la izquierda.

Si observamos el mapa de fases de la figura 2.11 (correspondiente a x0 =

0.75) podemos distinguir estas dos situaciones:
1. Si se sobrepasa un umbral cercano (pero mayor) a
1 1 2
b = √ − ( √ )3 = √ ,
3 3 3 3

22
 que es el valor donde la variedad lenta toma un máximo local de b con
respecto de x, las órbitas pasan a una zona de mayor contracción de x
(sı́stole) y vuelven después al equilibrio de diástole.

2. En cambio, si no se sobrepasa dicho umbral en la variable b con el

músculo sin contraer, las órbitas van directamente al equilibrio sin pasar
por la sı́stole.

Veamos ahora por qué Zeeman introdujo el cambio del ejemplo 2 al ejem-
plo 3.
Si hubiéramos tenido en el sistema ḃ = −(b − b0 ) como ecuación lenta, el
sistema habrı́a sido
(
ḃ = −(b − b0 ),
(2.1)
ẋ = −(x3 − x + b).

Figura 2.12: Ejemplo del sistema (2.1) con b0 = 0.1

23
 √
En este caso habrı́a tres equilibrios si |b0 | < 2/(3 3), puesto que la recta
b = b0 cortarı́a a la variedad lenta en tres puntos, siendo dos de ellos sumide-
ros y el otro una fuente. Para probar esto consideramos la matriz jacobiana
asociada al sistema anterior
 
−1 0
J(b, x) =  1 −3x2 + 1  .
−
 
Los autovalores asociados a dicha matriz son
p
−( − 1 + 3x2 ) ± ( − 1 + 3x2 )2 − 4(3x2 − 1)
λ= .
2
Simplificando las cuentas, obtenemos que los autovalores son
1 − 3x2
λ1 = −1, λ2 = ,

donde x es el valor de la ordenada de cada punto de equilibrio del sistema.

Figura 2.13: Ejemplo del sistema (2.1) con b0 = 1

24
 √ √
Se tiene que un autovalor es positivo y otro negativo si x ∈ (−1/ 3, 1/ 3),
por lo que el punto de equilibrio correspondiente
√ es√un punto de silla, es decir,
es inestable. En cambio, si x ∈ (∞, −1/ 3) ∪ (1/ 3, ∞) se tiene que λ2 < 0
y por tanto los puntos de equilibrio correspondientes son nodos estables. Esto
contradecirı́a la cualidad 3, ya que habrı́a trayectorias que se quedarı́an en
sı́stole, como de hecho puede observarse en la figura 2.12.
√
Si b0 > 2/(3 3), solo habrı́a un equilibrio (b0 , x0 ), que corresponderı́a a la
sı́stole. En este caso hay dos autovalores negativos, por el razonamiento ante-
rior, por lo que el equilibrio serı́a un nodo estable, como se ve en la figura 2.13.
√
Zeeman, sin embargo, no dice que si b0 < −2/(3 3), el mapa de fases no
diferirı́a mucho del ejemplo 6, habrı́a solo un equilibrio estable, diástole, al
que tenderı́an todas las órbitas, como se ve en la Figura 2.14.

Figura 2.14: Ejemplo del sistema (2.1) con b0 = −1

Aunque Zeeman propone el modelo correspondiente al ejemplo 6 inten-

25
tando imitar la cualidad 2, creemos que no es un buen modelo puesto que
en cualquier caso las soluciones van al nodo estable, sin producirse ninguna
oscilación, lo cual quiere decir que el corazón se pararı́a y provocarı́a un fallo
cardı́aco.

Por tanto, pensamos que el mejor modelo para representar las carac-
terı́sticas buscadas corresponde al ejemplo 5. Consideramos a continuación
el sistema de dicho ejemplo para un valor de la tensión a > 0 genérico.
(
ḃ = x,
(2.2)
ẋ = −(x3 − ax + b),
Cuando a = 0, los autovalores asociados al sistema tienen parte real nula,
por tanto, el teorema de Hartman-Grobman no nos dice qué tipo de equilibrio
hay en el (0, 0), aunque como se ve numéricamente en la figura 2.15 parece
que hay un foco estable y el corazón pararı́a en diástole. Este ejemplo se
corresponde con el experimento de Rybak presentado en el capı́tulo anterior.

Figura 2.15: Ejemplo del sistema (2.2) con a = 0

26
 Para valores de a > 0 pequeños, las órbitas tienden a un ciclo lı́mite con
amplitud en x pequeña, por lo que prácticamente el movimiento se queda
en diástole como se aprecia en la figura 2.16. Más adelante, aplicaremos el
teorema de Liénard para justificar la existencia de dicha órbita periódica.
A pesar de que Zeeman propone un modelo con un equilibrio estable, como
dice la cualidad (i), consideramos que es mejor que el equilibrio sea inestable.
De esta manera, aparece un ciclo lı́mite estable, lo que permite representar
mejor la periodicidad del ciclo cardı́aco.

Figura 2.16: Ejemplo del sistema (2.2) con a = 0.2

Para a grande, se ve en la figura 2.17 que la oscilación en x es más amplia

y aunque el equilibrio sea inestable como querı́amos, eso hace que la exten-
sión muscular del corazón sea mayor. De esta manera, la ley de Starling se
cumplirı́a con este modelo. De hecho, para valores de a demasiado grandes,
la extensión puede no ser soportada por las fibras del corazón, provocando
un fallo cardı́aco.

27
 Figura 2.17: Ejemplo del sistema (2.2) con a = 2

Nuestra propuesta como modificación del modelo de Zeeman

Por otro lado, si queremos que el equilibrio se encuentre más cerca de la
zona de diástole que de la de sı́stole, podemos considerar el siguiente sistema
de ecuaciones, que es una generalización del modelo 6 de Zeeman, para una
tensión a general,
(
ḃ = x − x0 ,
(2.3)
ẋ = −(x3 − ax + b),
y con la principal diferencia de que
r
a
0 < x0 < . (2.4)
3
Aunque no podemos aplicar el teorema de Liénard a este sistema y, por
tanto, no podemos justificar la existencia de un ciclo lı́mite estable, numéri-
camente podemos ver en la figura 2.18 que también hay una órbita cerrada.

28
En este caso, el punto de equilibrio es (b0 , x0 ) = (−x30 + ax0 , x0 ) y la matriz
jacobiana resultante del sistema anterior evaluada en dicho punto es
 
0 1
J(b0 , x0 ) =  1 −3x20 + a  .
−
 
Los autovalores correspondientes son
p
−(3x20 − a) ± (3x20 − a)2 − 4
λ= .
2
p
Puesto que estamos considerando el caso 0 < x0 < a/3, se tiene que
3x20 − a < 0. Por tanto, ambos autovalores son reales positivos y el punto de
equilibrio resulta ser un nodo inestable.

Figura 2.18: Ejemplo del sistema (2.3) con x0 = 0.5 y a = 1

Si aumenta la tensión, aumenta también la amplitud de extensión de las

fibras musculares. Por tanto, el latido serı́a más fuerte, explicando de nuevo

29
la ley de Starling. Esto se ve gráficamente en la figura 2.19.

Figura 2.19: Ejemplo del sistema (2.3) con x0 = 0.5 y a = 2

La diferencia más importante a nivel cualitativo respecto del caso anterior

(2.2) es que cuando la tensión a disminuye, x0 deja de cumplir la condición
(2.4) y el equilibrio (ax0 − x30 , x0 ) pasa a ser un sumidero donde el corazón
acabarı́a parando, como ocurre en la figura 2.20. En ese sentido, el valor
a = 3x20 se puede interpretar como la tensión umbral, a partir de la cual el
corazón late de forma periódica.

30
Figura 2.20: Ejemplo del sistema (2.3) con x0 = 0.5 y a = 0.5

31
Capı́tulo 3

Teoremas de existencia o no de
soluciones periódicas en el
plano

A continuación se exponen una serie de teoremas que justifican o niegan

la existencia de órbitas cerradas en un sistema de ecuaciones diferenciales.
Más adelante estos teoremas serán aplicados explı́citamente a los modelos
estudiados en el capı́tulo anterior.

Se considera el siguiente sistema

dx


 = F (x, y),
dt (3.1)
 dy = G(x, y)

dt
1 2
con F, G ∈ C (R ).

Los siguientes teoremas son debidos a Poincaré-Bendixson. Uno de ellos

proporciona un criterio positivo para la existencia de órbitas periódicas [3] y
otro un criterio negativo [5]. Daremos solo la demostración del que propor-
ciona el criterio negativo.
Teorema 3.1. Si ∂F/∂x + ∂G/∂y es siempre positivo o siempre negativo
en una cierta región del plano de fases, entonces el sistema (3.1) no tiene
ninguna trayectoria cerrada en esa región.

Demostración. Vamos a razonar por reducción al absurdo. Es decir, su-

ponemos que una determinada región tiene una trayectoria cerrada Γ =

32
[x(t), y(t)], cuyo interior es R. Aplicando el teorema de Green y la hipótesis
se tiene que
Z Z Z
∂G ∂F
(F dy − Gdx) = ( + )dxdy 6= 0.
Γ R ∂y ∂x
Pero a lo largo de Γ se sabe que

dx = F dt,
dy = Gdt.
Por tanto
Z Z T
(F dy − Gdx) = (F G − GF )dt = 0.
Γ 0
Se llega a una contradicción, ası́ que en la región anterior no hay trayec-
torias cerradas.

Teorema 3.2 (Teorema generalizado de Poincaré-Bendixson). Sean F, G ∈

C 1 (E), donde E es un subconjunto abierto de R2 . Sea Γ = [x(t), y(t)] una
trayectoria del sistema (3.1) contenida en un subconjunto compacto F de E
para t ≥ t0 . Entonces, si (3.1) tiene un número finito de puntos crı́ticos en
F , resulta que el conjunto de puntos de acumulación de Γ es un punto crı́tico
de (3.1), es una órbita periódica de (3.1) o bien consiste en un número finito
de puntos crı́ticos {p1 , ..., pn } de (3.1) y en un conjunto numerable de órbitas
lı́mite de (3.1) cuyos puntos de acumulación cuando t → ∞ y t → −∞
pertenecen a {p1 , ..., pn }.
A continuación, para ejemplificar dicho teorema, consideramos el siguiente
sistema no lineal [5],
dx


 = −y + x(1 − x2 − y 2 ),
dt (3.2)
 dy = x + y(1 − x2 − y 2 ),

dt
ejemplo de un sistema que tiene una trayectoria cerrada aislada, es decir, no
hay más trayectorias cerradas próximas a ella.

Hacemos un cambio de coordenadas a polares r y θ, con x = rcosθ e

y = rsenθ. Se tiene que x2 + y 2 = r2 y θ = tg −1 (y/x), ası́ que obtenemos las
siguientes fórmulas:
dx dy dr
x +y =r
dt dt dt
33
y
dy dx dθ
x −y = r2 .
dt dt dt
Multiplicando la primera ecuación del sistema (3.2) por x, la segunda por y
y sumando, se obtiene
dr
r = r2 (1 − r2 ).
dt
Del mismo modo, si ahora se multiplica la segunda por x, la primera por y
y se restan, se llega a
dθ
r2 = r2 .
dt
El único punto crı́tico se encuentra para r = 0. Como lo que nos interesa
son las trayectorias, podemos suponer que r > 0. De esta manera, se tiene
que las igualdades anteriores se convierten en
dr


 = r(1 − r2 ),
dt (3.3)
dθ
= 1.


dt
El sistema anterior es fácil de resolver por separado y su solución es

r = √ 1

,
1 + ce−2t (3.4)
 θ = t + t0 .


Por tanto, la solución general del sistema (3.2) es


cos(t + t0 )
x = √ ,


1 + ce−2t
(3.5)
sin(t + t0 )
y = √

 .
1 + ce−2t
Analizando geométricamente (3.5) se tiene que para c = 0, la solución se
corresponde en (3.3) con r = 1 y θ = t + t0 , las cuales describen la trayectoria
circular cerrada x2 + y 2 = 1. Si c < 0 entonces r > 1 y además se tiene que
r → 1 cuando t → ∞. Por otro lado, si c > 0, resulta r < 1 y de nuevo r → 1
cuando t → ∞. Esto demuestra que existe una sola trayectoria cerrada (que
es la que corresponde a r = 1) a la que tienden el resto de trayectorias en
forma de espiral cuando t → ∞, ya sea por dentro o por fuera de ella.

34
 En general, hacer este tipo de razonamientos no resulta sencillo, y en este
sentido los teoremas 3.1 y 3.2 pueden ayudar a probar la existencia o no de
trayectorias cerradas.
Como hemos dicho antes, el sistema (3.2) tiene un único punto crı́tico en
el punto (0, 0). Además, es fácil ver que la región R limitada por los cı́rculos
r = 1/2 y r = 2 no contiene puntos crı́ticos. Antes hemos visto que
dr
= r(1 − r2 )
dt
para r > 0. Esto quiere decir que sobre el cı́rculo interior r = 1/2 se tiene
dr/dt > 0 y en cambio, en el cı́rculo exterior r = 2 se tiene que dr/dt < 0. Por
tanto, el vector que define el campo del sistema (3.2) apunta hacia la región
R en todos los puntos de la frontera de R. Es decir, cualquier trayectoria
que pase por un punto frontera entrará en R y permanecerá dentro cuando
t → ∞. Por el teorema 3.2 de Poincaré-Bendixson, se deduce que, como en
R no hay puntos crı́ticos, los puntos de acumulación de cualquier órbita en
R deben corresponder con los de una órbita periódica. Por tanto, existe una
trayectoria cerrada C0 en R. Concretamente C0 se corresponde con el cı́rculo
r = 1 como ya hemos visto anteriormente.

El teorema 3.2 es, en general, difı́cil de aplicar, a pesar de ser satisfactorio

a nivel teórico. Existe un criterio más práctico que asegura la existencia de
órbitas cerradas para cualquier ecuación de la forma

ẍ + f (x)ẋ + g(x) = 0, (3.6)

conocida como ecuación de Liénard.
Rx
Si denotamos F (x) = 0
f (s)ds, esta ecuación es equivalente al siguiente
sistema

ẋ = y − F (x),
(3.7)
ẏ = −g(x).
Por tanto, hablar de una trayectoria cerrada de (3.7) equivale a hablar
de una solución periódica de (3.6).

Consideramos ahora el siguiente teorema, conocido como teorema de

Liénard, que garantiza la existencia de una solución periódica para este sis-
tema [3].
Teorema 3.3 (Teorema de Liénard). Supongamos que las funciones F y g
anteriores satisfacen las siguientes propiedades:

35
 (i) F, g ∈ C 1 (R).
(ii) g es una función impar en x que cumple que xg(x) > 0 para x 6= 0 y
F es también una función impar tal que F 0 (0) < 0.
(iii) F tiene un único cero positivo en x = a y crece monótonamente hacia
el infinito para x ≥ a cuando x → ∞.
Entonces el sistema (3.7) correspondiente a la ecuación de Liénard tiene
exactamente una órbita periódica, que es además un ciclo lı́mite estable.
Demostración. El único punto crı́tico del sistema (3.6) es el origen. Además,
el sistema es invariante al cambio de variables (x, y) ↔ (−x, −y), ya que si
(x(t), y(t)) describe una trayectoria, también lo hará (−x(t), −y(t)). Como
para x > 0, la variable y decrece, no puede haber una órbita cerrada en
ese semiplano y, por el mismo motivo, tampoco en x < 0, donde y solo
crece. Por tanto, cualquier posible órbita cerrada debe cortar en dos puntos
al eje y. Además, como en la semirrecta positiva del eje y el flujo va hacia la
derecha, mientras que en la negativa va hacia la izquierda, la órbita cerrada
necesariamente contiene el origen.
Como su simétrica respecto del origen también debe ser trayectoria y
ambas no pueden cortarse, si Γ es una trayectoria cerrada, necesariamente es
simétrica respecto del origen.
Por otro lado, a lo largo de la curva y = F (x), el flujo es vertical hacia
abajo en el semiplano x > 0 y vertical apuntando hacia arriba en el semiplano
x < 0. Debido a lo anterior, cualquier trayectoria que empiece en un punto
P0 = (0, y0 ) en el eje positivo y, cortará a la curva y = F (x) verticalmente en
algún punto P2 = (x2 , y2 ) y cortará también el eje negativo y horizontalmente
en un punto P4 = (x4 , y4 ). Debido a la simetrı́a descrita anteriormente, Γ es
una trayectoria cerrada si y solo si −y0 = y4 .
Consideramos la función de energı́a
y2
u(x, y) = + G(x),
2
Rx
donde G(x) = 0 g(s)ds. Entonces la propiedad anterior equivale a que
u(0, y0 ) = u(0, y4 ).
Sea A el arco P0 P4 de Γ, definimos
Z
φ(α) := du = u(0, y4 ) − u(0, y0 ),
A
donde α = x2 , es decir, la abscisa del punto P2 . Por tanto, por la necesaria
simetrı́a respecto al origen de las posibles órbitas cerradas, Γ es cerrada si y
solo si φ(α) = 0.

36
 Nuestro objetivo ahora es probar que φ(α) tiene exactamente un cero en
α = α0 y que α0 > a, siendo a el punto de corte de la curva y = F (x) con el
eje x.
Sabemos que a lo largo de Γ,
du = g(x)dx + ydy = − [y − F (x)] dy + ydy = F (x)dy. (3.8)
Si 0 < x < a, se tiene F (x) < 0 y dy = −g(x)dt < 0, ası́ que du > 0.
Esto implica que φ(α) > 0, es decir, u(0, y4 ) > u(0, y0 ) si α ≤ a. Entonces,
cualquier trayectoria Γ que corte la curva y = F (x) en un punto P2 con
0 < x2 = α ≤ a no es cerrada.
Veamos ahora que para α ≥ a, la función φ(α) es monótona decreciente.
Y decrece desde el valor positivo φ(a) hasta −∞ según α crece en el intervalo
[a, ∞), con lo que quedarı́a probado que φ(α) tiene un único cero α0 > a y
por tanto el teorema.
Para α > a, dividimos el arco A en tres partes A1 = P0 P1 , A2 = P1 P3 y
A3 = P3 P4 , y definimos las siguientes funciones
Z Z Z
φ1 (α) := du, φ2 (α) := du, φ3 (α) := du,
A1 A2 A3
de forma que
φ(α) = φ1 (α) + φ2 (α) + φ3 (α).

37
 A lo largo de Γ se tiene
   
dy yg(x) −g(x)F (x)
du = g(x) + y dx = g(x) − dx = dx. (3.9)
dx y − F (x) y − F (x)

A lo largo de A1 y de A3 tenemos que F (x) < 0 y g(x) > 0, lo cual implica

que −g(x)F (x) > 0 y por otro lado
dx
= dt > 0.
y − F (x)

Como consecuencia φ1 (α), φ3 (α) son positivos.

En cambio, a lo largo de A2 se tiene que tanto F (x) como g(x) son posi-
tivos y dx/(y − F (x)) = dt > 0, ası́ que φ2 (α) < 0.

Como las trayectorias no se cortan, si aumentamos α, el arco A1 sube y

el arco A3 baja.
A lo largo de A1 , los lı́mites de integración en x quedan fijos en x0 = 0 y
en x1 = a; ası́ que para cada x fijo en [0, a], aumentar α eleva el arco A1 lo
que provoca un incremento en la y y eso hace que el integrando (3.9) descrito
anteriormente decrezca, es decir, φ1 (α) decrece.
De forma similar, a lo largo de A3 , los lı́mites de integración permanecen
fijos en x3 = a y en x4 = 0. Para cada x ∈ [0, a], aumentando α baja el
arco A3 y decrece y. El integrando (3.9) decrece y por tanto decrece también
φ3 (α), ya que
Z 0 Z a 
−g(x)F (x)  −g(x)F (x) 
φ3 (α) = dx =  y − F (x)  dx.
 
a y − F (x) 0

A lo largo de A2 podemos considerar du = F (x)dy. Como las trayectorias

no se cruzan, aumentar α hace que el arco A2 se mueva hacia la derecha. Los
lı́mites de integración en y son y1 (α) e y3 (α), con y1 (α) aumentando cuando
α crece y y3 (α) disminuyendo. Para cada y fijo, aumentar α significa que el
valor de x que le corresponde es mayor y por tanto mayor es el valor de F (x)
y, como
Z y1 (α)
φ2 (α) = − F (x)dy,
y3 (α)

entonces φ2 (α) decrece.

Por tanto, hemos probado que para α ≥ a, la función φ es monótona
decreciente en α. Falta ver que φ2 (α) −→ −∞ cuando α → ∞.

38
 Como hemos dicho antes, a lo largo de A2 , du = F (x)dy = −F (x)g(x)dt <
0. Ası́ que considerando el módulo de φ2 (α) se tiene
 Z  Z
 y1 (α)  y1 (α) Z y1 (α)
|φ2 (α)| = − F (x)dy  = F (x)dy > F (x)dy,
 
 y3 (α)  y3 (α) y3 (ᾱ)

para un cierto ᾱ fijo con a < ᾱ < α, ya que si ᾱ < α, y3 (ᾱ) > y3 (α) y el
intervalo de integración se hace más pequeño. Además, como las trayectorias
no se cortan, si α → ∞ el punto P1 estará más arriba y y1 (α) → ∞, ya que
y1 (α) > y2 (α) = F (α) y F crece monótonamente hacia ∞. Por tanto, queda
probado que φ2 (α) −→ −∞ si α → ∞.
Puesto que φ(α) es una función continua que decrece monótonamente
desde el valor positivo φ(a) hasta −∞ cuando α crece en el intervalo [a, ∞],
se sigue que φ(α) = 0 en un solo valor α = α0 en (a, ∞).
Es decir, el sistema (3.7) tiene una única trayectoria cerrada Γ0 que pasa
por (α0 , F (α0 )).
Dado que φ(α) < 0 para α > α0 y por la simetrı́a descrita al principio de
la demostración, veremos que si α 6= α0 , los puntos de intersección de Γ con
el eje y se acercan a Γ0 , es decir, Γ0 es un ciclo lı́mite estable.

Para probar esto, notemos que, para α > α0 , |y4 | debe ser menor que y0 ,
pero por la simetrı́a de la ecuación, la órbita continuarı́a por el semiplano

39
izquierdo de forma que el siguiente corte con el eje y deberı́a quedar por
debajo de |y4 | para que no corte a la parte simétrica respecto del origen
(lı́nea discontinua en la figura) de la parte de la órbita ya dibujada en el
semiplano x > 0. Procediendo recursivamente la órbita da vueltas en torno a
Γ0 . Por otro lado, en el interior, si g 0 (0) > 0 también las órbitas tienden a Γ0
por el teorema de Poincaré-Bendixson, ya que el equilibiro (0, 0) es inestable
al ser los autovalores asociados a la matriz del sistema
p
−F 0 (0) ± F 0 (0)2 − 4g 0 (0)
λ1,2 = ,
2
ambos valores positivos por la hipótesis F 0 (0) < 0. En el posible caso en que
g 0 (0) = 0, se podrı́a razonar como para la parte exterior del ciclo lı́mite.

3.1. Aplicación de los teoremas a los ejemplos

del capı́tulo 2
En esta sección, vamos a aplicar los teoremas anteriores al sistema de
ecuaciones siguiente, que se corresponde al ejemplo 5,
(
ḃ = x,
(3.10)
ẋ = −(x3 − ax + b)
para un valor de tensión a > 0.

En primer lugar aplicamos el teorema 3.1, donde las funciones F y G del

teorema serán ahora F (x, b) = −(x3 − ax + b)/ y G(x, b) = x. Por tanto,

∂F ∂G −3x2 + a
+ = .
∂x ∂b 
Es decir, el teorema justifica que no puedephaber trayectoriaspcerradas
en las regiones correspondientes a x ∈ (−∞, − a/3) ni en x ∈ ( a/3, ∞)
puesto que en ellas
∂F ∂G
+ < 0.
∂x ∂b
De la misma manera,
p tampoco
p hay órbitas periódicas en la región correspon-
diente a x ∈ (− a/3, a/3), ya que en esa región

∂F ∂G
+ > 0.
∂x ∂b

40
 A continuación aplicaremos el teorema de Liénard al sistema (3.10).

Derivando la segunda ecuación y sustituyendo la primera se obtiene la

ecuación de segundo orden
3x2 − a 1
ẍ + ẋ + x = 0, (3.11)
 }
| {z 
|{z}
f (x) g(x)

con forma de ecuación de Liénard, donde f (x) = (3x2 − a)/ y g(x) = x/
serán las funciones f y g del teorema. Por tanto, F (x) = (x3 − ax)/.

Veamos que este ejemplo cumple las condiciones de la hipótesis:

(i) Es claro que F, g ∈ C 1 (R).
(ii) Para comprobar que g es una función impar, calculamos
1
−g(x) = − x = g(−x).

Además, xg(x) = x2 / > 0 para x 6= 0.
En el caso de F ,
x3 − ax
−F (x) = − = F (−x)

y queda probado que es también una función impar. Por último,
3x2 − a
F 0 (x) =

y por tanto F 0 (0) = −a/ < 0, ya que a,  > 0.
(iii) Como F se puede descomponer de esta forma
√ √
x3 − ax x(x − a)(x + a)
F (x) = = ,
 
√
F tiene un único cero positivo en x =p a > 0.pDado que F 0 (x) =
(3x2 − a)/ es p para x ∈ (−∞, − a/3) ∪ ( a/3, ∞) y negativa
ppositiva √
para x ∈ (− a/3, a/3), se tiene que F (x) es creciente para x ≥ a
como pedı́a la condición. Además crece hacia el infinito cuando x → ∞.
Por tanto, se cumplen todas las condiciones del teorema y podemos con-
cluir que el sistema (3.10) tiene un único ciclo lı́mite estable.

41
3.2. Aplicación de los teoremas a la ecuación
de Van der Pol
Por último, comentaremos que para  = 1 y a = 3, la ecuación (3.11) es
un caso particular de la bien conocida ecuación de Van der Pol. La ecuación
de Van der Pol es una ecuación diferencial de segundo orden homogénea de
la siguiente forma

ẍ + µ(x2 − 1)ẋ + x = 0, (3.12)

con µ una constante positiva que representa el coeficiente de amortiguamien-
to.
El sistema equivalente a la ecuación anterior es

ẋ = y,
(3.13)
ẏ = −x − µ(x2 − 1)y.

La ecuación de Van der Pol surgió para modelar un circuito eléctrico de

serie RLC. Pero más adelante, ha sido usada en otros campos cientı́ficos,
dentro de la fı́sica, la ingenierı́a o la biologı́a. Por ejemplo, para describir el
potencial de acción de las neuronas, en sismologı́a para modelar el compor-
tamiento de dos placas tectónicas en una falla y también se usó este tipo de
oscilador eléctrico no lineal como precursor de las primeras radios comercia-
les.

El oscilador de Van der Pol es un caso especial de oscilador con amorti-

guamiento no lineal. La fuerza de amortiguamiento a la que está sometida
el sistema es µ(x2 − 1)dx/dt. Esto quiere decir que el sistema favorece las
oscilaciones pequeñas y amortigua las oscilaciones grandes.

La principal aplicación del teorema 3.3 de Liénard se produce con la

ecuación de Van der Pol. El único punto crı́tico del sistema es el punto (0, 0).
Por otro lado, se tiene que la matriz jacobiana correspondiente al sistema
(3.13) evaluada en el punto crı́tico es
 
0 1
J(0, 0) = .
−1 −µ

Por tanto, los autovalores asociados son

p
µ ± µ2 − 4
λ= .
2

42
Como µ > 0 se tiene que los dos autovalores tienen parte real positiva, es
decir, el sistema tiene un equilibrio inestable.

Se puede aplicar el teorema 3.3 de Liénard a la ecuación anterior, consi-

derando f (x) = µ(x2 − 1) y g(x) = x. De esta manera, es fácil ver que se
cumplen las condiciones (i) y (ii). Como
1 1
F (x) = µ( x3 − x) = µx(x2 − 3),
3 3
√ √
tiene un solo cero positivo
√ x = 3. Por tanto, F (x) < 0 si x ∈ (0, 3), y
F (x) > 0 para x > 3. Además, F (x) → ∞ cuando x → ∞. Por otro lado,
F 0 (x) = µ(x2 − 1), ası́ que F 0 (0) = −µ < 0. Se cumplen todas las condiciones
del teorema y se concluye que la ecuación de Van der Pol (3.12) tiene una
única solución periódica que además es un ciclo lı́mite.

43
Capı́tulo 4

Simulación de un
electrocardiograma

Un electrocardiograma (ECG) es una prueba que representa la actividad

eléctrica que se produce en cada latido del corazón en función del tiempo. En
dicha prueba se obtiene una gráfica donde aparece una onda que representa
los estı́mulos cardı́acos. En un electrocardiograma correspondiente a un ciclo
cardı́aco normal deben aparecer representadas: la despolarización auricular,
la repolarización auricular, la despolarización ventricular y la repolarización
ventricular que se producen en cada latido del corazón.

Figura 4.1: Componentes de un electrocardiograma normal

El impulso cardı́aco se genera en el nodo sinoauricular (SA), situado entre

44
la aurı́cula derecha y la vena cava. Dicho nodo genera impulsos periódicamen-
te, es decir, actúa como un marcapasos fisiológico del corazón. Estos impulsos
se propagan por todas las células del corazón y la despolarización se propaga
hacia el nodo auriculoventricular (AV). El nodo AV actúa como sustituto del
nodo SA si éste deja de funcionar. Desde el nodo AV, el impulso se propaga
por el haz de His y las fibras de Purkinje que generan la contracción de los
ventrı́culos [6].

Figura 4.2: Esquema de los tres osciladores en el caso normal

La señal que genera un ECG normal está compuesta por una onda P, un
intervalo PR, un complejo QRS, un intervalo ST y otra onda T. En la figura
4.1 se ven estos elementos. Cada pico y cada valle de esta figura representan
cierta actividad eléctrica del corazón. La onda P corresponde a la despolari-
zación auricular y coincide con el primer momento de la diástole. El principio
de la sı́stole se correponde con el complejo QRS de despolarización ventricu-
lar. El intervalo ST se corresponde con parte del proceso de repolarización
ventricular. La onda T se genera por la repolarización ventricular.

Volviendo al sistema (2.3) que hemos considerado más adecuado para

modelar el movimiento del corazón,
(
ḃ = x − x0 ,
ẋ = −(x3 − ax + b),

con r
a
0 < x0 < ,
3

45
es posible hallar su solución y dibujar la variable b como función de t.

De esta forma, para x0 = 0.5,  = 1/100 y a = 1, se obtiene la gráfica

de la figura 4.3, correspondiente a una aproximación del electrocardiograma
que generarı́a dicho modelo.

Figura 4.3: Electrocardiograma generado a partir del sistema (2.3)

En la gráfica generada se puede ver la periodicidad de la variable b aso-

ciada al control electroquı́mico en el latido del corazón, pero no se pueden
apreciar los picos y valles caracterı́sticos de un electrocardiograma normal.

A continuación estudiaremos un modelo mucho más elaborado que pro-

pone Gois en [1], y que se describe también en [4], que simula mejor el
comportamiento del ECG. Gois propone un modelo matemático que descri-
be el ritmo cardı́aco considerando tres osciladores modificados de Van der
Pol acoplados con un retardo de tiempo.

46
 Gois hace uso de la ecuación de Van der Pol para describir el modelo,
ya que el comportamiento de dicha ecuación se ajusta a las caracterı́sticas
que debe cumplir el entorno biológico del corazón, como por ejemplo y como
hemos visto en otros capı́tulos, la existencia de un ciclo lı́mite.

El uso de tres osciladores se debe a que cada uno de ellos representa

una región diferente del corazón. Un oscilador se corresponde con el nodo
sinoauricular, otro con el nodo auriculoventriuclar y el tercero se correspon-
de con el complejo His-Purkinje. Más concretamente, nos centraremos aquı́,
por simplicidad, en el modelo sin retardo, que ya da una buena aproxima-
ción del electrocardiograma. Gois propone el siguiente sistema de ecuaciones
diferenciales para describir la dinámica del corazón,

x˙1 = x2 ,
x˙2 = −aSA x2 (x1 − wSA1 )(x1 − wSA2 ) − x1 (x1 + dSA )(x1 + eSA )
+ ρSA sin(ωSA t) − kSA−AV (x1 − x3 ) − kSA−HP (x1 − x5 ),
x˙3 = x4 ,
x˙4 = −aAV x4 (x3 − wAV1 )(x3 − wAV2 ) − x3 (x3 + dAV )(x3 + eAV )
+ ρAV sin(ωAV t) − kAV −SA (x3 − x1 ) − kAV −HP (x3 − x5 ),
x˙5 = x6 ,
x˙6 = −aHP x6 (x5 − wHP1 )(x5 − wHP2 ) − x5 (x5 + dHP )(x5 + eHP )
+ ρHP sin(ωHP t) − kHP −SA (x5 − x1 ) − kHP −AV (x5 − x3 ),

donde las variables que describen las oscilaciones son (x1 , x2 ), (x3 , x4 ) y
(x5 , x6 ) y, los parámetros son aSA , aAV , aHP , wSA1 , wAV1 , wHP1 , wSA2 , wAV2 ,
wHP2 =, dSA , dAV , dHP , eSA , eAV , eHP , ωSA , ωAV , ωHP , ρSA , ρAV , ρHP , α0 ,
α1 , α3 , α5 , siendo kSA−AV , kSA−HP , kAV −SA , kAV −HP , kHP −SA y kHP −AV los
parámetros que indican el acoplamiento entre los distintos osciladores.

La combinación de las ondas procedentes de cada región del corazón da

lugar a la forma del ECG. Bajo estas consideraciones, parece lógico que los
osciladores acoplados correspondientes a la señal que genera cada región sean
capaces de representar la dinámica general del corazón.

47
 Con los siguientes parámetros que aparecen en [1], y que se muestran en
la siguiente tabla,

aSA = 3 aAV = 3 aHP = 5

wSA1 = 0.2 wAV1 = 0.1 wHP1 = 1
wSA2 = −1.9 wAV2 = −0.1 wHP2 = −1
dSA = 3 dAV = 3 dHP = 3
eSA = 4.9 eAV = 3 eHP = 7
ωSA = 0 ωAV = 0 ωHP = 0
ρSA = 1 ρAV = 1 ρHP = 20
kSA−AV = 0 kAV −SA = 5 kHP −SA = 0
kSA−HP = 0 kAV −HP = 0 kHP −AV = 20
α0 = 1 α1 = 0.1 α3 = 0.05
α5 = 0.4

hemos generado la gráfica correspondiente a un electrocardiograma nor-

mal representando
ECG = α0 + α1 x1 + α3 x3 + α5 x5
en función del tiempo, como se propone en [1].

Figura 4.4: Electrocardiograma normal

48
 La gráfica resultante es la de la figura 4.4. En dicha figura se aprecian
cualitativamente las diferentes partes que debe presentar un electrocardiogra-
ma normal, la señal contiene las ondas más importantes P, QRS y T. Esto
demuestra cómo el modelo es capaz de describir el comportamiento general
de un ECG real.

Por otro lado, el modelo anterior es capaz de simular ciertas patologı́as

cardı́acas que se pueden identificar estudiando el ECG. Si se considera el caso
en el que no existe comunicación entre el nodo SA y el nodo AV, la constante
kAV −SA deberı́a ser nula, eliminando ası́ el acoplamiento entre el primer y el
segundo oscilador. Considerando ahora kAV −SA = 0 y manteniendo el valor
del resto de parámetros como la tabla, el sistema resultante está dirigido
por el nodo AV. Esta patologı́a se denomina aleteo ventricular y presenta la
forma de la figura 4.5.

Figura 4.5: Electrocardiograma aleteo ventricular

Si comparamos las figuras 4.3 y 4.5, se puede ver que el comportamiento

de ambos electrocardiogramas es muy similar a nivel cualitativo. Esto parece

49
lógico, puesto que en el modelo de Zeeman modificado que hemos considerado
para representar la figura 4.3 no hay acoplamiento entre en nodo SA y el nodo
AV, como ocurre en el caso del aleteo ventricular.

50
 Apéndice A

Programas de MATLAB

En este apéndice se muestran las funciones de MATLAB usadas para ge-

nerar las gráficas presentadas a lo largo del trabajo y que han servido como
apoyo para estudiar los sucesivos modelos.

1 f u n c t i o n F = f u n c i o n s i s t e m a ( ˜ ,X)
2 global funcionf funciong
3 x = X( 1 ) ;
4 b = X( 2 ) ;
5 F = [ eval ( funcionf ) ; eval ( funciong ) ] ;

1 % Dibuja e l campo v e c t o r i a l a s o c i a d o a un s i s t e m a
autonomo
2 % p lano db/ dt= g ( x , b ) dx/ dt = f ( x , b ) y s u s s o l u c i o n e s .
3
4 % CAMBIAR l neas 1 7 , 18 y 40 dependiendo d e l s i s t e m a .
5
6 clear
7 global funcionf funciong
8 clf reset
9
10 puntos = 1 6 ;
11 xmin = −2; xmax = 2 ; bmin = −2; bmax = 2 ;
12 i n c r x = ( xmax−xmin ) / puntos ; i n c r b = (bmax−bmin ) / puntos ;
13
14 [ x , b ] = meshgrid ( xmin : i n c r x : xmax , bmin : i n c r b : bmax) ;
15 [ n ,m] = s i z e ( x ) ;

51
16
17 f u n c i o n f = ’ −10∗(x .ˆ3 −2.∗ x+b ) ’ ;
18 f u n c i o n g = ’ x −0.5 ’ ;
19 Heading = [ ’ \n\n\n ’ , ’ P a r metros para e l s i s t e m a
a u t nomo \n\n ’ , ’ dx/ dt = ’ , . . .
20 f u n c i o n f , ’ db/ dt = ’ , f u n c i o n g , ’ \n\n ’ ] ;
21 f p r i n t f ( Heading )
22 f p r i n t f ( ’ \n xmin = %g xmax = %g bmin = %g bmax = %g \n ’
,...
23 xmin , xmax , bmin , bmax)
24
25 f p r i n t f ( ’ \n i n c r x = %g i n c r b = %g \n\n ’ , i n c r x , i n c r b )
26 evalf = eval ( funcionf ) ;
27 evalg = eval ( funciong ) ;
28 L = s q r t ( e v a l f .ˆ2+ e v a l g . ˆ 2 ) ;
29 q u i v e r ( b , x , e v a l g . / L , e v a l f . / L , 0 . 5 ) % Dibuja e l campo
vectorial
30 a x i s ( [ xmin , xmax , bmin , bmax ] )
31 g r i d on
32 x l a b e l ( ’ b v a r i a b l e de c o n t r o l e l e c t r o q u mico ’ )
33 y l a b e l ( ’ x l o n g i t u d de l a f i b r a muscular ’ )
34 T i t l e = [ ’ db/ dt = ’ , f u n c i o n g , ’ dx/ dt = ’ , f u n c i o n f ] ;
35 t i t l e ( Title )
36 hold
37
38 h o l d on
39 a=xmin : 1 / 2 0 0 : xmax ;
40 y = −a . ˆ 3 + 2 . ∗ a ; % D e f i n i r l a v a r i e d a d l e n t a
41 p l o t ( y , a , ’ r ’ ) % Dibuja l a v a r i e d a d l e n t a
42
43 p r e g = i n p u t ( ’ D i b u j o una s o l u c i n con e l r a t n ? ( s /n ) ’
, ’s ’);
44 w h i l e preg ( 1 ) == ’ s ’
45
46 f p r i n t f ( ’ \n L o c a l i c e con e l r a t n la c o n d i c i n i n i c i a l
en l a ventana \n ’ )
47 [ i n i c i a l b , i n i c i a l x ] = g i n p u t ( 1 ) ; % Toma como c o n d i c i n
i n i c i a l donde p i n c h e s con e l r a t n
48
49 [ tiempo , s o l ] = ode15s ( ’ f u n c i o n s i s t e m a ’ , [ 0 , 3 0 ] , [
inicialx , inicialb ]) ;

52
50 h o l d on % El comando ode15s r e s u e l v e l a e c u a c i o n
d i f e r e n c i a l para t =0 ,30 y en e s e punto i n i c i a l
51 p l o t ( s o l ( : , 2 ) , s o l ( : , 1 ) , ’ k ’ , ’ LineWidth ’ , 1 . 5 )
52
53 [ tiempo , s o l ] = ode15s ( ’ f u n c i o n s i s t e m a ’ , [ 0 , − 3 0 ] , [
inicialx , inicialb ]) ;
54 h o l d on % f u n c i o n s i s t e m a toma l a cadena de c a r a c t e r e s
de l a s f u n c i o n e s y l a pasa a ’ numeros ’
55 p l o t ( s o l ( : , 2 ) , s o l ( : , 1 ) , ’ k ’ , ’ LineWidth ’ , 1 . 5 )
56 p r e g = i n p u t ( ’ D i b u j o o t r a s o l u c i n ? ’ , ’ s ’ ) ;
57 end
58
59 p r e g = i n p u t ( ’ D i b u j o una s o l u c i n con v a l o r inicial
por t e c l a d o ? ( s /n ) ’ , ’ s ’ ) ;
60 w h i l e preg ( 1 ) == ’ s ’
61
62
63 fprintf ( ’ \n coordenada x i n i c i a l en l a ventana \n ’ )
64 inicialx = input ( ’ x 0 = ’ ) ;
65 fprintf ( ’ \n coordenada b i n i c i a l en l a ventana \n ’ )
66 inicialb = input ( ’ b 0 = ’ ) ;
67
68 [ tiempo , s o l ] = ode15s ( ’ f u n c i o n s i s t e m a ’ , [ 0 , 3 0 ] , [
inicialx , inicialb ]) ;
69 h o l d on
70 p l o t ( s o l ( : , 2 ) , s o l ( : , 1 ) , ’ k ’ , ’ LineWidth ’ , 1 . 5 )
71 [ tiempo , s o l ] = ode15s ( ’ f u n c i o n s i s t e m a ’ , [ 0 , − 3 0 ] , [
inicialx , inicialb ]) ;
72 h o l d on
73 p l o t ( s o l ( : , 2 ) , s o l ( : , 1 ) , ’ k ’ , ’ LineWidth ’ , 1 . 5 )
74 preg = input ( ’ D i b u j o otra s o l u c i n ? ’ , ’ s ’ ) ;
75 end

1 % Dibuja e l e l e c t r o c a r d i o g r a m a a s o c i a d o a l modelo de
Zeeman m o d i f i c a d o .
2
3 clear all
4 % Condiciones i n i c i a l e s
5 x0 = [ 0 , 0 ] ;
6 t f =40; % Tiempo f i n a l

53
 7 f=@( t , x ) [ −100∗( x ( 1 ) .ˆ3 −x ( 1 )+x ( 2 ) ) ; x ( 1 ) − 0 . 5 ] ; %
Ecuaciones del sistema
8 [ t , x]= ode15s ( f , [ 0 , t f ] , x0 ) ;
9 plot ( t , x (: ,2) )
10 a x i s ( [ 0 t f −0.6 0 . 6 ] )
11 g r i d on
12 xlabel ( ’ t ’ )
13 ylabel ( ’b ’ ) ;
14 t i t l e ( ’ Electrocardiograma ’ )

1 % Genera l a g r a f i c a c o r r e s p o n d i e n t e a un
e l e c t r o c a r d i o g r a m a normal o con
2 % a l e t e o v e n t r i c u l a r en f u n c i o n de como d e f i n a m o s l o s
parametros .
3
4 x0 = [ 0 , 0 . 7 , 0 , 0 . 2 , 0 , 0 . 7 ] ;
5 % D e f i n i c i o n de c o n s t a n t e s comunes
6 a0 =1; a1 = 0 . 1 ; a3 = 0 . 0 5 ; a5 = 0 . 4 ;
7 as a =3; wsa1 = 0 . 2 ; wsa2 = −1.9; dsa =3;
8 aav =3; wav1 = 0 . 1 ; wav2= −0.1; dav =3; eav =3;
9 ahp =5; whp1=1; whp2=−1; dhp=3; ehp =7;
10 ksaav =0; kavhp =0; ksahp =0; khpsa =0;
11 r o s a =1; roav =1; rohp =20;
12
13 % Seleccionar tipo a continuacion
14 t i p o= ’ normal ’ ;
15 e s a = 4 . 9 ; kavsa =5; khpav =20; wsa =0; wav=0; whp=0; t f =40;
16
17 % t i p o =’ a l e t e o v e n t r i c u l a r ’ ;
18 % e s a = 4 . 9 ; kavsa =0; khpav =20; wsa =0; wav=0; whp=0; t f
=40;
19
20 f=@( t , x ) [ x ( 2 ) ;
21 −asa ∗x ( 2 ) ∗( x ( 1 )−wsa1 ) ∗( x ( 1 )−wsa2 )−x ( 1 ) ∗( x ( 1 )+
dsa ) ∗( x ( 1 )+e s a )+r o s a ∗ s i n ( wsa∗ t )−ksaav ∗( x ( 1 )−x ( 3 ) )−
ksahp ∗( x ( 1 )−x ( 5 ) ) ;
22 x(4) ;
23 −aav ∗x ( 4 ) ∗( x ( 3 )−wav1 ) ∗( x ( 3 )−wav2 )−x ( 3 ) ∗( x ( 3 )+
dav ) ∗( x ( 3 )+eav )+roav ∗ s i n ( wav∗ t )−kavsa ∗( x ( 3 )−x ( 1 ) )−
kavhp ∗( x ( 3 )−x ( 5 ) ) ;

54
24 x(6) ;
25 −ahp∗x ( 6 ) ∗( x ( 5 )−whp1 ) ∗( x ( 5 )−whp2 )−x ( 5 ) ∗( x ( 5 )+
dhp ) ∗( x ( 5 )+ehp )+rohp ∗ s i n (whp∗ t )−khpsa ∗( x ( 5 )−x ( 1 ) )−
khpav ∗( x ( 5 )−x ( 3 ) ) ] ;
26 [ t , x]= ode45 ( f , [ 0 , t f ] , x0 ) ;
27 e c g =(a0+a1 ∗x ( : , 1 )+a3 ∗x ( : , 3 )+a5 ∗x ( : , 5 ) ) ;
28 p l o t ( t , ecg )
29 a x i s ( [ 6 t f −0.5 2 ] )
30 g r i d on
31 xlabel ( ’ t ’ )
32 y l a b e l ( ’ECG ’ ) ;
33 Titulo = [ ’ Electrocardiograma ’ , tipo ] ;
34 t i t l e ( Titulo )

55
Bibliografı́a

[1] Gois, Sandra R.F.S.M. y Marcelo A. Savi An analysis of heart

rhythm dynamics using a three-coupled oscillator model, Chaos, Solitons
and Fracctals, Vol. 41, 2009, pp. 2553-2565.

[2] Jones, D.S., M.J. Plank y B.D. Sleeman, Differential equations and
mathematical biology, Chapman and Hall/CRC, segunda edición, Boca
Ratón, 2009.

[3] Perko, L., Differential equations and dynamical sistems, Springer, ter-
cera edición, Nueva York, 2001.

[4] Quiroz Juarez, M.A., Jimenez Ramirez, O. y Vazquez Medi-

na, R. Análisis e implementación de un modelo matemático para gene-
rar señales de electrocardiograma sintéticas, Conference MATECOMPU
2015, Varadero, Matanzas, Cuba.

[5] Simmons, G.F., J.S. Robertson y B.D. Sleeman, Ecuaciones dife-

renciales con aplicaciones y notas, McGraw-Hill, segunda edición, 1993.

[6] Tresguerres, Jesús A. F., Anatomı́a y fisiologı́a del cuerpo humano,

Mc Graw-Hill, primera edición, 2009.

[7] Zeeman, E.C., Differential equations for the heartbeat and nerve impul-
se, Towards a Theoretical Biology, Vol. 4, 1972, pp. 8-67.

56