Excitabilidad y recuperación en sistemas lentos y rápidos

Dos escalas temporales en un mismo estado

El capítulo 26 presentó por primera vez ecuaciones con escalas temporales distintas. Allí el interés era numérico, porque un modo rápido podía obligar a un método explícito a usar pasos pequeños aunque la evolución visible pareciera lenta. También señalamos que cambiar el integrador y reducir el modelo eran operaciones diferentes. En este capítulo retomamos esa estructura desde la dinámica no lineal.

Cuando la separación de escalas interactúa con una geometría no lineal, una trayectoria puede avanzar lentamente durante un intervalo prolongado, realizar después una transición rápida y regresar a otra región lenta. El valor que alcanza la variable lenta modifica, a su vez, la respuesta posible ante una nueva perturbación. De esta combinación surgen la excitabilidad, la recuperación y las oscilaciones de relajación que el capítulo 42 encontró en el oscilador de Van der Pol con $\mu$ grande.

Consideremos el sistema

La variable $x$ cambia a una velocidad de orden uno en el tiempo $t$, mientras que $y$ lo hace a una velocidad de orden $\varepsilon$. Con $\varepsilon=0{,}01$, una variación de orden uno en $y$ requiere típicamente un intervalo unas cien veces mayor que una variación comparable de $x$, aunque la razón efectiva depende también de $f$ y $g$. La separación de escalas es una propiedad dinámica. Dos variables pueden tener unidades distintas y evolucionar a velocidades parecidas, o tener las mismas unidades y tiempos característicos muy diferentes. Lo que importa es que, después de una elección adecuada de escalas, aparezca un parámetro pequeño que multiplica una de las velocidades.

Con el tiempo lento $\tau=\varepsilon t$ el mismo sistema se escribe

Las dos formas describen el mismo sistema mientras $\varepsilon>0$. La diferencia aparece en el límite singular $\varepsilon\to0$. En tiempo rápido, $y$ queda congelada. En tiempo lento, la primera ecuación deja de ser diferencial y se convierte en una restricción algebraica. La teoría geométrica de perturbaciones singulares aprovecha esa dualidad para analizar el límite y relacionarlo, bajo hipótesis adicionales, con el sistema de $\varepsilon$ positivo (Kuehn, 2015).

El subsistema rápido y la variedad crítica

En tiempo rápido y con $\varepsilon=0$ se obtiene $\dot x=f(x,y;\mu)$ y $\dot y=0$. La variable $y$ se comporta como un parámetro y los equilibrios del subsistema rápido forman la variedad crítica,

El adjetivo crítica significa aquí que el conjunto está formado por equilibrios de la dinámica rápida, algunos atractores y otros repulsores. Con una sola variable rápida, la estabilidad transversal se lee en $f_x$. Donde $f_x<0$ la variedad atrae en la dirección rápida, donde $f_x>0$ repele y donde $f_x=0$ la linealización rápida pierde hiperbolicidad.

La condición $f_x\neq0$ expresa en este caso escalar la hiperbolicidad normal. El teorema de Fenichel establece que una porción compacta y normalmente hiperbólica de $C_0$ persiste, para $\varepsilon$ positivo y pequeño, como una variedad lenta invariante cercana que conserva su carácter atractor o repulsor (Fenichel, 1979). El resultado es local en torno a la porción considerada y excluye justamente los puntos donde se pierde la hiperbolicidad normal.

El problema reducido

En tiempo lento y con $\varepsilon=0$ quedan la restricción $f(x,y;\mu)=0$ y la ecuación $dy/d\tau=g(x,y;\mu)$. La restricción obliga a permanecer sobre $C_0$ y la ecuación diferencial describe el desplazamiento a lo largo de ella. Si la variedad puede escribirse localmente como $y=h(x)$, la derivación respecto de $\tau$ da

Siempre que $h'(x)\neq0$ se obtiene una ecuación para el movimiento lento. En un pliegue, donde $h'(x)=0$ y se pierde la hiperbolicidad normal, la reducción se vuelve singular. La aproximación lenta falla allí porque la propia descripción asintótica cambia de carácter. Cerca de un pliegue, términos que parecían despreciables deciden si la trayectoria salta hacia otra rama, permanece un tiempo cerca de una región repulsora o atraviesa una transición más delicada.

El modelo de FitzHugh-Nagumo

Usaremos una versión adimensional del modelo propuesto por FitzHugh (1961) y llevado a un circuito eléctrico por Nagumo y colaboradores (1962),

con $a=0{,}7$, $b=0{,}8$ y $\varepsilon=0{,}08$. No lo trataremos como la descripción literal de una neurona. Nos interesa porque contiene una variable rápida, una variable lenta, una nulclina cúbica con pliegues y un régimen en el que un equilibrio estable responde a perturbaciones grandes con una excursión amplia antes de regresar.

Con $I=0$ la nulclina rápida es $y=x-x^3/3$ y la lenta, $y=(x+a)/b$. El equilibrio está en su intersección, que cumple $bx^3+(3-3b)x+3a=0$, es decir, $0{,}8x^3+0{,}6x+2{,}1=0$. Su derivada, $2{,}4x^2+0{,}6$, es positiva, de modo que la cúbica es estrictamente creciente y tiene una única raíz real,

El jacobiano es

con traza $-0{,}502580$ y determinante $0{,}108069$ en el equilibrio. Los valores propios son $-0{,}251290\pm0{,}211949\,i$, por lo que el equilibrio es un foco estable. La excitabilidad que estudiaremos no requiere, por tanto, que el reposo sea inestable. En el equilibrio, la tasa rápida $1-x^{*2}\approx-0{,}439$ es unas 6,9 veces mayor en valor absoluto que la tasa lenta $-\varepsilon b=-0{,}064$. Con $\varepsilon=0{,}08$ la separación de escalas es clara, aunque no extrema.

Pliegues y ramas de la cúbica

Sobre la variedad crítica $y=x-x^3/3$ se tiene $f_x=1-x^2$. Las ramas exteriores, con $|x|>1$, atraen en la dirección rápida y la rama central, con $|x|<1$, repele. Los pliegues están en $x=\pm1$, con $y=\pm2/3$ (figura 47.1). Lejos de ellos, una trayectoria cercana a una rama atractora se acerca rápidamente a una variedad lenta y luego se desplaza sobre ella mucho más despacio. Al llegar a una zona donde la atracción transversal desaparece, el movimiento lento deja de sostenerse y se produce una transición rápida hacia otra rama.

El límite singular sugiere dónde está el umbral. Si una perturbación instantánea desplaza $x$ con $y$ congelada, el subsistema rápido la devuelve a la rama izquierda cuando el nuevo valor de $x$ queda a la izquierda de la rama central. Si queda a su derecha, la lleva a la rama derecha. Con $y=y^*$, la rama central pasa por $x\approx-0{,}786$, por lo que el límite singular predice un umbral $\Delta x\approx0{,}413$. Con $\varepsilon=0{,}08$ la integración numérica da $0{,}555$, porque la variable lenta se mueve durante la respuesta. El umbral numérico se acerca a la predicción cuando $\varepsilon$ disminuye y vale $0{,}459$, $0{,}426$ y $0{,}416$ con $\varepsilon=0{,}02$, $0{,}005$ y $0{,}001$.

La regla de llegar al pliegue y saltar tampoco es exacta. Para $\varepsilon>0$ la variedad lenta no coincide con $C_0$. Cerca de la pérdida de hiperbolicidad aparecen, además, retrasos y órbitas especiales, llamadas canards, que siguen durante un tiempo la rama repulsora. El límite singular organiza la geometría, aunque lo que ocurre cerca de los pliegues exige analizar el sistema con $\varepsilon$ positivo.

Excitabilidad y umbral de un protocolo

Un sistema excitable tiene un estado de reposo estable, las perturbaciones pequeñas decaen directamente hacia él y algunas perturbaciones suficientemente grandes desencadenan una excursión macroscópica antes del retorno. En los experimentos siguientes, los parámetros $a$, $b$ y $\varepsilon$ permanecen fijos y solo cambia la entrada $I(t)$. Las trayectorias diferentes provienen de una misma organización dinámica, sin que el sistema atraviese ninguna bifurcación.

Apliquemos, desde el equilibrio, un pulso $I(t)=A$ entre $t=5$ y $t=6$, con $I=0$ fuera de ese intervalo. Con $A=0{,}60$ la variable rápida aumenta hasta $-0{,}461$ y la perturbación decae sin excursión. Con $A=0{,}65$ la trayectoria alcanza $x\approx1{,}708$, recorre una amplia vuelta en el plano de fases y regresa al mismo equilibrio (figura 47.2a). Con el criterio operativo $\max x>1$, el cambio de respuesta ocurre en $A\approx0{,}60653$ para este protocolo.

Ese valor no es una constante material del modelo. Si cambian la duración del pulso, el estado inicial, la dirección de la perturbación o el criterio con que se define una excursión grande, cambia también el umbral. La pregunta que responde el cálculo es qué amplitud de un pulso de duración uno hace que $\max x>1$. Hablar del umbral intrínseco del sistema requiere especificar antes qué perturbaciones están permitidas y qué respuesta se considera grande. Por eso el umbral medido debe informarse junto con el protocolo que lo produjo.

Memoria en la variable lenta y período refractario

La excursión no devuelve de inmediato al sistema a su condición inicial. Mientras $x$ cambia rápidamente, $y$ también evoluciona. Al terminar la fase activa, la variable lenta puede estar lejos de $y^*$. Dos instantes con valores parecidos de $x$ no son, por tanto, estados equivalentes. Esta es una forma elemental de memoria dinámica, en la que el estado completo conserva una consecuencia lenta del recorrido sin que ninguna variable almacene explícitamente el pasado.

Un desplazamiento instantáneo $\Delta x=1$ aplicado en reposo produce una excursión con máximo $x\approx1{,}793$. Si durante la recuperación se aplica el mismo desplazamiento por segunda vez en $t=24$, no hay una segunda excursión y el máximo posterior es $-0{,}632$. En $t=32$ el mismo desplazamiento produce de nuevo una respuesta grande, con máximo $1{,}712$. En $t=50$ la respuesta, $1{,}794$, es prácticamente la de reposo.

La geometría del límite singular explica estos resultados. En $t=24$ el estado es $(x,y)\approx(-1{,}632;\,-0{,}221)$. Tras el desplazamiento $x$ pasa a $-0{,}632$, que está a la izquierda de la rama central en esa altura, $x\approx-0{,}225$, así que la dinámica rápida devuelve el estado a la rama izquierda. En $t=32$ el estado es $(-1{,}382;\,-0{,}527)$. El desplazamiento lleva $x$ a $-0{,}382$, a la derecha de la rama central, que ahora está en $x\approx-0{,}599$. Se produce entonces la excursión. El impulso aplicado no cambió. Cambió la altura $y$ del estado sobre el que actúa y, con ella, la posición de la rama central, que hace de umbral. La figura 47.2b muestra la respuesta en función del instante del segundo impulso. Hay un período refractario absoluto, durante el cual no se produce ninguna excursión, que termina cerca de $t=30$, seguido de un período en que la respuesta es algo menor que la de reposo.

Recuperación del estado completo

La palabra recuperación puede sugerir un temporizador externo. El modelo muestra una interpretación más precisa. El sistema se recupera porque su estado vuelve gradualmente a una región desde la cual una perturbación dada puede producir de nuevo una excursión amplia. En la trayectoria estudiada, $x$ regresa a valores negativos y pasa incluso por debajo de su valor de reposo mucho antes de que $y$ vuelva al equilibrio (figura 47.3a). El observable rápido parece haber terminado cuando el sistema completo todavía no se ha restablecido.

Que la variable rápida recupere su valor basal no implica que la lenta haya recuperado el suyo. Aun cuando ambas estén cerca del equilibrio, nada garantiza que la respuesta a una nueva perturbación sea la misma que en reposo. En un sistema no lineal, una diferencia pequeña en la variable lenta puede modificar mucho la respuesta a una entrada posterior.

Una excursión recorrida por fases

La excursión completa puede reconstruirse desde la geometría. El equilibrio está sobre la rama izquierda, que atrae en la dirección rápida. Una perturbación débil desplaza $x$, pero el campo rápido la devuelve al entorno de esa rama casi sin que cambie $y$. Una perturbación mayor lleva el estado más allá de la rama central, donde la dinámica rápida es repulsora. Desde ahí $x$ cambia en un intervalo corto comparado con la evolución de $y$ y la trayectoria cruza el plano casi en horizontal, porque $y$ se comporta como un parámetro congelado.

Al llegar a la rama derecha vuelve a dominar la separación de escalas y el estado se desplaza lentamente a lo largo de ella mientras $y$ aumenta, porque la ecuación lenta sigue actuando. Cerca del pliegue derecho la rama deja de sostener la atracción normal y se produce una transición rápida hacia la izquierda. La trayectoria llega a la rama izquierda con un valor de $y$ muy distinto del inicial y la recuperación consiste en una deriva lenta a lo largo de esa rama hasta que $y$ se acerca a $y^*$. El cuadro 47.1 resume el recorrido según la escala que domina en cada tramo. Las fronteras entre fases dependen de la tolerancia con que se decida que un movimiento es rápido o lento, por lo que no son estados discretos del sistema.

Cuadro 47.1. Fases de una excursión del modelo de FitzHugh-Nagumo según la escala temporal dominante.

Oscilaciones de relajación y período singular

Si la dinámica lenta conduce a cada rama hacia un pliegue sin detenerse en un equilibrio, la trayectoria repite los saltos y cierra una órbita periódica, una oscilación de relajación (Strogatz, 2015, sec. 7.5). Con $I=0{,}5$ el equilibrio de FitzHugh-Nagumo está en $x\approx-0{,}805$, sobre la rama central repulsora. Ninguna de las ramas exteriores contiene equilibrios. El problema reducido permite calcular el período en el límite singular. Con $I$ constante, la variedad crítica es $y=x-x^3/3+I$ y la ecuación reducida es

El denominador explica por qué $x=\pm1$ exige cuidado. La velocidad reducida diverge allí, aunque el sistema original no tiene velocidad infinita. Lo que diverge es la aproximación construida al imponer que la trayectoria permanezca sobre la variedad crítica. Una singularidad de la reducción no es, sin más, una singularidad del modelo.

En el límite singular, el ciclo recorre la rama izquierda desde $x=-2$ hasta el pliegue $x=-1$, salta horizontalmente hasta $x=2$ sobre la rama derecha, la recorre hasta el pliegue $x=1$ y salta de vuelta a $x=-2$. Los puntos de llegada provienen de igualar las alturas de los pliegues, porque $x-x^3/3=\pm2/3$ tiene las raíces dobles $\pm1$ y las simples $\mp2$. Como los saltos no consumen tiempo lento, el período en tiempo lento es

y el período en tiempo rápido es aproximadamente $T_0/\varepsilon$. El bloque siguiente evalúa $T_0$ por cuadratura y lo compara con el producto $\varepsilon T$ medido en simulaciones. Usa el método LSODA de SciPy (Virtanen et al., 2020), que cambia a un método implícito cuando detecta rigidez.

import numpy as np
from scipy.integrate import quad, solve_ivp

a, b, I = 0.7, 0.8, 0.5
N = lambda x: x + a - b * (x - x**3 / 3 + I)    # velocidad lenta sobre la cúbica
T0 = (quad(lambda x: (1 - x**2) / N(x), -2, -1)[0]
      + quad(lambda x: (1 - x**2) / N(x), 2, 1)[0])
print(f"período singular en tiempo lento  T0 = {T0:.5f}")

cruce = lambda t, z: z[0]                       # x = 0 recorrido hacia arriba
cruce.direction = 1
print("   eps    período T    eps*T")
for eps in (0.08, 0.02, 0.005):
    campo = lambda t, z: [z[0] - z[0]**3 / 3 - z[1] + I,
                          eps * (z[0] + a - b * z[1])]
    sol = solve_ivp(campo, (0, 24 / eps + 100), [2.0, 0.0], events=cruce,
                    method="LSODA", rtol=1e-10, atol=1e-12)
    T = np.diff(sol.t_events[0])[-1]
    print(f"{eps:6.3f}   {T:10.4f}   {eps * T:7.4f}")

Al ejecutarlo se obtiene esta salida.

período singular en tiempo lento  T0 = 1.95196
   eps    período T    eps*T
 0.080      39.4744    3.1580
 0.020     122.1134    2.4423
 0.005     430.7446    2.1537

El producto $\varepsilon T$ se acerca a $T_0\approx1{,}952$, aunque despacio. La diferencia $\varepsilon T-T_0$ vale $1{,}21$, $0{,}49$ y $0{,}20$. Dividida por $\varepsilon^{2/3}$ da números entre 6,5 y 6,9. La corrección decrece, entonces, aproximadamente como $\varepsilon^{2/3}$, la misma escala del paso lento por un pliegue que el capítulo 43 obtuvo a partir de la ecuación de Airy. El tiempo que la trayectoria pasa cerca de cada pliegue es el origen de la corrección. La figura 47.3b muestra que con $\varepsilon=0{,}005$ el ciclo ya sigue de cerca las ramas de la cúbica y los saltos horizontales del límite singular.

Excitabilidad, biestabilidad y oscilación

Un sistema excitable puede tener un único atractor estable. En el ejemplo con $I=0$, todas las trayectorias estudiadas vuelven al mismo equilibrio. La excitabilidad se distingue así de la biestabilidad, en la que dos atractores coexisten y una perturbación traslada el estado de una cuenca a otra. También se distingue de una oscilación autosostenida, en la que el sistema recorre un ciclo sin necesidad de nuevas perturbaciones, como ocurre con $I=0{,}5$.

La misma nulclina cúbica participa en los tres fenómenos sin volverlos equivalentes. Lo que decide el régimen es la posición del equilibrio respecto de los pliegues y el resultado de las bifurcaciones que la reorganizan al cambiar los parámetros, algo que debe demostrarse con el análisis correspondiente. Un pico transitorio grande, un salto entre atractores y un ciclo repetitivo pueden producir registros parecidos durante un intervalo corto, aunque su organización en el espacio de estados sea distinta.

Del límite singular al sistema con ε positivo

El límite $\varepsilon=0$ da una geometría poderosa que no debe leerse como una igualdad exacta para $\varepsilon>0$. El teorema de Fenichel conecta ambos en las regiones normalmente hiperbólicas, donde las variedades lentas se mantienen cerca de las porciones correspondientes de la variedad crítica (Kuehn, 2015). En los pliegues desaparece la separación hiperbólica transversal que sostenía el argumento. Allí se concentra parte de la dinámica más interesante de estos sistemas.

Tampoco hay un significado universal de $\varepsilon$ pequeño. La validez de una aproximación depende del intervalo temporal considerado y de la distancia a las regiones no hiperbólicas. Con $\varepsilon=0{,}08$ la separación de escalas es visible, aunque el umbral de excitación difiere en un 35 % de su valor singular y el período de la oscilación de relajación, en un 62 %.

Las simulaciones del capítulo usaron el método DOP853 (Hairer et al., 1993) con tolerancia relativa $10^{-10}$, tolerancia absoluta $10^{-12}$ y paso máximo $0{,}02$. Para el pulso de amplitud $0{,}65$, el máximo de $x$ obtenido con el método implícito Radau (Hairer y Wanner, 1996) difiere en $3\times10^{-6}$ del calculado con DOP853. El equilibrio y la posición de los pliegues se obtuvieron analíticamente y no dependen de esos cálculos.

El capítulo 48 introducirá otra forma de dependencia temporal. Hasta aquí, el estado presente $(x,y)$ contiene toda la información necesaria para determinar el futuro. En los sistemas con retardo, la derivada depende explícitamente de valores anteriores de la trayectoria. Esa modificación obliga a ampliar la idea de estado.

Ejercicios

Ejercicio 47.1. En el modelo de FitzHugh-Nagumo con $I=0$, (a) calcule en el equilibrio la tasa rápida $f_x$, la tasa lenta $-\varepsilon b$ y su cociente. (b) ¿Cómo cambia ese cociente si $\varepsilon=0{,}01$? (c) Explique por qué el cociente mide mejor la separación de escalas que el valor de $\varepsilon$ por sí solo.

Ejercicio 47.2. Considere $I$ constante con $a=0{,}7$, $b=0{,}8$ y $\varepsilon=0{,}08$. (a) ¿Para qué valor de $I$ el equilibrio coincide con el pliegue izquierdo $x=-1$? (b) La bifurcación de Hopf ocurre cuando la traza del jacobiano se anula. Calcule el valor de $x$ correspondiente y el valor de $I$. (c) Explique por qué los dos valores difieren y qué ocurre con esa diferencia cuando $\varepsilon\to0$.

Ejercicio 47.3. Con $I=0$, el problema reducido sobre la rama izquierda es $dx/d\tau=N(x)/(1-x^2)$, con $N(x)=x+a-b(x-x^3/3)$. (a) Compruebe que $N(x^*)=0$. (b) Calcule la derivada del segundo miembro en $x^*$ y conviértala a tiempo rápido. (c) Compárela con la parte real de los valores propios del sistema completo.

Ejercicio 47.4. Para la oscilación de relajación con $I=0{,}5$, (a) calcule por separado el tiempo lento que se pasa en cada rama si las integrales valen $1{,}219$ y $0{,}733$. (b) Estime el período en tiempo rápido con $\varepsilon=0{,}02$ usando solo el límite singular y compárelo con el valor numérico, $122{,}1$. (c) Explique el signo de la diferencia.

Ejercicio 47.5. Use el criterio geométrico del texto para decidir si un desplazamiento $\Delta x=1$ produce una excursión en los estados $(x,y)=(-1{,}632;\,-0{,}221)$ y $(-1{,}382;\,-0{,}527)$. Para cada estado, calcule la posición de la rama central a la altura $y$ correspondiente.

Ejercicio 47.6. Clasifique cada comportamiento como excitabilidad, biestabilidad u oscilación autosostenida y justifique la clasificación. (a) Un registro muestra picos aislados de gran amplitud que aparecen solo después de estímulos externos. (b) Un sistema permanece indefinidamente en uno de dos niveles y cambia de nivel tras perturbaciones grandes. (c) Un registro muestra picos regulares sin estímulos externos.

Referencias

Fenichel, N. (1979). Geometric singular perturbation theory for ordinary differential equations. Journal of Differential Equations, 31(1), 53–98. https://doi.org/10.1016/0022-0396(79)90152-9

FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal, 1(6), 445–466. https://doi.org/10.1016/S0006-3495(61)86902-6

Hairer, E., Nørsett, S. P. y Wanner, G. (1993). Solving ordinary differential equations I: Nonstiff problems (2.ª ed. rev.). Springer. https://doi.org/10.1007/978-3-540-78862-1

Hairer, E. y Wanner, G. (1996). Solving ordinary differential equations II: Stiff and differential-algebraic problems (2.ª ed. rev.). Springer. https://doi.org/10.1007/978-3-642-05221-7

Kuehn, C. (2015). Multiple time scale dynamics. Springer. https://doi.org/10.1007/978-3-319-12316-5

Nagumo, J., Arimoto, S. y Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proceedings of the IRE, 50(10), 2061–2070. https://doi.org/10.1109/JRPROC.1962.288235

Strogatz, S. H. (2015). Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering (2.ª ed.). Westview Press.

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., … SciPy 1.0 Contributors. (2020). SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nature Methods, 17(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2