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

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.

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.

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 × 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 fbras 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 laFigura 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.

Figura 1.

Gráficas de las variables υ (sólido) y ω (punteado) calculadas con las ecuaciones 1 en el artículo original (líneas gruesas). También mostramos el resultado de las modificaciones sugeridas en Rogers & McCullock^[13] (líneas delgadas).

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 u y υ son los morfógenos con coeficientes de difusión D_u 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 u y v 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 (u, υ), el modelo BVAM tiene los tres puntos fijos mostrados en la ecuaciónVIII.15. Si h = −1 todos ellos colapsan en el punto (u₀, υ₀) = (0, 0). En el caso h ≠ 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.

Figura 2.

Diagrama de bifurcación del sistema BVAM, sin difusión, en el plano (g, f), para h = −2.5, C = 0 y η = 0.5. La línea continua es una bifurcación de Hopf y la línea punteada es una bifurcación de silla-nodo. Al cruzar la línea, el sistema pasa de tener un ciclo límite en la región izquierda a ser biestable en la región derecha.

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.

Figura 3.

Ejemplos de patrones numéricos en puntos cercanos a la bifurcación de silla-nodo. Observe que los frentes de onda en la región verde viajan, ya sea con velocidad proporcional a la curvatura, o en modo desordenado. El patrón cambia de forma en este último caso, por lo cual se puede esperar la presencia de oscilaciones caóticas como en el punto rojo. En todos los cálculos se usaron los parámetros de la Fig. 2 con D_u = Dυ = 1, excepto en el punto rojo, donde se usó C = 0.2 y D_u = 0.77Dυ. Estos valores se propusieron sin que mediara análisis alguno, sólo para observar el comportamiento cuando el diagrama de bifurcación (calculado para C = 0) ya no es válido.

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 × 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 u 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.

Figura 4.

Resultados de un cálculo usando los parámetros D_u = Dυ = 1, f = f_c + 0.01, g = 0.175, C = 0, h = -2.5, dt = 0.06, dx = 0.2 y η = 10, en un dominio plano discreto con una matriz de 120 × 120, con condiciones periódicas a la frontera. Estos valores de los parámetros nos situán cerca de la línea crítica punteada de la Fig. 2, o sea, en el punto marcado “solitones regulares” en la Fig. 3. (A) Historia en el tiempo de la concentración de u en el sitio central del dominio espacial. (B) Transformada de Fourier de u(t); nótese la periodicidad y los armónicos que nos dan la forma de onda. (C) Retrato de fase, u(t) - v(t), en donde son evidentes las trayectorias heteroclínicas que comunican 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.

Figura 5.

(A) Propagación de ondas tipo solitón en una superficie paramétrica cerrada. (B) Ondas espirales en la misma superficie.

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 u 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 u que en las superficies que los acotan. En un corazón sano se presentan dos 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.

Figura 6.

Ondas solitarias producidas por el modelo BVAM en un dominio tridimensional. Se muestra una instantánea del cálculo en dos perspectivas para que se aprecie la forma del dominio.

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.

Figura 7.

(A) Electrocardiograma típico de una persona sana. (B) Retrato de fase extraído de (A). Compare este retrato de fase con el de la Fig. 4C, en el que se muestra la curva cerrada.

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.

Figura 8.

(A) Historia temporal de u y patrón de u obtenido con los parámetros de la Fig. 4 pero modificando f = f_c− 0.0001 y g = 0.175. (B) Lo mismo, pero ahora modificando g = 0.165. En las gráficas de la columna central se muestran las transformadas de Fourier de los patrones.

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.

Figura 9.

Cálculos numéricos variando los parámetros del modelo para imitar el latido de un corazón afectado por anormalidades. Mostramos la historia temporal de u y una instantánea de la concentración de u en el dominio cuadrado.

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_c + 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_c 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 u y v 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.