Presentamos varios ejemplos en los que aparecen ciclos límite atrayentes.
Presas y depredadores cíclicos: el modelo Rosenzweig-McArthur
Volvemos a las presas y los depredadores. Anteriormente estudiamos un modelo básico para la interacción entre tiburones y sardinas, a saber
$$ \begin{cases} \dot{x}= x(1-x) - xy\\ \dot{y}=\beta y ( x-\alpha) \end{cases} $$
donde $x$ representa la densidad de sardinas, $y$ la densidad de tiburones, y $\alpha,\beta$ son parámetros positivos. Observamos y demostramos que no puede haber ciclos límite en este modelo, lo que significa que las poblaciones de tiburones y sardinas no pueden oscilar periódicamente de forma robusta.
Una pregunta natural es: ¿qué factores biológicos crean ciclos límite en un modelo presa-predador? Uno de esos factores es la inclusión de lo que los ecólogos denominan "respuesta funcional". La respuesta funcional es la velocidad a la que cada depredador captura a su presa. En el modelo anterior, la respuesta funcional es una función linealmente creciente de la densidad de presas, lo que significa que los depredadores no se sacian. Una forma sencilla de remediarlo es tomar una respuesta funcional de la forma
$$ \frac{cx}{a+x} $$
donde $a,c$ son parámetros positivos.
Esta modificación conduce al siguiente modelo, que ponemos en forma adimensional:
$$ \begin{cases} \dot{x}= x\left( 1-\frac{x}{\gamma}\right) - \frac{xy}{1+x}\\ \dot{y}=\beta y \left( \frac{x}{1+x}-\alpha\right) \end{cases} $$
donde $\alpha,\beta,\gamma$ son parámetros positivos. Este modelo fue introducido por Rosenzweig y McArthur.
Invarianza del cuadrante positivo. Sólo nos interesa el cuadrante positivo por la razón obvia de que las densidades de población son números no negativos. Para que el modelo esté bien definido, necesitamos comprobar que partiendo del cuadrante positivo, tenemos $(x(t),y(t))$ que permanecen en él para todo $t>0$. Tenemos
$$ \frac{\dot{x}}{x} = \frac{\text{d}}{\text{d} t}\ln x=1-\frac{x}{\gamma} -\frac{y}{1+x}. $$
Integrando ambos lados desde el tiempo $0$ hasta el tiempo $t$, obtenemos
$$ x(t)=x_0 \, \exp\left( \int_0^t \left(1-\frac{x(s)}{\gamma} -\frac{y(s)}{1+x(s)}\right) \mathrm{d}s\right). $$
De aquí se deduce que si $x_0=0$ entonces $x(t)=0$ para todo $t>0$, y si $x_0>0$, entonces $x(t)>0$ para todo $t>0$. Obtenemos la misma conclusión para $y$. Esto demuestra algo más que la invariancia del cuadrante positivo: demuestra la invariancia de su interior y la de su frontera.
Antes de jugar con el modelo, determinemos las líneas nulas y los puntos fijos.
Nullclines.
Fijando $\dot{x}=0$, encontramos las $x$-nullclines:
$$ x=0\quad\text{and}\quad y=g(x)=(1+x)\left( 1-\frac{x}{\gamma}\right) $$
donde introducimos la función $g(x)$ por comodidad posterior. Esta función es una parábola cuyo pico se alcanza para
$$ \frac{\gamma-1}{2}. $$
Fijando $\dot{y}=0$, encontramos las $y$-nulas:
$$ y=0\quad\text{and}\quad x=\frac{\alpha}{1-\alpha}. $$
Puntos fijos.
Para todos los valores de los parámetros, $(0,0)$ y $(\gamma,0)$ son puntos fijos. El primero corresponde a la ausencia de ambas poblaciones. El segundo corresponde a la ausencia de depredadores, mientras que las presas tienen una densidad $\gamma$. Es posible que exista un tercer punto fijo:
$$ (\bar{x}, g(\bar{x})) $$
donde
$$ \bar{x}=\frac{\alpha}{1-\alpha}. $$
Dado que las densidades de población son números positivos, este punto fijo sólo es relevante si $\bar{x}\leq \gamma$. Esta condición significa simplemente que las líneas nulas $x=alpha/(1-\alpha)$ e $y=g(x)$ se cruzan dentro del cuadrante positivo.
Suponemos además que $\alpha<1$, porque de lo contrario todas las soluciones que parten de densidades iniciales estrictamente positivas tienden a $(\gamma,0)$.
Una bifurcación de Hopf supercrítica. Aunque no hay duda de que en este modelo existe un ciclo límite, apliquemos Andronov-Hopf bifurcation theorem
, con $\alpha$ como parámetro de bifurcación. (También podríamos aplicar
Teorema de Poincaré-Bendixson
.)
Nos contentaremos con comprobar que la traza de la matriz jacobiana en el punto fijo $(\bar{x},g(\bar{x}))$ cambia de signo a medida que $\alpha$ pasa por el valor $(\gamma-1)/(\gamma+1)$, mientras que el determinante es siempre positivo. De paso, obsérvese que para que la bifurcación se produzca en el cuadrante positivo, $\gamma$ tiene que ser estrictamente mayor que $1$, lo que concuerda con el experimento numérico. Por tanto, supondremos que $\gamma>1$.
La matriz jacobiana en $(\bar{x},g(\bar{x}))$ es
$$ A= \begin{pmatrix} \alpha g’(\bar{x}) & -\alpha \\ \frac{\beta}{(1+\bar{x})^2} g(\bar{x}) & 0 \end{pmatrix}. $$
Dado que $g(x)>0$ siempre que $x\in [0,\gamma[$, y $\bar{x}\in\thinspace ]0,\gamma[$, tenemos
$$ \text{det}(A)=\frac{\alpha\beta}{(1+\bar{x})^2} g(\bar{x})>0. $$
Tenemos
$$ \text{tr}(A)=\alpha g’(\bar{x})=\frac{\gamma-1-2\bar{x}}{\gamma}. $$
Por lo tanto
$$ \text{tr}(A) \begin{cases} <0 & \text{if}\quad \bar{x}>\frac{\gamma-1}{2}\\ =0 & \text{if}\quad \bar{x}=\frac{\gamma-1}{2}\\ >0 & \text{if}\quad \bar{x}<\frac{\gamma-1}{2}. \end{cases} $$
El parámetro de bifurcación $\alpha^*$ se encuentra resolviendo la ecuación
$$ \frac{\alpha}{1-\alpha}=\frac{\gamma-1}{2} $$
que da
$$ \alpha^*=\frac{\gamma-1}{\gamma+1}. $$
Un modelo de glucólisis
Examinamos un modelo sencillo de un proceso bioquímico fundamental denominado glucólisis, en el que las células obtienen energía descomponiendo azúcar. En forma adimensional, las ecuaciones son
$$ \begin{cases} \dot{x}= -x+ay+x^2 y\\ \dot{y}=b-ay-x^2 y \end{cases} $$
donde $x$ e $y$ son, respectivamente, las concentraciones de ADP (adenosín difosfato) y F6P (fructosa-6-fosfato), y $a,b>0$ son parámetros cinéticos. Existe un único punto fijo
$$ \left(b\, ,\frac{b}{a+b^2}\right) $$
situada en la intersección de las líneas nulas
$$ y=\frac{x}{a+x^2}\quad\text{and}\quad y=\frac{b}{a+x^2}. $$
Podemos observar que aparece un ciclo límite para ciertos valores de $(a,b)$. (Explicamos a continuación la forma de la región del plano $ab$ en la que tenemos un ciclo límite).
Teorema de Poincaré-Bendixson en acción. Vamos a aplicar este teorema para demostrar rigurosamente que existe un ciclo límite en este sistema. Para ello, construimos una región trampa. Afirmamos que la ‘caja’ delimitada por la línea roja que se muestra en la figura es una región de este tipo.
¿Podemos aplicar ahora el teorema de Poincaré-Bendixson y concluir que existe un ciclo límite en la caja? Todavía no, ¡porque hay un punto fijo dentro! Pero supongamos que el punto fijo es un repulsor (fuente nodal o espiral), entonces estamos en buena forma considerando una versión ‘perforada’ de la caja, como se muestra en la figura.
El repulsor conduce todas las soluciones vecinas a la región sombreada, y la región está libre de puntos fijos: Se aplica el teorema de Poincaré-Bendixson.
Ahora veremos en qué condiciones el punto fijo es un repulsor. La matriz jacobiana en el punto fijo es
$$ A= \begin{pmatrix} -1+\frac{2b^2}{a+b^2} & a+b^2 \\ -\frac{2b^2}{a+b^2} & -(a+b^2) \end{pmatrix}. $$
Tenemos $\text{det}(A)=a+b^2>0$, por lo tanto el punto fijo nunca es una silla de montar. A continuación encontramos
$$ \text{tr}(A)=-\frac{b^4+(2a-1)b^2+a+a^2}{a+b^2} $$
El punto fijo es repulsivo para $\text{tr}(A)>0$ y atractivo para $\text{tr}(A)<0$. La línea divisoria $\text{tr}(A)=0$ en el plano $ab$ está formada por las dos curvas
$$ b=\sqrt{\frac{1}{2}\left(1-2a\pm\sqrt{{\scriptsize 1-8a}}\right)}\, . $$
Modelo de reacción química oscilante
La reacción Belousov-Zhabotinsky marcó un hito en la historia de la química, ya que demostró que determinadas reacciones químicas podían oscilar durante largos periodos de tiempo antes de alcanzar el equilibrio. Fue descubierta por Belousov en la década de 1950, pero sus trabajos no se difundieron hasta finales de los años sesenta. Antes, los químicos creían que todas las reacciones químicas tendían monotónicamente al equilibrio.
En la actualidad, existe una gran clase de reacciones químicas que muestran oscilaciones. Una de las más sencillas es la reacción dióxido de cloro-ácido yodo-malónico (ClO$_2$-I$_2$-MA). Las ecuaciones diferenciales exactas que modelan esta reacción son extremadamente complicadas. Afortunadamente, existe una forma de aproximar la reacción porque la concentración de algunos de los reactantes varía mucho más lentamente que la de otros. La idea es tratar la concentración de los reactantes lentos como constante. De este modo, el modelo no tendrá en cuenta la aproximación final al equilibrio, pero sí las oscilaciones. Tras una adecuada adimensionalización, se obtiene el siguiente modelo bidimensional:
$$ \begin{cases} \dot{x}= a-x-\frac{4xy}{1+x^2}\\ \dot{y}=bx\left( 1-\frac{y}{1+x^2}\right) \end{cases} $$
donde $x$ e $y$ son las concentraciones adimensionales de $\text{I}^-$ y $\text{ClO}_2^-$, respectivamente, y $a,b$ son parámetros positivos.
Tenemos $\dot{x}=0$ en la curva
$$ y=\frac{(a-x)(1+x^2)}{4x} $$
y $\dot{y}=0$ en el eje $y$ y en la parábola
$$ y=1+x^2. $$
Estas dos curvas se cruzan en un único punto, por lo que sólo tenemos un punto fijo:
$$ \left( \frac{a}{5}\, , \,1+\left(\frac{a}{5}\right)^2\right). $$
Se puede observar que para ciertos parámetros existe un ciclo límite de atracción. También se puede observar una bifurcación de Hopf supercrítica.
El lector puede buscar una región de atrapamiento, como en el ejemplo anterior, para aplicar el teorema de Poincaré-Bendixson. Se puede demostrar que la curva
$$ b=\frac{3a}{4}-\frac{25}{a} $$
divide el plano $ab$ en dos regiones: para $(a,b)$ por debajo de la curva, tenemos un ciclo límite atrayente que rodea al punto fijo que es un repulsor (fuente espiral), y por encima de la curva, no hay ciclo límite y el punto fijo es atractivo (es un sumidero espiral).
Un modelo presa-predador con efecto Allee para las presas
Consideramos un modelo presa-predador que se obtiene de nuevo a partir del modelo
$$ \begin{cases} \dot{x}= x(1-x) - xy\\ \dot{y}=\beta y ( x-\alpha). \end{cases} $$
Esta vez, no modificamos el término de interacción $xy$ (como hicimos antes para obtener Rosenzweig-McArthur model), sino el término logístico $x(1-x)$. Suponemos que existe un umbral de crecimiento para la población de presas (en ausencia de depredadores). La forma más sencilla de hacerlo es establecer
$$ \dot{x} = x\left(\frac{x}{\gamma}-1\right)(1-x) $$
donde $0<\gamma<1$. Este modelo es fácil de analizar gráficamente.
Llegamos así al siguiente modelo:
$$ \begin{cases} \dot{x} = x\left(\frac{x}{\gamma}-1\right)(1-x)-xy\\ \dot{y} = \beta y (x-\alpha) \end{cases} $$
donde $\alpha,\beta,\gamma$ son parámetros positivos con $\gamma<1$. Procediendo como para el modelo de Rosenzweig-McArthur, se puede comprobar que el cuadrante positivo es invariante (de hecho, tanto su interior como su frontera son invariantes).
Puntos fijos. Hay tres puntos fijos que existen para todos los valores de los parámetros:
$$ (0,0),\; (\gamma,0),\; (1,0). $$
Una cuarta puede entrar en juego si los nuloclinos $x=\alpha$ e $y=g(x)=(\frac{x}{\gamma}-1)(1-x)$ se cruzan dentro del cuadrante positivo, a saber
$$ (\alpha,g(\alpha)), $$
que es el caso si y sólo si $\gamma\leq \alpha\leq 1$. Corresponde a la coexistencia de ambas poblaciones. A medida que $\alpha$ varía, vamos a observar una bifurcación de Hopf (supercrítica).
El comportamiento observado de los puntos fijos puede confirmarse mediante linealización. También se puede comprobar el mecanismo básico de la bifurcación de Hopf.
Neurodinámica: el modelo FitzHugh-Nagumo
El modelo FitzHugh-Nagumo es un modelo del modelo Hodgkin-Huxley. Aunque este modelo no es tan biológicamente exacto como el modelo original, sin embargo capta el comportamiento esencial de los impulsos nerviosos, incluyendo el fenómeno llamado excitabilidad, que se describirá más adelante. El modelo está definido por las ecuaciones
$$ \begin{cases} \dot{x}= x-\frac{1}{3} x^3 -y +I \\ \dot{y}=\frac{1}{\tau}(x+\alpha-\beta y). \end{cases} $$
En estas ecuaciones, $x$ es similar al voltaje y representa la excitabilidad del sistema; la variable $y$ representa una combinación de otras fuerzas que tienden a devolver el sistema al reposo (retroalimentación). Hay cuatro parámetros positivos: $\alpha,\beta,\tau,I$. El parámetro $I$ es un parámetro de estímulo que conduce a la excitación del sistema; $I$ es como una corriente de estímulo aplicada. El parámetro $\tau$ es una escala de tiempo que permite hacer que $x$ evolucione más rápido que $y$; a grandes rasgos $\dot{x}/\dot{y}approx \tau$.
La curva en la que $\dot{x}=0$ es la función cúbica
$$ y=x-\frac{1}{3}x^3+I $$
y $\dot{y}=0$ en la recta curva
$$ y=\frac{x+\alpha}{\beta}. $$
Existe un único punto fijo $(\bar{x},\bar{y})$ que es la solución de una ecuación cúbica.
Los valores típicos para los parámetros $\alpha,\beta,\tau$ son $0,7, 0,8,13$, respectivamente. Nosotros sólo afinaremos el parámetro $I$. Observaremos lo que se denomina el fenómeno de bloque de excitación.
Un ejemplo con infinitos ciclos límite
Consideramos el sistema
$$ \begin{cases} \dot{x}= y + (x^2+y^2)^2 x\sin\left(\frac{1}{x^2+y^2}\right)\\ \dot{y}= -x + (x^2+y^2)^2 y\sin\left(\frac{1}{x^2+y^2}\right). \end{cases} $$
El experimento digital sugiere claramente que existen infinitos ciclos límite circulares.
Pasando a coordenadas polares, podemos analizar el sistema con bastante facilidad:
$$ \begin{cases} \dot{r}=r^5 \sin\left(\frac{1}{r^2}\right)=f(r)\\ \dot{\theta}=-1. \end{cases} $$
Obtenemos ecuaciones desacopladas. La segunda es simplemente una rotación en sentido contrario a las agujas del reloj a velocidad constante. La ecuación para $r(t)$ tiene puntos fijos en
$$ \frac{1}{r^2}=n\pi,\quad n\thinspace\text{any positive integer}. $$
Apliquemos la prueba de la derivada al sistema unidimensional $\dot{r}=f(r)$. Tenemos
$$ f’(r) = r^2\left[ 5 r^2 \sin\left(\frac{1}{r^2}\right)-2 \cos\left(\frac{1}{r^2}\right)\right]. $$
En los puntos $\frac{1}{r^2}=n\pi$, tenemos $\sin\left(\frac{1}{r^2}\right)=0$ y
$$ f’(r)=-2r\cos(n\pi)\thinspace\thinspace\text{is}\thinspace \begin{cases} \text{positive} & \text{if}\thinspace\thinspace n\thinspace\text{odd}\\ \text{negative} & \text{if}\thinspace\thinspace\thinspace n\thinspace \text{even}. \end{cases} $$
Por lo tanto, $\frac{1}{sqrt{2\pi}}, \frac{1}{sqrt{4\pi}}$,... son puntos fijos estables, y $\frac{1}{sqrt{\pi}}, \frac{1}{sqrt{3\pi}$,... son puntos fijos inestables. Volviendo a las coordenadas cartesianas, esto significa que tenemos infinitos ciclos límite que son círculos de radio $\frac{1}{\sqrt{n\pi}}$, donde los que tienen $n$ pares son ciclos límite que atraen, los que tienen $n$ impares son ciclos límite que repelen.
No todos los ciclos límite son curvas convexas
Todos los ejemplos de ciclos límite en el plano que hemos visto hasta ahora son curvas convexas. Una curva es convexa si se encuentra completamente a un lado de todas y cada una de sus rectas tangentes. El siguiente ejemplo muestra que no tiene por qué ser así, ¡aunque no es tan fácil encontrar ejemplos así!