4
$\begingroup$

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]:

bifurcation diagram

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:

trace-determinant diagram

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$:

pre-critical behavior

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:

oscillation analysis

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.

$\endgroup$
5
  • 2
    $\begingroup$ Are you sure this Hopf bifurcation is supercritical? (I haven't checked this myself.) What you are showing is consistent with a subcritical bifurcation of an unstable limit cycle that connects to a large amplitude stable limit cycle through a fold bifurcation at some smaller $I$. It's this large amplitude cycle that you would be observing here. To test this, look for a large amplitude cycle close to the one you are observing for $I < I_{crit}$, close to this value. $\endgroup$ Aug 9, 2020 at 21:25
  • $\begingroup$ @HansEngler following the textbook I'm using, I applied the Poincare-Bendixon theorem to show that there is a stable limit cycle at around $I=1.5$. But there is a wrinkle I left out -- if you continue the analysis there is another bifurcation where the limit cycle disappears around $I_{in}=2.03$. I left this out b/c I thought it was not relevant, but maybe I was wrong? But I have thought maybe I'm observing another limit cycle that I didn't know about, I even posted this as a comment then deleted though it certainly smells like a saddle node bifurcation. I'm a bit out of my depth. $\endgroup$
    – eric
    Aug 9, 2020 at 22:06
  • $\begingroup$ @HansEngler I should note the textbook says the system undergoes two Hopf Bifurcations (it leaves this as an exercise -- note this is not homework I am doing this for self study :) ). I have calculated two critical points: one at I = 0.97 and another at I=2.03. Folks can find the full pdf of the book here (this part is Chapter 8) if curious -- it is the author's web page so I think it is OK to post: cvr.yorku.ca/webpages/wilson.htm#book $\endgroup$
    – eric
    Aug 9, 2020 at 22:29
  • $\begingroup$ @HansEngler Yikes I will check over my calculations, will add edit to show what I did. $\endgroup$
    – eric
    Aug 10, 2020 at 0:31
  • $\begingroup$ @HansEngler I added my calculations to the end of my question, and also discuss what happens at I=0.866. Just to add to that: when you simulate above I=0.866, if that were the critical value we would expect it have an unstable eqm, but it doesn't I've done a lot of Jacobians above 0.866, below my calculated $Icrit$, and they always come back negative reals. Now this isn't to say I haven't found an interlocking system of mistakes all pointing to the same wrong answer, but I'd need to see the calculation to believe it at this point. :) $\endgroup$
    – eric
    Aug 10, 2020 at 3:04

3 Answers 3

2
+100
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Beautiful! Frankly this is far enough beyond my abilities right now that I will probably have to take a few hours tonight after work to study it. But it would make sense that I was flat-out wrong about the type of bifurcation, as I was going on intuition there not proof. But if this holds up my guess is this is the answer. $\endgroup$
    – eric
    Aug 10, 2020 at 17:22
  • $\begingroup$ I've been working through this and it seems promising. I think I got it to work if we change it so that $a=0.8$, $b=12.5$. Only then could I get the $\dot{x_2}$ half of things to replicate the original system, as $\dot{x_2}=10\dot{R}$. That would change the denominator b-a to 11.7, which means we have $-20+31$ so the sign is still positive, and the bifurcations are indeed subcritical. Does this seem right? $\endgroup$
    – eric
    Aug 11, 2020 at 4:54
  • $\begingroup$ Yes, it appears that this is a rescaling and the sign of the Lyapunov exponent would not be affected. $\endgroup$ Aug 19, 2020 at 21:44
2
$\begingroup$

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}]

enter image description here

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]

enter image description here

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]

enter image description here

$\endgroup$
3
  • $\begingroup$ Good to see independent confirmation of the critical points -- that was one discussion point above. $\endgroup$
    – eric
    Aug 10, 2020 at 13:46
  • $\begingroup$ Nice, I can certainly confirm this. Note also that these limit cycles persist for $I < I_{crit,1}$ and for $I > I_{crit,2}$, as an easy modification of this code shows. $\endgroup$ Aug 10, 2020 at 16:18
  • $\begingroup$ One thing this response is really elegant and makes me wish I used Mathematica and not just Python/Matlab. It really is elegant code and the plots are beautiful. $\endgroup$
    – eric
    Aug 10, 2020 at 17:23
1
$\begingroup$

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):

subcritical bifurcation

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

Reference
Strogatz (1994) Nonlinear Dynamics and Chaos.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .