Banana Peels and Action Potentials: a Detailed Look at the Mathematics of Action Potentials in the Fitzhugh-Nagumo System



This post has been edited to answer some questions asked on Facebook.

Professors Dr. Paul Carter of Arizona University and Dr. Bjorn Sandstede of Brown University wrote quite a large paper in the journal of SIAM  Journal on Applied Dynamical Systems  titled Unpeeling a Homoclinic Banana in the Fitzhugh-Namugo System. The paper was 114 pages long, and that doesn’t even include the references! What possibly could they have talked about for 114 pages? It must have been quite a complicated question for 114 pages!

Actually it is not really that hard of a question to ask. When considering the simplest mathematical models of an axon, why do some axons always fire two spikes instead of one? This seemingly simple questions took 114 pages of math to mathematically prove, and it was not as easy as unpeeling a banana as the title suggests. Mathematicians salivate over questions like these, simple to ask, difficult to prove. Obviously, I won’t be explaining in all of the gory mind-melting details, but I wanted to give you a flavor of why this question was even asked, what are the main ideas behind the proof, and what biologists can learn from such pure application of math.

As I said,  this is a mathematically simple model of an axon and uses the Fitzhugh-Nagumo model for its neural dynamics, and is a well-known relaxation oscillator model for action potentials.  It is a mathematical abstraction of the real biology, much like the Izhikevich model. This is perhaps the simplest model of an action potential we can write down mathematically. Here, the authors consider their axon in an one-dimensional domain along the axon. For a while now, mathematicians have noticed that for some parameters, this model fires two spikes, not one.  No matter how small the input into the axon, it always fires two spikes. This is a bit different than most biologists would expect. Normally, a small input should only give one action potential, not two. Some biologists more fluent in electrophysiology may point out that many neurons are capable of small phasic bursts, (a cluster of a few spikes), but those neuron models usually rely on slower ion channels (calcium activated potassium) or refractory periods (slow sodium inactivation) to produce phasic bursts.  These dynamics pieces are missing in this model. The short phasic burst of two spikes arises solely from the fact the axon has length and is a partial differential equation! This gives the normally two-dimensional Fitzhugh-Nagumo system an extra spatial dimension to work with. This extra room allows the model to fire an extra spike. Looking at the figure I’ve modified from their paper,  we see a stereotypical spike traveling down an axon.


Figure 1: This is a modified version of figure 2 from the paper. Here I have relabeled the axis to be more readable. Also, notice that after an final action potential their is a ringing or subthreshold oscillation. This turns out to be a key ingredient for the two action potentials in the Fitzhugh-Nagumo Oscillator. 


So now that we’ve poised the question let’s see how they set up the proof. They start by making the traditional change of coordinates. Here we care about traveling wave solutions, so they introduce a new  variable z= x +c t Where x is position on the axon, t is time, and c is the action potential speed along the axon. In this frame of reference, the system is now a three dimensional ordinary differential equation, and has a single equilibrium point. The equilibrium point is well characterized and is a special type of saddle point that has a homoclinic orbit. A homoclinic orbit is any saddle point such that when you leave the equilibrium point in the unstable direction you come back to the same equilibrium point along the stable direction. In this case this orbit represents the stable waveform as it travels down an axon. Figure 2, taken from Wikipedia, shows the stereotypical homoclinic orbit.  Now, normally these homoclinic orbits only go around once before returning back to rest. But not always. Because this is a three dimensional ODE (thanks to the spatial dimension), they can oscillate upon return to the equilibrium. This is known as ringing and can be seen in figure 1 (circled in yellow).  This turns out to be a major key to the proof. Only axons who have a ringing return to equilibrium can have two action potentials.

Figure 2: A Homoclinic Equilibrium point. Notice, how the unstable flow out of the equilibrium saddle returns back through the stable flow.  This cartoon is in the creative commons, and was made by Wikipedia user  mathmoclair.


However, there is a second idea behind this proof, the canard explosion. this can be a bit tricky to explain so let’s look at figure 3, which is another modified figure from the  paper. Normally, action potentials can be explained by these cubic shaped manifolds, which support relaxation oscillations. The key point here is there are two stable manifolds in solid red and one unstable manifold in dashed red. A manifold is an underlying surface (or this case line) where solutions tend to follow, or be repelled away from. In this case the homoclinic solution is attracted to the stable manifold branches and tends to follow them. The first action potential (blue) follows the traditional path, and as expected the homoclinic orbit does not like to follow  unstable manifold.  A canard explosion happens when the solution follows the unstable manifold for a while. This is of course unstable and eventually will fly to one of the stable branches. when it flies to the far stable branch this corresponds to firing the second action potential (dashed black). This only happens when the ringing is too large. It follows the unstable manifold for too long and fires a second action potential.

Figure 3: This is a modified version of figure 10 from Carter and Sandstede. Again, I have relabeled everything to be clearer. Solid red lines are the stable manifold, while the dashed line is the unstable manifold. The blue line is the first action potential (AP). The Black dashed line is the second AP. Notice how the second AP follows the unstable manifold in a canard explosion. Green is the ringing oscillation after the second action potential. The Airy point a point on the unstable manifold that determines canard expositions versus ringing oscillations.  


Thus the two main ingredients of the proof are to look at under what conditions does this happens. While conceptually this seems simple enough to do, mathematically this is no easy matter, and requires a lot of singular perturbation theory in order to accomplish. The trick here is actually tracking how close the solutions get to the manifolds. To do this often you consider sub regions and track when a solution enters and exits. Take example of the author’s figure 13, as an example. I have modified it a bit to make my particular point clearer.  Here they zoom in on the part of the cubic manifold with the equilibrium point on it. The red solid line stable manifold. The red dashed line is the unstable manifold. The blue line is the path of the beginning of the second action potential. One will notice the blue path exits on the right.  This is the canard explosion. The green line the ringing, after the second Action potential. Notice that these fail to exit from the right. They reenter across the blue dashed line, and decay after each cycle. Thus they never explode to a second action potential. Solving this mathematically usually requires solving a lot of boundary layer value problems.  Each time the path enters a region they track the solution, estimating the highly nonlinear equations sufficiently enough to know approximately where they will exit. Furthermore, they track them just enough to know if the exit/entrance point is less than or greater than the previous exit/entrance point.  At the end this gives them the proof they are looking for, two action potentials or one.
Figure 4 is a modified version of figure 13 from Carter and Sandstede 2018. I have relabeled many of the lines so that each line has a name rather than a symbol. As we can see depending on where a solution enters the purple region, also determines where it exits. The whole crux of the proof is tracking where these lines enter and exit, and this takes 114 pages of math to prove. 



Obliviously there is a lot going on, and above is just one of the many pieces of many layers they have to peel away to solve the problem. But the mathematically inclined biologists can learn a few things from this study. Perhaps the easiest, and most philosophical point is there are a lot of mechanisms that can cause a phenomenon in biology, and just because one mechanism can cause something doesn’t mean that a completely different mechanism is mathematically possible. As an example most electrophysiologists know most phasic bursts are caused by calcium activated potassium channels. Here we see it could be caused, or at least reinforced by the physical existence of the axon as well.

The second point, and more mathematically difficult, is that having a good geometrical understanding of the dynamics of the system will help computational biologists find different phenomena in their model. This proof was first conceived as geometric proof, thinking about how the homoclinic orbit was peeled by the unstable manifold. The math just glues the geometric jigsaw puzzle together.

The third lesson is singular perturbation is a very powerful tool. It is impossible to solve these equations analytically. You need to estimate the solution so that its analysis is simpler. These ideas don’t just hang out in the depths of the applied math department. They are making headways into actual experiments, and having at least a passing understanding of the main results can help motivate new and exciting experiments.




**Note, I have modified each figure from the paper to make things a bit clearer to those who haven’t read the paper. All credit for the simulations and figures go to Dr. Carter, and Dr. Sandstede. I only relabeled things using the power of Photoshop.


------------------------------------------------------------------------------------------------------

Over the past few weeks there have been several questions people have asked me over facebook, or in person. I would like to address them below. The questions are great questions, and I feel that we all can learn from them. 



What is the injected current protocol used in the paper?

The professors used a mixture of delta impulses, and choosing different initial conditions. It depends on which figure you’re referring to. For my modified figure 1 they used delta impulses. However, for some figures, modified figures 3 and 4, they just choses different initial conditions. When they just chose initial conditions they use initial conditions near the equilibrium, but not on it. Because the equilibrium point is a homoclinic saddle, the only way for it to return to rest is to spike (once or twice) then return to rest. For biologists this may seem very weird, but this is mathematically equivalent to a delta impulse, which is an infinitely strong and infinitesimally short pulse of current. Mathematically delta impulses are the same as resetting the initial conditions, which is what the professors did interchangeably. As for experimental biologists, they can think of delta impulses as a mathematical idealization of brief current pulses of about 1 or 2 milliseconds of length. It is just enough current to get the axon to spike, then turns off as soon as the first spike occurs.

Does the Hodgkin-Huxley model have a ringing phenomenon?



Yes, the Hodgkin-Huxley model does have this ringing effect but it decays much faster than the ringing in the Fitzhugh-Nagumo model. I’ve attached a picture of this ringing below. The below picture was generated using the model parameters for Dr. Ermentrout's  book. These are the exact same parameters used in the original Hodgkin-Huxley model. I used initial conditions to insure there would be a single spike, with V= -62, n =0, m = 0, h = 0. I set the background current to 2 nA. This isn’t necessary, it just makes the ringing easier to see. The ring still exists with the background current set to 0 nA. This after spike ringing is known as subthreshold oscillations, and is also known as resonation. Many neurons are indeed resonators, but many more neurons are not resonators. Whether or not a neuron is a resonator is dependent on whether the rest equilibrium point has complex eigenvalues. You can see this blog article for a more in-depth discussion of resonance. As a final point the Fitzhugh model is sometimes a resonator, and sometimes it is not a resonator. It depends on the parameters. Here the professors focus on parameter sets where the Fitzhugh-Nagumo model is a resonator, as it is the only set of parameters capable of generating two spikes.


Figure 5: The ringing in the Hodgkin-Huxley model. Notice that the amplitude of the subthreshold oscillations are smaller than the Fitzhugh Nagumo model. Top shows the full spike, while bottom zooms in on the subthreshold range between -50 mV to -70 mV.


Is the ringing phenomenon, and the one spike, two spike a side-effect of the Fitzhugh-Nagumo models reduction?

The answer here is both yes and no, depending on exactly how you define which neuron your trying to fit. Historically the Fitzhugh-Nagumo model was a mathematically reduced version of the Hodgkin-Huxley model of squid neurons. Later mathematicians changed the parameters of the Fitzhugh-Nagumo model to investigate all the different behaviors this simple model could exhibit. When doing this they found this one-spike, two-spike phenomenon that the authors study. This model is no longer representing the Hodgkin-Huxley model, so in this narrow sense the answer is, “Yes, it is a bad fit to squid neurons.” However, that is not the same as saying no real biological neuron could never behave like this. As an example, I can take this particular parameter set that the authors have and construct a neural model using the Hodgkin-Huxley formalism. The Hodgkin-Huxley formalism is a specific way to describe ion channels kinetic and is the gold standard when building biologically realistic neuron models. I am confident that someone can easily build a model using the HH-formalism to build a neuron that looks behaves the same as the models the authors study. This model could be biologically plausible, and could very well exist in nature. It is just a matter of finding the right set of parameters. So in the broader sense the answer is, “No, this is not a consequence of bad neural fitting. This phenomenon is mathematically possible in real biological systems.”

Is the phenomena found in real biological systems?

My guess is we will never really know. As I said in the blog article, slow calcium activated potassium channels are a much more reliable way to produce these short burst of a few action potentials. My educated guess is, if biology does use the phenomenon, it uses it at the same time as other different mechanisms, such as slow calcium activated potassium channels. Often Biology uses the many different mechanisms to do the same function. This makes biology more robust, and less likely to fail. It would also make any experiment to find this phenomenon very difficult to perform.

留言