Artículo original

 

Modelo de propagación de ondas solitarias en el corazón

 

Model for propagation of solitary waves in the heart

 

Ivonne Domínguez¹, Rafael A. Barrio²*, Carmen Varea² y José Luis Aragón³

 

¹ Facultad de Ciencias, Universidad Nacional Autónoma de México (UNAM). Ciudad Universitaria, Deleg. Coyoacán, C.P. 04510, México, D.F.

² Depto. de Física Química, Instituto de Física, Universidad Nacional Autónoma de México (UNAM). Apdo. Postal 20-364, C.P. 01000, México, D.F. E-mail: *barrio@fisica.unam.mx

³ Depto. de Nanotecnología, Centro de Física Aplicada y Tecnología Avanzada, Universidad Nacional Autónoma de México (UNAM). Apdo. Postal 1-1010, C.P. 76000, Querétaro, México.

 

Artículo recibido el 12 de marzo de 2013.
Aceptado el 31 de julio de 2013.

 

Resumen

En la actividad eléctrica cardiaca, distintos tipos de onda viajan a través del corazón. Presentamos un modelo de la actividad eléctrica del corazón, proponemos que los frentes de onda homogéneos que se propagan en el corazón son, de hecho, ondas solitarias, o "solitones". Usamos un conjunto general de ecuaciones de reacción-difusión conocido como el modelo de Barrio-Varea-Aragón-Maini (BVAM)^[1], el cual presenta una abundancia de bifurcaciones no lineales, siendo capaces de encontrar la ruta al caos usando un mapeo de las ecuaciones de amplitud a la dinámica de la ecuación compleja de Ginzburg-Landau. Estudiamos numéricamente la dinámica de los frentes de onda en el modelo BVAM para describir los mecanismos que conducen a la fibrilación en el corazón y comparamos los resultados con datos experimentales.

Palabras Clave: Actividad eléctrica del corazón, caos, sistemas de reacción difusión, solitones.

 

Abstract

In cardiac electrical activity, different types of waves meander through the heart. We present a model of the electrical activity of the heart that proposes that the homogeneous wave fronts propagating through the heart are in fact solitons. We use a general set of reaction-diffusion equations known as the Barrio-Varea-Aragón-Maini (BVAM) model^[1]that presents a wealth of non-linear bifurcations, and we are able to follow the route to chaos, using a mapping of the amplitude equations to the dynamics of the complex Ginzburg-Landau equation. We study the dynamics of wave fronts numerically in the BVAM model to describe the mechanisms leading to heart fibrillation and compare the findings with experimental data.

Key Words: Electrical activity of the heart, chaos, reaction-difusion systems, solitons.

 

Introducción

Durante las últimas décadas, la investigación sobre el corazón ha sido motivada, no sólo por el deseo de descubrir los secretos de este órgano vital, sino también por su importancia clínica.

Para el bombeo eficiente del corazón se requiere que las células que constituyen el tejido muscular cardiaco se contraigan coordinadas en el espacio, de tal forma que el órgano pueda cambiar su forma y cumplir su función de bombear la sangre. Este proceso de contracción está originado por un estímulo eléctrico que se propaga a través de todas las células cardiacas de manera coherente y ordenada. El ciclo rítmico de contracción y relajación de ~10¹⁰ células musculares es controlado por un patrón complejo de activación eléctrica.

El ritmo cardiaco normal puede ser resultado de muy distintas interacciones entre las células que dan origen al funcionamiento sin falla en más de 3 x 10⁹ latidos consecutivos a lo largo de una vida. Es evidente que el proceso dinámico tiene que ser muy robusto y sin pérdidas de energía o desgaste. Por lo tanto, podemos pensar que el proceso debe ser la consecuencia de una bifurcación no-lineal. Sin embargo, para discernir de qué mecanismo se trata, podemos examinar la información que podamos extraer de las fallas y alteraciones en su funcionamiento, ya que cada fenómeno no-lineal tiene una forma característica de seguir un árbol de bifurcaciones.

El estrecho vínculo entre las alteraciones de la actividad eléctrica del corazón y los problemas cardiacos es la base de la capacidad diagnóstica del electrocardiograma o ECG, que es la más antigua herramienta no invasiva para el diagnóstico de enfermedades del corazón. Tres categorías funcionales del corazón; automatismo, excitabilidad y conductividad, se relacionan con las propiedades eléctricas de la masa y de las fibras cardiacas^[2]. Son estas propiedades las que, registradas gráficamente mediante electrodos colocados en la superficie corporal, dan origen al electrocardiograma, el cual se genera debido a que las fibras musculares deben accionar sincrónicamente por zonas, para que el corazón cumpla adecuadamente con su función de bombear la sangre. Esto implica el registro de señales eléctricas del orden de algunos milivolts, sobre cualquier punto de su volumen corporal. El electrocardiograma es entonces la historia temporal de la suma de la actividad eléctrica de las células en un volumen, registrada desde la superficie corporal.

Existen muchos modelos que pueden explicar la actividad oscilatoria coherente del sistema, e incluso reproducir los datos experimentales de un ECG. Sin embargo, casi ninguno puede reproducir las anomalías que se presentan en la muerte cardiaca, como por ejemplo la aparición del estado de fibrilación usualmente observado antes de la pérdida de la función cardiaca.

En este artículo propondremos un modelo simple que cumple con estos requisitos y además revela el hecho sorprendente, pero aceptable a posteriori, de que la excitación eléctrica en el corazón se propaga en el espacio como frentes de onda solitarias ("solitones"). En efecto, este tipo de ondas no lineales tienen la particularidad de que su propagación es no disipativa, hecho muy relevante cuando se considera que el corazón no puede dejar de funcionar, ni por unos momentos, durante toda una vida. Además, debido a la refractoriedad de las células excitables, cuando dos frentes chocan no se suman como las ondas lineales, sino que el estado excitado se comparte con los dos frentes que chocan y cuando se alejan se recupera la forma de ambos frentes. Esto sucede tanto en imágenes experimentales^[3], como en simulaciones dinámicas de los solitones en nuestro modelo (ver animaciones en Barrio^[4]). Demostraremos que este tipo de excitación es susceptible de presentar caos cuando se varían los parámetros del sistema no-lineal. Además, realizaremos cálculos numéricos cuyos resultados pueden ser comparados con datos experimentales que validen las conclusiones de nuestro modelo.

En la siguiente sección revisaremos los modelos que actualmente sí pueden simular tanto el funcionamiento normal, como las fallas que mencionamos, lo cual será necesario para presentar lo original de nuestro modelo.

 

Modelos actualmente usados

Al proponer un modelo de la actividad eléctrica del corazón se tienen que tomar en cuenta ciertos hechos biológicos de suma importancia^[5]: 1) Debe existir un solo marcapasos que coordine en el tiempo las excitaciones que se propagan en el espacio. 2) Estas excitaciones deben estar coordinadas en una cavidad, y se deben acoplar a las modificaciones de forma causadas por la contracción de las células. 3) Las contracciones deben responder a señales externas químicas y mecánicas que sean capaces de modificar su frecuencia y amplitud. 4) El mecanismo debe ser robusto ante ruido y perturbaciones externas. 5) El corazón tiene una forma muy precisa de dejar de funcionar, relacionada con una ruta específica al caos^[6].

Podemos clasificar los modelos existentes en tres categorías: 1) Los modelos fenomenológicos que consideran una interacción no lineal cúbica, como el modelo de FitzHugh-Nagumo^[7,8] y otros^[9]. 2) Los modelos de primera generación, que explican el comportamiento y descripción de las corrientes iónicas, todos ellos basados en el modelo del potencial de acción de Hodgking-Huxley (modelo H-H)^[10]. Por ejemplo, los modelos de Beeler- Reuter'^[11] y Low-Rudy^[12]. Y 3) Los modelos de segunda generación, que hacen énfasis en la descripción detallada de la fisiología celular, como la concentración de iones y funcionamiento de bombas, y están íntimamente ligados a datos experimentales. Describimos a continuación con más detalle el modelo fenomenológico y el de primera generación, que han sido la base de casi todos los actualmente usados.

 

El modelo de FitzHugh-Nagumo

El modelo de FitzHugh-Nagumo^[7] es una simplificación del modelo H-H a dos variables, de tal forma que presenta una bifurcación de Hopf subcrítica que provoca un cambio repentino en el valor de la variable que representa el potencial de membrana. En su formulación original el modelo está dado por:

donde a, b, c₁, c₂ y c₃ son parámetros ajustables. El valor de estos parámetros da un potencial de acción normalizado. La corriente aplicada i_app también se escala para que coincida con este modelo normalizado. Imágenes de υ y ω se muestran en la Figura 1. El potencial de acción mostrado en la figura resulta de aplicar una corriente i_app con resistencia de 0.05 Ohms, y una duracion de t = 50 s hasta t = 60 s.

 

Se han hecho algunas variaciones de este modelo para superar las limitaciones del modelo original. Por ejemplo, un detalle de la formulación original que no es consistente con los datos fisiológicos es que la neurona hiperpolariza en la fase de repolarización. Esto puede verse en la Fig. 1, donde υ alcanza valores muy por debajo del potencial de reposo antes de volver al estado de reposo. Una modificación de las ecuaciones que elimina este problema fue sugerida por Rogers y McCulloch^[13].

 

El modelo de Hodgking-Huxley

Este modelo ha sido el más exitoso para tratar todo tipo de células excitables, incluyendo neuronas, miocitos, células musculares, etc. Fue concebido originalmente para estudiar la propagación de señales eléctricas en el axón gigante de un calamar. El modelo original publicado en 1952, y que concedió el Premio Nobel a sus autores, es un sistema de al menos tres ecuaciones diferenciales parciales acopladas en las que los potenciales de acción se producen debido a la no linealidad de las interacciones entre las corrientes de diversos iones.

El principio físico usado es muy elemental y consiste en la Ley de Ohm. En el funcionamiento normal de una célula, varias corrientes iónicas contribuyen a los cambios en el potencial de membrana. Hodgkin y Huxley^[2] tomaron como premisa el modelo circuital de membrana, donde cada corriente iónica está representada como:

donde I_x es la corriente asociada al ión x, g_x es la conductancia del canal iónico correspondiente, V_m es el potencial de membrana y E_x es el potencial de Nernst del ión x.

El modelo H-H fue desarrollado para una célula aislada; en estas condiciones, ninguna corriente puede fluir hacia el medio intracelular, por tanto I_m= 0, y la evolución del potencial de membrana puede escribirse, a partir de la Ecuación 2 como:

A partir de sus experimentos, observaron que cuando la membrana de una célula nerviosa se despolariza más allá de un valor umbral, se produce en la célula una corriente de entrada transitoria producida por iones de sodio, seguida por una corriente de salida, formada por iones de potasio. Caracterizaron la dependencia con el voltaje de la conductancia de los canales iónicos del axón gigante postulando la existencia de ciertas partículas cargadas que, unidas al canal iónico correspondiente, permitirían o no el flujo de iones por el canal. Cada una de estas partículas controlaría las compuertas de un canal, de forma que la acción individual de cada partícula correspondería a la probabilidad de apertura (z) o cierre (1 - z) de cada compuerta (0 ≤ z ≤ 1). Cada compuerta puede pasar del estado abierto al cerrado o viceversa, siguiendo una ley dependiente del voltaje de primer orden:

donde la variable α_z es el intervalo de tiempo en el que el canal iónico está abierto y β_z el intervalo de tiempo en el que está cerrado. Éstas se obtienen empíricamente y tienen una dependencia no lineal de funciones exponenciales de V_m. Las Ecuaciones 2, 3 y 4 son, por lo tanto, un sistema de ecuaciones diferenciales acopladas, cuya solución, bajo condiciones de estado estacionario, se caracteriza por sus valores asintóticos z_∞ y su constante de tiempo τ_z.

Finalmente, los modelos de segunda generación, proveen una liga entre los mecanismos celulares y el comportamiento colectivo, pero son complicados y particulares, sin embargo, daremos una breve descripción de los avances obtenidos paralelamente a los obtenidos experimentalmente. Denis Noble^[14] publicó el primer modelo matemático de célula cardiaca, en el que se describía el potencial de acción de las fibras de Purkinje, a partir de la adaptación de las ecuaciones del modelo H-H. Más de una década después, en 1975, McAllister et al.^[15] mejoraron el modelo de Noble mediante la incorporación de nuevas evidencias experimentales sobre las propiedades de la membrana celular. Seguidamente, Beeler y Reuter_[11] desarrollaron en 1977 el primer modelo de miocito ventricular, reformulado después por Luo y Rudy en 1994^[12,16]. Si bien la mayor parte de los modelos celulares cardiacos estaban dirigidos a representar el comportamiento de las células ventriculares, unos pocos modelos fueron desarrollados también para las células auriculares. Las células del nodo sinoauricular fueron las primeras células auriculares en ser descritas matemáticamente. En1998, Nygren et al.^[17] publicaron un modelo de miocito auricular, reformulado en 1999 por Courtemanche et al.^[18,19].

 

Modelo BVAM

El modelo de reacción-difusión que estudiamos en este trabajo es el llamado BVAM (por Barrio-Varea-Aragón- Maini)'^2]. Este modelo, a diferencia de otros, fue concebido no para modelar un sistema o mecanismo específico, sino que se obtuvo como la forma más general de un sistema de reacción-difusión con dos variables, que respeta la conservación de masa.

Consideremos la expresión más general de un sistema de reacción-difusión^[20]:

donde υ y ν son los morfógenos con coeficientes de difusión D_υ y D_ν, respectivamente. F y G son funciones que describen la cinética química y son, en general, no lineales y no existe un criterio general para establecerlas, sus expresiones dependen del problema químico o biológico bajo estudio.

En el modelo BVAM, los términos que representan a la cinética química, F y G, se obtienen suponiendo que son funciones no lineales analíticas y se desarrollan en serie de Taylor hasta orden tres. Varios coeficientes de la expansión en serie se eliminan para garantizar la conservación de masa y el resultado se puede escribir, en forma adimensional, como:

donde η es un parámetro que impone la escala espacio temporal del sistema. Las tasas de crecimiento de υ y ν están descritas en términos de los parámetros f, g y h, que son convenientes para explorar la región de parámetros de interés en este trabajo. La razón entre las intensidades de los términos cuadráticos y cúbicos está dada por C.

Los criterios para obtener una inestabilidad de Turing^[21], imponen restricciones sobre los coeficientes de difusión y los parámetros f, g, h y C. Se ha mostrado también que, además de los patrones espaciales estacionarios (patrones de Turing), que han sido usados en múltiples ocasiones para modelar procesos de pigmentación en la piel de los animales y otros sistemas biológicos (ver por ejemplo Barrio et al.^[1]), en el espacio de parámetros se pueden encontrar también varias bifurcaciones de Hopf y de Hopf-Turing^[22]. En otra región de parámetros se encuentra una bifurcación que produce patrones espaciales de puntos de radio constante, que no es una bifurcación de Turing^[23].

De gran interés para nuestro trabajo se encuentra una región en el espacio de parámetros, cercana a una bifurcación de Hopf, en donde se obtienen frentes de ondas viajeras, cuya velocidad depende de la curvatura local^[24]. Debido a este hecho, en un dominio bidimensional plano los frentes giran alrededor de los puntos donde la curvatura cambia de signo.

Otro aspecto interesante del modelo es que existe una región de parámetros que produce soluciones con caos. Usando el parámetro C , que indica la razón entre los términos cuadráticos y cúbicos en la cinética, se encuentran patrones espaciales de Turing cuya amplitud oscila en el tiempo^[25]. Variando el parámetro de bifurcación, las oscilaciones doblan su período, subsecuentes variaciones del parámetro de bifurcación, producen una región de oscilaciones cuasiperiódicas (hay por lo menos dos frecuencias de oscilación inconmesuradas entre sí^[26]), también llamada bifurcación de Toro, y eventualmente las oscilaciones se vuelven caóticas. Esta ruta al caos se conoce como el mecanismo de Ruelle-Takens-Newhouse^[27].

En este artículo analizaremos más profundamente esta ruta al caos, empezando con un sistema que presente ondas solitarias, con el objeto de comparar los resultados teóricos con los procesos de fibrilación que producen la muerte cardiaca. Para esto, comenzaremos por analizar las propiedades del sistema y la posibilidad de obtener una ruta al caos. En el Apéndice mostramos el método matemático que seguimos para este análisis. La clave es mostrar que las ecuaciones de amplitud del modelo BVAM se pueden transformar en una ecuación de Landau-Ginzburg con coeficientes complejos (Eq. VIII.25), la cual presenta caos.

 

Aplicación a la actividad eléctrica cardiaca Región en el espacio de parámetros para Ondas Solitarias

En ausencia de difusión, en el espacio fase (υ, ν), el modelo BVAM tiene los tres puntos fijos mostrados en la ecuación VIII.15. Si h = -1 todos ellos colapsan en el punto (υ₀ , ν₀) = (0, 0). En el caso h 4 1, diferentes elecciones de parámetros pueden llevarnos a diferentes configuraciones de los puntos de equilibrio en el espacio fase, por ejemplo un punto de silla en el origen (0, 0) rodeado de dos puntos fijos estables, con una trayectoria heteroclínica entre ellos. Esto significa que el sistema no lineal puede tener una bifurcación a un sistema biestable (dos puntos de equlibrio estable) en que algunas regiones del dominio físico serán atraídas por uno de los puntos fijos estables y otras serán atraídas hacia el otro. En esta situación se producen interfases entre estas regiones en las que la concentración de los morfógenos cambia abruptamente, provocando un perfil sigmoidal al atravesar dichas regiones.

La bifurcación puede estudiarse con un análisis lineal de las cinéticas F y G de la ecuación 6 (por ejemplo vea Varea et al.^[24]) y el resultado se muestra en la Fig. 2.

La región de parámetros de interés es el área entre la curva punteada y la continua de la Fig. 2. La región presenta oscilaciones regulares de alta frecuencia (taquicardia) causadas por la existencia de un ciclo límite que encierra a los puntos fijos, los cuales son inestables. En la frontera azul, los puntos fijos se vuelven estables y la cuenca de atracción crece al cruzar la curva, donde se tiene un sistema biestable con frentes de ondas viajeras. Por ejemplo, un punto de la frontera es (g_c, f_c) = (0.165, 0.5043), para C = 0, h = -2.5 y η= 0.5.

Estos solitones se mueven con velocidad proporcional a la curvatura local del perfil y, como consecuencia tiene un movimiento giratorio alrededor de los puntos donde cambia la curvatura del frente. Además, muy cerca de la línea punteada la distancia entre frentes es constante, pero si uno se aleja de la frontera de la transición entonces uno obtiene franjas delgadas y franjas anchas alternadas, como se muestra en el recuadro inferior de la Fig. 3.

En el apartado siguiente mostraremos los resultados de cálculos numéricos efectuados variando los parámetros alrededor de la región de interés.

 

Cálculos numéricos en un plano

Para explorar la región de parámetros de interés, integramos las ecuaciones (6) usando una variación del método de gradiente conjugado (Conjugate Gradient Method), con condiciones periódicas a la frontera y con intervalo de tiempo Δt = 0.01 en un entramado bidimensional (matriz) de 120 x 120. Este valor del paso del tiempo asegura la estabilidad del proceso^'28] y la confiabilidad de los resultados. En la Fig. 3 mostramos los patrones obtenidos numéricamente para varios puntos en la región de interés.

Los recuadros muestran la concentración de υ en un código de color en donde el rojo indica concentración alta y azul baja. Esto puede asociarse al estado de polarización de las células en una región. Observe que para valores grandes de f aparecen ondas espirales, que se mueven lentamente, pero cambian de color rápidamente en el tiempo, lo cual resulta en un espectro con dos frecuencias distintas. A la izquierda de la línea azul se obtienen ondas espirales del tipo Belousov-Zhabotinsky^[21] con diversas frecuencias. Al aproximarse a la bifurcación de Hopf, se observa que aparecen gotas y "ojos", cuya descripción y análisis está fuera del propósito de este artículo. Si en esta región se cambia la razón entre los coeficientes de difusión, entonces los patrones dejan de oscilar y se obtiene un patrón espacial fijo, del tipo Turing.

La región en que se presentan solitones está a la derecha de la línea azul para valores de g ~ 0.165. Al disminuir el valor de f las oscilaciones dejan de ser regulares. Típicamente se observa que aparece un doblamiento de período, la forma de onda cambia y además pueden aparecer algunas componentes de Fourier inconmensuradas con las ya existentes en el espectro (cuasiperiodicidad).

El funcionamiento de un corazón latiendo regularmente puede ser simulado por latidos regulares y robustos. En este trabajo proponemos que la actividad eléctrica en un corazón sano puede simularse con solitones regulares. En la Fig. 4(A) mostramos la forma de onda que se obtiene con este modelo. Esta forma de onda contiene un número finito (y conmensuradas entre sí) de componentes de Fourier, como se puede constatar en el espectro de potencias mostrado en la Fig. 4(B). El retrato de fase en la Fig. 4(C) muestra claramente las trayectorias heteroclínicas que unen a los dos puntos fijos estables.

Todos los cálculos fueron hechos con el programa COMSOL Multiphysics (COMSOL), que usa el método del elemento finito, y con MATLAB, usando un método de Euler simple. Todos los cálculos concuerdan.

 

Cálculos numéricos en una superficie con curvatura

No es evidente que los patrones obtenidos en un dominio bidimensional plano con condiciones periódicas se obtengan en otras geometrías. Para aplicar nuestro modelo al funcionamiento de un corazón real, al menos es necesario hacer los cálculos en una superficie cerrada convexa. En la Fig. 5 mostramos que las ondas viajeras se presentan en una superficie paramétrica simple que se asemeja a la forma de una cavidad cardiaca. En esta figura hemos asociado una contracción del radio de curvatura local de la superficie proporcional al valor de u(x, y) y la barra de colores señala en amarillo las regiones de contracción. Ésta es una representación burda de la actividad muscular del corazón debida a la actividad eléctrica modelada.

Los solitones que se muestran en la Fig. 5(A) rotan alrededor de los puntos de curvatura cero, es decir, tienen el mismo comportamiento que en un dominio plano. En la Fig. 5(B) se muestran estos solitones pero para un cálculo en el que se obtienen ondas espirales del tipo de las mostradas en la Fig 3, del lado izquierdo de la línea azul.

El tamaño del dominio fue regulado variando η de tal forma que sólo presentará dos centros de rotación, en una cavidad modelo de un "corazón sano".

 

En un dominio tridimensional

Es evidente que el músculo cardiaco es un dominio tridimensional, por lo tanto debemos probar nuestro modelo de la actividad eléctrica del corazón en tres dimensiones. Para esto usamos el programa COMSOL Multiphysics, que es capaz de integrar ecuaciones diferenciales no lineales en un dominio de tres dimensiones usando el método de elemento finito.

En la Fig. 6 se muestra una instantánea de la concentración de υ en un dominio acotado por dos semi-elipses de diferente excentricidad y una cubierta inferior esférica, en un cálculo con los mismos parámetros de la Fig 4. Esto simula la forma de una cavidad cardiaca. La escala de colores es igual que en las demás figuras. Debido a la transparencia de la figura, observamos que el patrón de ondas solitarias sigue un comportamiento igual que en el dominio plano y se transmite en la tercera dimensión, lo cual significa que todos los puntos del interior presentan la misma concentración de υ que en las superficies que los acotan. En un corazón sano se presentandos nodos principales, por lo que, para obtener dos nodos fijos, en este cálculo se fijó el valor de u = 0 en los extremos superior e inferior del dominio, lo cual tiene el efecto de fijar el punto de rotación de los solitones, simulando así la excitación externa del corazón.

 

Comparación con el funcionamiento del corazón

Un dato clave experimental es el electrocardiograma que, en términos generales, simula la historia temporal del potencial de membrana de las células cardiacas. La Fig. 7 muestra un electrocardiograma típico de una persona sana.

Como puede observarse, el retrato de fase es comparable con el mostrado por los solitones, en el sentido de formar una trayectoria cerrada entre dos puntos fijos. Sin embargo, la forma de onda no es similar. Esto se debe obviamente a que en el ECG no se registra la actividad en una única célula. De todas formas podemos ver que la forma de onda obtenida con nuestro modelo es mucho más estructurada que las del modelo de FiztHugh-Nagumo, mostradas en la Fig. 1.

Nuestro modelo es capaz de reproducir historias temporales de excitación similares a las formas de onda registradas en un ECG, si variamos los parámetros adecuadamente. En la Fig. 8(A) mostramos un cálculo que presenta "wavelets", en los que los frentes de u máxima son angostos, la forma de onda es muy similar a la mostrada en la Fig. 7. También mostramos el espectro de Fourier de las oscilaciones.

En la Fig. 8(B) mostramos el patrón de ondas espirales que se obtiene justo antes de cruzar la línea de bifurcación azul de la Fig. 3, que se atraviesa al disminuir g. Observe que la historia temporal en la Fig. 8(B) muestra oscilaciones mucho más rápidas que en la Fig. 8(A), simulando taquicardia. Además, la forma de onda de las ondas espirales es prácticamente senoidal como se puede constatar en el espectro de frecuencias mostrado.

Uno de los puntos importantes que pueden validar este modelo es investigar si es capaz de simular las fallas y anormalidades del funcionamiento de un corazón. Una de las afecciones más frecuentes, que implican la muerte cardiaca, es la llamada fibrilación ventricular, en la que después de una taquicardia, se presentan arritmias severas y finalmente la muerte, ya que en estas condiciones el corazón no puede bombear sangre eficientemente.

En la Fig. 9 mostramos los cambios producidos en los patrones y en el comportamiento temporal de u cuando se cambian los parámetros del modelo.

Empezamos mostrando en (a) los solitones usando los valores de f y g de la Fig. 8, que corresponden a la simulación de un corazón sano. En (b) se muestra la aparición de un doblamiento de período cuando se aumenta el valor de f=f + 0.012 y g = g_c. Observe la aparición de otros nodos. En (c) obtenemos la aparición de "wavelets" rotando a alta frecuencia cuando se baja el valor de f para acercarse al valor f y g_c. En (d) sigue la aproximación a la línea crítica, que se alcanza en el cálculo mostrado en (f ). Observe la sucesión de bifurcaciones al hacer esto. En ese punto se modifica el valor de C = 0.2 y se cumple la condición de caos obtenida en el apéndice (ver Ec. VIII.31) poniendo D_u/D_ν = 0.77. Observe la irregularidad de las oscilaciones y la aparición de múltiples frecuencias.

Este panorama es muy similar al mostrado en la Fig. 5 de Garfinkel et al.^[6], que muestran los electrocardiogramas de corazones humanos y caninos con fibrilación. Es conveniente enfatizar que los retratos de fase mostrados en esa referencia son similares a los obtenidos con nuestro modelo, tanto para un corazón sano, como con fibrilación.

 

Conclusiones

En este trabajo hemos examinado las propiedades de un modelo lo suficientemente rico para ser usado en diversos sistemas complejos. El trabajo fue guiado por un análisis lineal del sistema de ecuaciones, hecho alrededor de un punto fijo diferente al usado para obtener patrones de Turing. La riqueza de patrones que se obtienen es debida a la proximidad de dos líneas críticas en el espacio de parámetros. Esto permite tener una región de transición entre una bifuración de Hopf y otra de silla-nodo, la primera exhibe oscilaciones regulares en un ciclo límite, y la segunda presenta una situación biestable, que puede presentar frentes de onda viajeras del tipo solitón.

Todas estas características hacen de este modelo un sistema apropiado para simular la actividad eléctrica en el corazón, ya que en este órgano se presentan situaciones muy peculiares en el momento en que falla. Si esto es así, podríamos proponer que las ondas observadas en el músculo cardiaco son en realidad solitones, lo cual tiene sentido a posteriori, si se toma en cuenta que las oscilaciones tienen que ser lo suficientemente robustas para que el corazón lata sin fallar más de 10⁹ veces. Las ondas espirales no tienen esta robustez ni producen una forma de onda comparable a la medida en los ECG. Además, los solitones encontrados en nuestro modelo se pueden comportar de distinta forma al chocar con un obstáculo o una pared, dependiendo de las condiciones a la frontera. Pueden rebotar, como es usual, o anularse en la pared, o bien pueden proseguir el estado de excitación del bulto en la pared^[24], lo cual se observa en el corazón.

Si consideramos la forma de onda, podemos ver que la forma de onda del comportamiento colectivo en el dominio es comparable a la registrada en un ECG, mientras que otros modelos no lo hacen. Además, es de vital importancia que la frecuencia cardiaca pueda ser regulada en forma directa y sencilla. En nuestro modelo esto se puede hacer de dos maneras, disminuyendo f, o aumentando g. Si interpretamos a los morfógenos υ y ν como las substancias cuya reacción provoca el latido cardiaco, por ejemplo en la fosforilación, el tiempo de reacción debe depender de otra substancia, por ejemplo la adrenalina, que modifica los valores de los parámetros externos, los cuales tienen precisamente este significado.

Hemos mostrado además que las ecuaciones de amplitud del modelo BVAM pueden transformarse en un sistema de Landau-Ginzburg, el cual es susceptible de presentar caos bajo ciertas condiciones de los parámetros (ver Apéndice). Esto es de suma importancia, ya que trabajos experimentales^[6,29] muestran que cuando el corazón muere, lo hace siguiendo una ruta específica al caos, que se manifiesta como taquicardia, arritmia y fibrilación. La sugerencia experimental es que la ruta al caos sea siguiendo el mecanismo de Ruelle-Takens-Newhouse^[27]. En un artículo anterior hemos mostrado que el modelo BVAM presenta patrones de Turing con caos temporal siguiendo precisamente esa ruta^[25].

En este trabajo sugerimos que las fallas cardiacas pueden ser modeladas encontrando una ruta al caos similar a la observada experimentalmente. Es alentador que el tipo de fenómenos que el modelo BVAM presenta cerca de los estados críticos se asemeje mucho a los datos experimentales. Sin embargo, una demostración de la aparición de caos y una comparación detallada con los datos experimentales está actualmente en proceso y será el objeto de un trabajo futuro.

 

Agradecimientos

Este trabajo fue parcialmente financiado por los proyectos de CONACYT Nos. 179616 y 167244. Agradecemos a Gertrudis Hortensia González Gómez por sus excelentes comentarios en la revisión del manuscrito.

 

Referencias

1. Barrio, R.A., Varea, C., Aragón, J.L. & Maini, P.K. A two dimensional numerical study of spatial pattern formation in interactiong turing systems. Bull. Math. Biol. 61, 483 (1999).         [ Links ]

2. Barrio, R.A. Morphogenesis in Introducción a la Física Biológica. Vol. 3. Colín, L.G., Dagdug, L., Picquart, M. & Vázquez, E. (El Colegio Nacional, México, D.F., 2010).         [ Links ]

3. Bub, G., Shrier, A. & Glass, L. Spiral wave generation in heterogeneous excitable media. Phys. Rev. Lett. 85, 058101 (2002).         [ Links ]

4. Barrio, R.A. Aplicaciones del modelo bvam a sistemas complejos. Revista Digital Universitaria,http://www.revista.unam.mx/vol.ll/num6/art55, 11(6) (2010).         [ Links ]

5. Sundnes, J. et al. in Computing the Electrical activity of the Heart, Monographs in Computational Science and Engineering (eds. Bart, T.J. et al.) (Springer, Berlin, 2006).         [ Links ]

6. Garfinkel, A. et al. Quasiperiodicity and chaos in cardiac fibrillation. The Journal of Clinical Investigation 99(2) (1997) . Págs. 305-314.         [ Links ]

7. FitzHugh, R.A. Impuses and physiological states in theoretical models of nerve membrane. Biophysical Journal 1, 465 (1961). Págs. 445-466.         [ Links ]

8. Nagumo, J., Arimoto, S. & Yoshizawa, S. An active pulse transmission line simulating nerve axon. Proceedings of the IRE 50(10), 2061-2070 (1962).         [ Links ]

9. Echebarría, B. & Karma, A. Instability and spatiotemporal dynamics of alternants in paced cardiac tissue. Phys. Rev. Lett. 88(20), 208101 (2002).         [ Links ]

10. Hodgkin, A.L. & Huxley, A.F. A quantitative description of membrane current and its application to conduction and excitation in nerve. Journal of Physiology 117, 500-544 (1952).         [ Links ]

11. Beeler, G.W. & Reuter, H. Reconstruction of the action potential of ventricular myocardial fibres. J. Physiol. 268, 177-210, 1975).         [ Links ]

12. Luo, C.H. & Rudy, Y. A dynamic model of the cardiac action potential. ii. Afterdepolarizations, triggered activity, and potentiation. Circ. Res. 74, 1097-113 (1994).         [ Links ]

13. Rogers, J.M. & McCullock, A.D. A collocation-galerkin fe of cardiac action potential propagation. IEEE Trans. on Biomedical Engineering 41, 743-757 (1994).         [ Links ]

14. Noble, D. A modification of the hodgkin-huxley equations applicable to Purkinje fibre action and pace-maker potentials. J. Physiol. 160, 317-52 (1962).         [ Links ]

15. McAllister, R.E., Noble, D. &Tsien, R.W. Reconstruction of the electrical activity of cardiac Purkinje fibres. J. Physiol. 251, 1-59 (1975).         [ Links ]

16. Luo, C.H. & Rudy, Y. A dynamic model of the cardiac ventricular action potential. i. Simulations of ionic currents and concentration changes. Circ. Res. 74, 1071-96 (1994).         [ Links ]

17. Nygren, A. et al. Mathematical model of an adult human atrial cell: the role of k+ cu- rrents in repolarization. Circ. Res. 82, 63-81 (1998).         [ Links ]

18. Courtemanche, M., Ramírez, R.J. & Nattel, S. Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am. J. Physiol. 275, 301-321 (1998) .         [ Links ]

19. Courtemanche, M., Ramírez, R.J. & Nattel, S. Ionic targets for drug therapy and atrial fibrillation-induced electrical remodeling: insights from a mathematical model. Cardiovas. Res. 42, 477-489 (1999).         [ Links ]

20. Turing, A.M. The chemical basis of morphogenesis. Phil.Trans. R. Soc. London B 237(641), 37-72 (1952).         [ Links ]

21. Murray, J.D. Mathematical Biology II: Spatial Models and Biomedical Applications. Vol. 2. Interdisciplinary Applied Mathematics. 3^rd ed. (Springer, 2003).         [ Links ]

22. Lepannen, T. Computational Studies of Pattern Formation in Turing Systems. Ph.D. Thesis, Helsinki University of Technology (2004).         [ Links ]

23. Wolley, T.E. et al. Analysis of stationary droplets in a generic turing reaction-diffusion system. Phys. Rev. E. 82, 051929 (2010).         [ Links ]

24. Varea, C., Barrio, R.A. & Hernández, D. Soliton behaviour in reaction-diffusion model. J. Math. Biol. 54, 797-813 (2007).         [ Links ]

25. Aragón, J.L., Barrio, R.A., Woolley, T.E., Baker, R.E. & Maini, P.K. Non-linear effects on turing patterns: Time oscillations and chaos. Phys. Rev. E., 86, 026201 (2012).         [ Links ]

26. Strogatz, S.H. Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering. Studies in Nonlinearity (Westview Press, Cambridge, MA, 2008).         [ Links ]

27. Newhouse, S., Ruelle, D. & Takens, F. Occurrence of strange axiom a attractors near quasi periodic flows on tm, m > 3. Commun. Math. Phys. 64, 35-40 (1978).         [ Links ]

28. Aragón, J.L., Torres, M., Gil, D., Barrio, R.A. & Maini, P.K. Turing patterns with pentagonal symmetry. Phys. Rev. E. 65, 051913 (2002).         [ Links ]

29. Herlin, A. & Jacquemet, V. Eikonal-base initiation of fribillatory activity in thin walled cardiac propagation models. Chaos 21, 043136 (2011).         [ Links ]

30. Kuramoto, Y. Chemical Oscillations, Waves, and Turbulence (Springer, Berlin, NewYork, 1984).         [ Links ]

 

Información sobre los autores

Ivonne Domínguez

Ivonne Domínguez cursó sus estudios de Licenciatura en Física en la Facultad de Ciencias de la UNAM y durante un semestre en la Fakultät für Physics de la Universität Wien, Viena, Austria.

Actualmente se desempeña como asistente de investigador en el Instituto de Física de la UNAM donde realiza su tesis de licenciatura. Su línea de investigación está dirigida a la Dinámica no lineal y Sistemas complejos. Ésta se centra en el estudio de la fibrilación cardiaca usando un modelo de reacción-difusión. Ha tomado varios talleres nacionales de investigación en Nanotecnología, Energías renovables, Matemáticas y Matemáticas aplicadas. Recientemente participó en la escuela internacional "Spring College on Physics of Complex Systems" organizada por el ICTP (International Centre for Theoretical Physics) en Trieste, Italia.

Responsable de varios proyectos DGPA, PRIDE y CONACYT. Se le otorgó el Premio SMF a la Investigación Científica.

Rafael A. Barrio

Rafael A. Barrio obtuvo su licenciatura y dos maestrías en la UNAM, su doctorado en Física Teórica en la Universidad de Oxford UK y ha trabajado en la UNAM por 39 años, primero en el Instituto de Investigaciones en Materiales y actualmente en el Instituto de Física, como Investigador Titular "C" de tiempo completo. Es Profesor Honorario (Adjunt Professor) del Centre of Excellence in Computational Complex Systems Research de la Universidad Aalto, Helsinki, Finlandia. Actualmente coordina el proyecto Patterning, Segregation and Differentiation in Complex Networks, en el que participan la UNAM, la Universidad de Barcelona, la Universidad de Oxford y la Universidad de Helsinki. Es nivel III del SNI. Ha publicado 171 artículos a nivel internacional, entre 15 libros y capítulos de libros; asesorado 19 tesis. Ha presentado 86 conferencias internacionales y 122 pláticas por invitación, también a nivel internacional.

Su campo de especialidad es muy amplio en Metalurgia, Estado sólido y Sistemas complejos. Ha sido miembro de varios comités. Responsable de varios proyectos DGPA, PRIDE Y CONACYT. Se le otorgó el Premio SMF a la Investigación Científica.

Carmen varea

Carmen Varea es Doctora en Ciencias por la Universidad de California en Berkeley. Actualmente, trabaja como Investigadora Titular "C" en el Instituto de Física de la UNAM. Sus principales líneas de investigación son: la mecánica estadística de sistemas inhomogéneos, la formación de patrones, las cinéticas en sistemas físicos, químicos, biológicos y humanos, y la física del estado sólido. Es miembro del SNI nivel III. Ha colaborado con artículos en diversas publicaciones internacionales y recibió, en 1996, el Premio de investigación de la Sociedad Mexicana de Física.

José Luis Aragón

José Luis Aragón es Ingeniero Físico egresado de la Universidad Autónoma Metropolitana, doctor en física de materiales por el CICESE e Investigador Titular "C" del Centro de Física Aplicada y Tecnología Avanzada de la UNAM. Su línea de investigación principal es Biomatemáticas y sus áreas de interés son Física del estado sólido y Biología matemática. Ha publicado 99 artículos, la mayoría en revistas indexadas; también capítulos de libros y artículos in extenso en memorias con refereo, así como artículos periodísticos y de divulgación. Otras publicaciones son 26 libros de texto, un libro traducido y 2 números de libros especiales.

Ha dirigido 5 tesis de licenciatura, 1 de maestría y 3 de doctorado. Ha impartido 57 cursos de licenciatura y 7 de doctorado. Ha sido Jefe del Depto. de Nanotecnología del Centro de Física Aplicada y Tecnología Avanzada y Coordinador de Licenciatura en Tecnología. Es miembro del SNI nivel III.