Analysis of frequency and amplitude at Hopf bifurcation

Question

I am analyzing the following system, where $I_{in}$ is a scalar parameter: $$ \begin{aligned} &\dot{V} = 10 \left( V - \frac{V^3}{3} - R + I_{in} \right) \\ &\dot{R} = 0.8 \left( -R +1.25V + 1.5 \right) \end{aligned} $$

It is a simplified version of the Fitzhugh-Nagumo equations for neuronal excitability (reference to book below).

There is a single equilibrium, that varies with $I_{in}$, so we need to calculate the Jacobian at those equilibrium values and perform a stability analysis. Such an analysis reveals that as $I_{in}$ increases from zero to around 1.5, the system undergoes a supercritical Hopf bifurcation [edit: it undergoes what I intuitively thought was a supercritical bifurcation]:

We go from a stable center to an unstable center at a critical value (zero real eigenvalue) at $I_{crit}=0.966064$. Note that I calculated the limit cycle boundaries in that diagram by just getting the minimum and maximum values of V for each loop through the limit cycle (examples of such loops are shown below in Figures 3 and 4).

(Edit: I added my derivation of $I_{crit}$ below).

You can see the nature of the transition in the trace-determinant diagram in the following Figure 2:

As $I_{in}$ increases, the equilibrium point turns from a sink to a (stable) spiral, and then we hit the critical point at $I_{crit}$, after which we have a spiral source surrounded by a (stable) limit cycle, and eventually a stable source also surrounded by a limit cycle.

So far, so good, I think. This all seems pretty straightforward.

So what is the problem? At this point I am confused about a couple of things. In my book it says the following two facts (corollaries of the Hopf Bifurcation theorem) should be true near $I_{crit}$:

1. The amplitude of oscillations will be very small.
2. The frequency of oscillation should be close to $\omega/2\pi$ Hz, where $\omega$ is the imaginary part of the eigenvalue.

It seems that neither of these facts is true here.

First, the oscillation amplitude starts out very large, as you can see in the bifurcation diagram in Figure 1. There is none of that textbook ramping up of the amplitude.

Indeed, even when $I_{in}$ is less than $I_{crit}$, there already a large, stable, limit-cycle-like orbit in this phase space! The following figure shows some full orbits in the phase space (left), and a couple of V trajectories on the right. This is for $I_{in}=I_{crit}-0.00874$:

That is, we have lots of large-amplitude stable orbits cycling around some stable centers (such damped oscillations occur only for orbits that are close to the equilibrium point). So not only does the limit cycle start out with a large amplitude past $I_{crit}$, it seems there is already a kind of harbinger of a limit cycle with a large amplitude even before the bifurcation.

That said, the above two facts do seem to apply to the damped spirals in Figure 2: the amplitude of the spiral is very small (tending toward zero), and its frequency is basically exactly $\omega/2\pi$ -- it is basically double the frequency of the large-amplitude pseudo-limit cycle that encloses it. Could that be what my text is referring to?

This brings me explicitly to the second fact above: at $I_{crit}$ the eigenvalues are $\pm 3.05i$. Hence, the frequency of oscillation should be about 0.5 Hz, a period of 2 s. But instead I see a period of 4 seconds (0.25 Hz), as the following diagram of V versus time for $I=I_{crit}+0.000001$ shows:

I calculate the period based on the distance between the red X's. I'll mention again, though, if we were to do the same analysis of the oscillations of the damped spirals (as in Figure 3) the frequencies of those damped oscillations would basically be right -- it is the full limit cycles that seem off (though their order of magnitude is right).

Overall, this system is supposed to be approachable because of its simplicity but I've already spent about a week hitting my head against it, still not sure of some of the most basic facts about Hopf Bifurcations.

Derivation of the critical value

Note the Jacobian of the system is: $$ J = \begin{bmatrix} \frac{\partial F_1}{\partial V} & \frac{\partial F_1}{\partial R}\\ \frac{\partial F_2}{\partial V} & \frac{\partial F_2}{\partial R} \end{bmatrix} = \begin{bmatrix} 10(1-V^2) & -10 \\ 1 & -0.8 \end{bmatrix} $$ Our task is essentially to determine the system's (V,R) equilibria for different values of $I_{in}$. Then, we can plug these equilibrium values into the Jacobian for our stability analysis, and find the coefficient matrix where the real part of the eigenvalues go to zero.

How do we find this? First, I found the equilibrium value of V that would yield purely imaginary eigenvalues, and I did this using the trace. That is, the sum of the eigenvalues is the same as the sum of the values in the coefficient matrix (the trace). From the equation for the Jacobian above, we know the trace is zero when:

$$ 9.2 - 10V^2 = 0 \implies V = \pm \sqrt{0.92} = \pm 0.959 $$

Focusing on the negative root for now, this implies that our critical value of $I_{in}$ will be the one that generates $V_{eqm}=-0.959$.

How do we find this value of $I_{in}$? I did it by substitution, using the nullcline equations of our system. Namely, the equations for the nullclines of our system are given by:

$$ \begin{aligned} &R = V - \frac{V^3}{3} + I_{input}\\ &R = 1.25V + 1.5 \end{aligned} $$

So, if given a value $V_{eqm}$ we can plug the second nullcline equation $R(V)$ into the first, and solve for $I_{in}$ in terms of $V$. Namely, given a value of $V_{eqm}$, the $I_{in}$ that produces that will be:

$$ I_{in}=\frac{V^3}{3} + 0.25V + 1.5 $$

So, looping back to our question, if we plug in $V_{eqm}=-0.959$ into this equation, that yields $I_{crit}=0.966$. Note also that plugging this $I_{crit}$ into the original system of equations and numerically solving for the equilibrium using Python's fsolve() yields the equilibrium point (V, R) = (-0.959, 0.301), which gives secondary confirmation of our result.

The Jacobian at this equilibrium point is:

$$ J = \begin{bmatrix} 0.8 & -10 \\ 1 & -0.8 \end{bmatrix} $$

This coefficient matrix has purely imaginary eigenvalues $\pm3.06i$, as expected. So it seems we have a critical value where the real part of the eigenvalues reaches zero, as originally claimed. QED, maybe?

To address a question from a comment: when $I=0.866$ the equilibrum point is $(V, R) = (-1.04, 0.20)$, and the eigenvalues of the Jacobian are $-0.8\pm3.16i$. This, coupled with the secondary confirmation of the calculations from the trace-discriminant curve (Figure 2 above), make me think there isn't a mistake in the calculation of $I_{crit}$ above. That said, I have definitely made worse mistakes in my lifetime, and thought I was right, so we definitely shouldn't exclude this possibility.

Different question about same equations Hopf bifurcation and limit cycles

Reference
Wilson (1999) Spikes, decisions, and actions: dynamical foundations of neuroscience. OUP.

Answer 1

Not an answer. I leave this MATHEMATICA script as a calculation process to the critical data.

f[i0_, v_, r_] := {10 (v - v^3/3 - r + i0), 0.8 (-r + 1.25 v + 1.5)};
sols = Quiet@Solve[f[i0, x, y] == 0, {x, y}];
J0 = Grad[f[i0, x, y], {x, y}] /. sols[[1]];
eig = Eigenvalues[J0];
Plot[Re[eig], {i0, 0, 3}]

With the plot feedback we can search for the solutions for $I_{crit}$

soli1 = NMinimize[Re[eig].Re[eig], i0]
soli2 = NMinimize[{Re[eig].Re[eig], i0 >= 1.5}, i0]
icrit1 = i0 /. soli1[[2]]
icrit2 = i0 /. soli2[[2]]

For $0.966064 \le I \le 2.03394$ we have a stable limit cycle

tmax = 10;
solxy = NDSolve[Join[Thread[D[{x[t], y[t]}, t] == f[icrit1, x[t], y[t]]], {x[0] == 0, y[0] == 1.2}], {x, y}, {t, 0, tmax}][[1]];
gr1 = ParametricPlot[Evaluate[{x[t], y[t]} /. solxy], {t, 0, tmax}, PlotStyle -> Red];
solxy = NDSolve[Join[Thread[D[{x[t], y[t]}, t] == f[icrit2, x[t], y[t]]], {x[0] == 0, y[0] == 1.2}], {x, y}, {t, 0, tmax}][[1]];
gr1b = ParametricPlot[Evaluate[{x[t], y[t]} /. solxy], {t, 0, tmax}, PlotStyle -> Blue];
Show[gr1, gr1b, PlotRange -> All]

and for $I$ outside this interval

eps = 0.01; 
solxy = NDSolve[Join[Thread[D[{x[t], y[t]}, t] == f[icrit1 - eps, x[t], y[t]]], {x[0] == 0, y[0] == 1.2}], {x, y}, {t, 0, tmax}][[1]];
gr1 = ParametricPlot[Evaluate[{x[t], y[t]} /. solxy], {t, 0, tmax}, PlotStyle -> Red];
solxy = NDSolve[Join[Thread[D[{x[t], y[t]}, t] == f[icrit2 + eps, x[t], y[t]]], {x[0] == 0, y[0] == 1.2}], {x, y}, {t, 0, tmax}][[1]];
gr1b = ParametricPlot[Evaluate[{x[t], y[t]} /. solxy], {t, 0, tmax}, PlotStyle -> Blue];
Show[gr1, gr1b, PlotRange -> All]

Answer 2

Using the approach in this Scholarpedia article, second row, I am obtaining positive Lyapunov exponents for both Hopf bifurcations. Specifically, setting $$ x_1 = V, \; x_2 = 10(R - 1.5), \; F(x_1) = 10(x_1 - x_1^3/3) + \lambda, \; a = 0.08, b = 1.25 $$ the problem takes the form in row 2 of the set of examples, using a suitable $\lambda$. We now must compute the sign of $$ F'''(x_1^\ast) + \frac{F''(x_1^\ast)^2}{b-a} = -20 + \frac{(-20x_1^\ast)^2}{1.17} $$ for the $x_1$ values at the Hopf bifurcation points. These are $x_1^\ast = \mp 0.959166$ and thus the sign is positive in both cases, implying that the Hopf bifurcations are subcritical.

This would explain all your observations, at least qualitatively. It would also predict that the limit cycles persist for $I < I_{crit,1}$ and $I > I_{crit,2}$, which can indeed be observed numerically.

Answer 3

My answer will be more intuitive, and largely a supplement to that of Hans Engler's.

Briefly: you are not seeing a supercritical Hopf Bifurcation but a subcritical Hopf bifurcation (as @Hans Engler pointed out in his answer), and the two recalcitrant facts you were trying to prove are only true for supercritical Hopf bifurcations, so that's why you were not observing them.

A few things to note going forward:

1. It is impossible to tell, from linearization results alone, what type of Hopf bifurcation you have. Both types will show the real parts of eigenvalues turning zero. Similarly, tracking your coefficient matrix through the trace-discriminant plane doesn't tell us which type of Bifurcation we are seeing either, only that one occurs.
2. An analytical criterion exists, which is the gold standard for determining which type of Hopf bifurcation you have. But it can be difficult to apply outside some tractable cases (see Hans Engler's answer). For more on this, see the excellent and daunting Exercise 8.2.12 of Strogatz (1994).
3. Thankfully, it is possible to use numerical techniques to get a sense (not a proof) for the type of bifurcation you have. Namely, if after you pass your critical value you get a small amplitude limit cycle, and it shrinks back down if you reverse the parameter, it is probably supercritical. Otherwise (if you see an irreversible jump) then it is probably subcritical.
4. Also, the behavior you observed (little stable spirals surrounded by stable limit cycles) are the hallmark of subcritical bifurcations: what you have there is a stable spiral surrounded by an unstable limit cycle surrounded by a stable limit cycle. This is pretty much the hallmarak of a subcritical Hopf bifurcation!

Points 3 and 4 are taken from Strogatz (1994) Section 8.2.

Indeed, the following diagram is from the chapter of the book after the one where you got the original system from. It shows exactly the type of behavior you observed (though for a different system of equations):

This is what your system is doing. All is well in the world.

Reference
Strogatz (1994) Nonlinear Dynamics and Chaos.