We present various examples where attracting limit cycles show up.
Cycling preys and predators: the Rosenzweig-McArthur model
We come back to preys and predators. We previously studied a basic model for the interaction of sharks and sardines, namely
$$ \begin{cases} \dot{x}= x(1-x) - xy\\ \dot{y}=\beta y ( x-\alpha) \end{cases} $$
where $x$ represents the density of sardines, $y$ the density of sharks, and $\alpha,\beta$ are positive parameters. We observed and proved that there cannot be limit cycles in this model, which means that the populations of sharks and sardines cannot oscillate periodically in a robust way.A natural question is: what biological factors create limit cycles in a prey-predator model? The inclusion of a more realistic of what ecologists call ‘functional response’ is one of such factor. The functional response is the rate at which each predator captures prey. In the above model, the functional response is a linearly increasing function of prey density, which means that there is no satiety for predators! A simple way to remedy this is to take a functional response of the form
$$ \frac{cx}{a+x} $$
where $a,c$ are positive parameters.
This modification leads to the following model, that we put in a nondimensional form:
$$ \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} $$
where $\alpha,\beta,\gamma$ are positive parameters. This model was introduced by Rosenzweig and McArthur.Invariance of the positive quadrant. We are only interested in the positive quadrant for the obvious reason that population densities are nonnegative numbers. In order that the model be well-defined, we thus need to check that starting in the positive quadrant, we have $(x(t),y(t))$ remaining in it for all $t>0$. We have
$$ \frac{\dot{x}}{x} = \frac{\text{d}}{\text{d} t}\ln x=1-\frac{x}{\gamma} -\frac{y}{1+x}. $$
Integrating both sides from time $0$ to time $t$, we get$$ 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). $$
From this we see that is $x_0=0$ then $x(t)=0$ for all $t>0$, and if $x_0>0$, then $x(t)>0$ for all $t>0$. We get the same conclusion for $y$. This proves more than theinvariance of the positive quadrant: it proves the invariance of its interior and that of its boundary.
Before playing with the model, let’s determine the nullclines and the fixed points.
Nullclines.
Setting $\dot{x}=0$, we find the $x$-nullclines:
$$ x=0\quad\text{and}\quad y=g(x)=(1+x)\left( 1-\frac{x}{\gamma}\right) $$
where we introduced the function $g(x)$ for later convenience. This function is a parabole whose peak is attained for$$ x=\frac{\gamma-1}{2}. $$
Setting $\dot{y}=0$, we find the $y$-nullclines:$$ y=0\quad\text{and}\quad x=\frac{\alpha}{1-\alpha}. $$
Fixed points.
For all parameter values, $(0,0)$ and $(\gamma,0)$ are fixed points. The former corresponds to the absence of both populations. The latter corresponds to the absence of predators, whereas preys have density $\gamma$. There is possibly a third fixed point:
$$ (\bar{x}, g(\bar{x})) $$
where$$ \bar{x}=\frac{\alpha}{1-\alpha}. $$
Since population densities are positive numbers, this fixed point is relevant only if $\bar{x}\leq \gamma$. This condition simply means that the nullclines $x=\alpha/(1-\alpha)$ and $y=g(x)$ intersect inside the positive quadrant.We further assume that $\alpha<1$, because otherwise all solutions starting from strictly positive initial densities tend to $(\gamma,0)$.
A supercritical Hopf bifurcation. Although there is no question that a limit cycle exists in this model, let us apply
Andronov-Hopf bifurcation theorem
, with $\alpha$ as a bifurcation parameter. (We could also apply
Poincaré-Bendixson theorem
.) We will content ourselves with verifying that the trace of the Jacobian matrix at the fixed point $(\bar{x},g(\bar{x}))$ changes sign as $\alpha$ passes through the value $(\gamma-1)/(\gamma+1)$, whereas the determinant is always positive. In passing, observe that in order to see the bifurcation occurring in the positive quadrant, $\gamma$ has to be strictly larger than $1$, which agrees with the numerical experimental. So, we assume that $\gamma>1$.
The Jacobian matrix at $(\bar{x},g(\bar{x}))$ is
$$ A= \begin{pmatrix} \alpha g’(\bar{x}) & -\alpha \\ \frac{\beta}{(1+\bar{x})^2} g(\bar{x}) & 0 \end{pmatrix}. $$
Since $g(x)>0$ whenever $x\in [0,\gamma[$, and $\bar{x}\in\thinspace ]0,\gamma[$, we have$$ \text{det}(A)=\frac{\alpha\beta}{(1+\bar{x})^2} g(\bar{x})>0. $$
We have$$ \text{tr}(A)=\alpha g’(\bar{x})=\frac{\gamma-1-2\bar{x}}{\gamma}. $$
Hence$$ \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} $$
The bifurcation parameter $\alpha^*$ is found by solving the equation$$ \frac{\alpha}{1-\alpha}=\frac{\gamma-1}{2} $$
which gives$$ \alpha^*=\frac{\gamma-1}{\gamma+1}. $$
A model of glycolysis
We look at a simple model of a fundamental biochemical process called glycolysis in which cells obtain energy by breaking down sugar. In dimensionless form, the equations are
$$ \begin{cases} \dot{x}= -x+ay+x^2 y\\ \dot{y}=b-ay-x^2 y \end{cases} $$
where $x$ and $y$ are, respectively, the concentrations of ADP (adenosine diphosphate) and F6P (fructose-6-phosphate), and $a,b>0$ are kinetic parameters. There is a single fixed point$$ \left(b\, ,\frac{b}{a+b^2}\right) $$
lying at the intersection of the nullclines$$ y=\frac{x}{a+x^2}\quad\text{and}\quad y=\frac{b}{a+x^2}. $$
We can observe that there is a limit cycle showing up for certain values of $(a,b)$. (We explain below the shape of the region of the $ab$-plane in which we have a limit cycle.)
Poincaré-Bendixson theorem in action. We are going to apply this theorem to rigorously prove that there exists a limit cycle in this system. To achieve this aim, we construct a trapping region. We claim that the ‘box’ bounded by the red line shown in the figure is such a region.
$$
\dot{x}-(-\dot{y})= -x+ay+x^2y+(b-ay-x^2y)=b-x.
$$ $$
-\dot{y}>\dot{x}\quad\text{if}\quad x>b.
$$ $$
\frac{\text{d}y}{\text{d}x}=\frac{\text{d}y/\text{d}t}{\text{d}x/\text{d}t}
=\frac{\dot{y}}{\dot{x}}
$$
We have to show that all the vectors on the boundary of the box point into it. On the horizontal and vertical sides, this is left to the reader. The tricky part is to deal with the diagonal line of slope $-1$ extending from the point $(b,b/a)$ to the nullcline $y=x/(a+x^2)$. Let us consider the quantity $\dot{x}-(-\dot{y})$. We find
Can we now apply Poincaré-Bendixson theorem and conclude there exists a limit cycle in the box? Not quite yet, because there is a fixed point inside! But suppose that the fixed point is a repeller (nodal or spiral source), then we are in good shape by considering a ‘punctured’ version of the box, as shown in the figure.
The repeller drives all neighboring solutions into the hatched region, and the region is free of fixed points: Poincaré-Bendixson theorem applies!
Now we look under which conditions the fixed point is a repeller. The Jacobian matrix at the fixed point is
$$ A= \begin{pmatrix} -1+\frac{2b^2}{a+b^2} & a+b^2 \\ -\frac{2b^2}{a+b^2} & -(a+b^2) \end{pmatrix}. $$
We have $\text{det}(A)=a+b^2>0$, hence the fixed point is never a saddle. Next we find$$ \text{tr}(A)=-\frac{b^4+(2a-1)b^2+a+a^2}{a+b^2}. $$
The fixed point is repulsive for $\text{tr}(A)>0$ and attractive for $\text{tr}(A)<0$. The dividing line $\text{tr}(A)=0$ in the $ab$-plane is formed by the two curves$$ b=\sqrt{\frac{1}{2}\left(1-2a\pm\sqrt{{\scriptsize 1-8a}}\right)}\, . $$
A model for a chemical reaction that oscillates
The Belousov-Zhabotinsky reaction was a major turning point in the history of chemistry, since it showed that certain chemical reactions could oscillate for long periods of time before settling to equilibrium. It was discovered by Belousov in the 1950s, but his work was only disseminated at the end of the 1960s. Before this, chemists believed that all chemical reactions tended monotonically to equilibrium.
There is by now a large class of chemical reactions displaying oscillations. One particularly simple one is given by a chlorine dioxide-iodine-malonic acid (ClO$_2$-I$_2$-MA) reaction. The exact differential equations modeling this reaction are extremely complicated. Fortunately, there is a way of approximating the reaction because the concentration of some of the reactants varies much more slowly than some others. The idea is to treat the concentration of slow reactants as constant. By doing so, the model will not account for the eventual approach to equilibrium, but it can account for oscillations. After suitable nondimensionalization, one end up with the following two-dimensional model:
$$ \begin{cases} \dot{x}= a-x-\frac{4xy}{1+x^2}\\ \dot{y}=bx\left( 1-\frac{y}{1+x^2}\right) \end{cases} $$
where $x$ and $y$ are the dimensionless concentrations of $\text{I}^-$ and $\text{ClO}_2^-$, respectively, and $a,b$ are positive parameters.We have $\dot{x}=0$ on the curve
$$ y=\frac{(a-x)(1+x^2)}{4x} $$
and $\dot{y}=0$ on the $y$-axis and on the parabola
$$ y=1+x^2. $$
These two curves intersect at a single point, so we have only one fixedpoint:
$$ \left( \frac{a}{5}\, , \,1+\left(\frac{a}{5}\right)^2\right). $$
You can observe that for certain parameters there is an attracting limit cycle. You can also observe a supercritical Hopf bifurcation.
The reader can look for a trapping region, as in the previous example, to apply Poincaré-Bendixson theorem. It can be shown that the curve
$$ b=\frac{3a}{4}-\frac{25}{a} $$
divides the $ab$-plane into two regions: for $(a,b)$ below the curve, we have an attracting limit cycle surrounding the fixed point which is a repeller (spiral source), and above the curve, there is no limit cycle and the fixed point is attractive (it is a spiral sink).A prey-predator model with Allee effect for preys
We consider a prey-predator model which is obtained again from the model
$$ \begin{cases} \dot{x}= x(1-x) - xy\\ \dot{y}=\beta y ( x-\alpha). \end{cases} $$
This time, we don’t modify the interaction term $xy$ (as we did before to get Rosenzweig-McArthur model), but the logistic term $x(1-x)$. We assume that there is a threshold to growth for the prey population (in the absence of predators). The simplest way to do this is to set$$ \dot{x} = x\left(\frac{x}{\gamma}-1\right)(1-x) $$
where $0<\gamma<1$. This model is easy to analyse graphically.
We thus arrive at the following model:
$$ \begin{cases} \dot{x} = x\left(\frac{x}{\gamma}-1\right)(1-x)-xy\\ \dot{y} = \beta y (x-\alpha) \end{cases} $$
where $\alpha,\beta,\gamma$ are positive parameters with $\gamma<1$. Proceeding as for Rosenzweig-McArthur model, one can check that the positive quadrant is invariant (in fact, its interior and its boundary are both invariant).Fixed points. There are three fixed points which exist for all parameter values:
$$ (0,0),\; (\gamma,0),\; (1,0). $$
A fourth one can come into play if the nullclines $x=\alpha$ and $y=g(x)=(\frac{x}{\gamma}-1)(1-x)$ intersect inside the positive quadrant, namely$$ (\alpha,g(\alpha)), $$
which is the case if and only if $\gamma\leq \alpha\leq 1$. It corresponds to coexistence of both populations. As $\alpha$ varies, we are going to observe a (supercritical) Hopf bifurcation.
The observed behavior of fixed points can be confirmed by linearization. One can also check the basic mechanism of the Hopf bifurcation.
$$
A=
\begin{pmatrix}
xg’(x)+g(x)-y & -x \\
\beta y & \beta(x-\alpha)
\end{pmatrix}.
$$
The Jacobian matrix at a point $(x,y)$ is
The fixed point $(\gamma,0)$ is a source if $\alpha<\gamma$, and a saddle if $\alpha>\gamma$.
The fixed point $(1,0)$ is a saddle if $\alpha<1$, and a sink if $\alpha>1$.
The fixed point $(\alpha,g(\alpha))$, which exists only if $\alpha\in [\gamma,1]$, is sink if $\alpha>\frac{\gamma+1}{2}$, and a source if $\alpha>\frac{\gamma+1}{2}$. One can check that $\text{det}(A)=\alpha\beta g(\alpha)>0$ for all $\alpha\in ]\gamma,1[$, whereas $\text{tr}(A)$ changes it sign from positive to negative when $\alpha$ increases
and pass through the value $(\gamma+1)/2$.
Neurodynamics: the FitzHugh-Nagumo model
The FitzHugh-Nagumo model is a model of the Hodgkin-Huxley model. Although this model is not as biologically accurate as the original model, it nevertheless does capture the essential behavior of nerve impulses, including the phenomenon called excitability, to be described below. The model is defined by the equations
$$ \begin{cases} \dot{x}= x-\frac{x^3}{3} -y +I \\ \dot{y}=\frac{1}{\tau}(x+\alpha-\beta y). \end{cases} $$
In these equations, $x$ is similar to the voltage and represents the excitability of the system; the variable $y$ represents a combination of other forces that tend to return the system to rest (feedback). There are four positive parameters: $\alpha,\beta,\tau,I$. The parameter $I$ is a stimulus parameter that leads to excitation of the system; $I$ is like an applied stimulus current. The parameter $\tau$ is a time scale that allows to make $x$ evolving faster than $y$; roughly speaking $\dot{x}/\dot{y}\approx \tau$.The curve on which $\dot{x}=0$ is the cubic function
$$ y=x-\frac{x^3}{3}+I $$
and $\dot{y}=0$ on the line curve$$ y=\frac{x+\alpha}{\beta}. $$
There is a single fixed point $(\bar{x},\bar{y})$ which is the solution of a cubic equation.Typical values for the parameters $\alpha,\beta,\tau$ are $0.7, 0.8,13$, respectively. We will only tune the parameter $I$. We shall observe what is called the excitation block phenomenon.
An example with infinitely many limit cycles
We consider the system
$$ \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} $$
The digital experiment clearly suggests that there are infinitely many circular limit cycles.
Passing to polar coordinates, we can analyse the system quite easily:
$$ \begin{cases} \dot{r}=r^5 \sin\left(\frac{1}{r^2}\right)=f(r)\\ \dot{\theta}=-1. \end{cases} $$
We get decoupled equations. The second one is simply a counterclockwise rotation at constant speed. The equation for $r(t)$ has fixed points at$$ \frac{1}{r^2}=n\pi,\; n\;\text{any positive integer}. $$
Let us apply the derivative test to the one-dimensional system $\dot{r}=f(r)$. We have
$$ f’(r) = r^2\left( 5 r^2 \sin\left(\frac{1}{r^2}\right)-2 \cos\left(\frac{1}{r^2}\right)\right). $$
At the points $\frac{1}{r^2}=n\pi$, we have $\sin\left(\frac{1}{r^2}\right)=0$ and
$$ 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} $$
Therefore, $\frac{1}{\sqrt{2\pi}}, \frac{1}{\sqrt{4\pi}}$,... are stable fixed points, and $\frac{1}{\sqrt{\pi}}, \frac{1}{\sqrt{3\pi}}$,... are unstable fixed points. Back to Cartesian coordinates, this means that we have infinitely many limit cycles which are circles of radius $\frac{1}{\sqrt{n\pi}}$, where those with $n$ even are attracting limit cycles, those with $n$ odd are repelling limit cycles.Not all limit cycles are convex curves
All examples of limit cycles in the plane we’ve seen so far are convex curves. A curve is convex if it lies completely on one side of each and every one of its tangent lines. The next example shows that this need not be the case, although such examples are not so easy to find out!