7. FitzHugh-Nagumo: Phase plane and bifurcation analysis

Book chapters

See Chapter 4 and especially Chapter 4 Section 3 for background knowledge on phase plane analysis.

Python classes

In this exercise we study the phase plane of a two dimensional dynamical system implemented in the module phase_plane_analysis.fitzhugh_nagumo. To get started, copy the following code block into your Jupyter Notebook. Check the documentation to learn how to use these functions. Make sure you understand the parameters the functions take.

%matplotlib inline
import brian2 as b2
import matplotlib.pyplot as plt
import numpy as np
from neurodynex.phase_plane_analysis import fitzhugh_nagumo


fixed_point = fitzhugh_nagumo.get_fixed_point()
print("fixed_point: {}".format(fixed_point))

trajectory = fitzhugh_nagumo.get_trajectory()
plt.plot(trajectory[0], trajectory[1])

7.1. Exercise: Phase plane analysis

We have implemented the following Fitzhugh-Nagumo model.

(1)\[\begin{split}\left[\begin{array}{ccll} {\displaystyle \frac{du}{dt}} &=& u\left(1-u^{2}\right)-w+I \equiv F(u,w)\\[.2cm] {\displaystyle \frac{dw}{dt}} &=& \varepsilon \left(u -0.5w+1\right) \equiv \varepsilon G(u,w)\, ,\\ \end{array}\right.\end{split}\]

7.1.1. Question

Use the function plt.plot to plot the two nullclines of the Fitzhugh-Nagumo system given in Eq. (1) for \(I = 0\) and \(\varepsilon=0.1\).

Plot the nullclines in the \(u-w\) plane, for voltages in the region \(u\in\left[-2.5,2.5\right]\).


For instance the following example shows plotting the function \(y(x) = -\frac{x^2}{2} + x + 1\):

x = np.arange(-2.5, 2.51, .1)  # create an array of x values
y = -x**2 / 2. + x + 1  # calculate the function values for the given x values
plt.plot(x, y, color='black')  # plot y as a function of x
plt.xlim(-2.5, 2.5)  # constrain the x limits of the plot

You can use similar code to plot the nullclines, inserting the appropriate equations.

7.1.2. Question

Get the lists t, u and w by calling t, u, w = fitzhugh_nagumo.get_trajectory(u_0, w_0, I) for \(u_0 = 0\), \(w_0= 0\) and \(I = 1.3\). They are corresponding values of \(t\), \(u(t)\) and \(w(t)\) during trajectories starting at the given point \((u_0,w_0)\) for a given constant input current \(I\). Plot the nullclines for this given current and the trajectories into the \(u-w\) plane.

7.1.3. Question

At this point for the same current \(I\), call the function plot_flow, which adds the flow created by the system Eq. (1) to your plot. This indicates the direction that trajectories will take.


If everything went right so far, the trajectories should follow the flow. First, create a new figure by calling plt.figure() and then plot the \(u\) data points from the trajectory obtained in the previous exercise  on the ordinate.

You can do this by using the plt.plot function and passing only the array of \(u\) data points:

u = [1,2,3,4]  # example data points of the u trajectory
plot(u, color='blue')  # plot will assume that u is the ordinate data

7.1.4. Question

Finally, change the input current in your python file to other values \(I>0\) and reload it. You might have to first define \(I\) as a variable and then use this variable in all following commands if you did not do so already. At which value of \(I\) do you observe the change in stability of the system?

7.2. Exercise: Jacobian & Eigenvalues

The linear stability of a system of differential equations can be evaluated by calculating the eigenvalues of the system’s Jacobian at the fixed points. In the following we will graphically explore the linear stability of the fixed point of the system Eq. (1). We will find that the linear stability changes as the input current crosses a critical value.

7.2.1. Question

Set \(\varepsilon=.1\) and \(I\) to zero for the moment. Then, the Jacobian of Eq. (1) as a function of the fixed point is given by

\[\begin{split}\begin{aligned} J\left(u_{0},w_{0}\right) & = & \left.\left(\begin{array}{cc} 1-3u_0^2 & -1\\[5pt] 0.1 & -0.05 \end{array}\right)\right.\end{aligned}\end{split}\]

Write a python function get_jacobian(u_0,w_0) that returns the Jacobian evaluated for a given fixed point \((u_0,v_0)\) as a python list.


An example for a function that returns a list corresponding to the matrix \(M(a,b)=\left(\begin{array}{cc} a & 1\\ 0 & b \end{array}\right)\) is:

def get_M(a,b):
        return [[a,1],[0,b]] # return the matrix

7.2.2. Question

The function u0,w0 = get_fixed_point(I) gives you the numerical coordinates of the fixed point for a given current \(I\). Use the function you created in the previous exercise  to evaluate the Jacobian at this fixed point and store it in a new variable J.

7.2.3. Question

Calculate the eigenvalues of the Jacobian J, which you computed in the previous exercise , by using the function np.linalg.eigvals(J). Both should be negative for \(I=0\).

7.3. Exercise: Bifurcation analysis

Wrap the code you wrote so far by a loop, to calculate the eigenvalues for increasing values of \(I\). Store the changing values of each eigenvalue in seperate lists, and finally plot their real values against \(I\).


You can start from this example loop:
import numpy as np
list1 = []
list2 = []
currents = np.arange(0,4,.1) # the I values to use
for I in currents:
    # your code to calculate the eigenvalues e = [e1,e2] for a given I goes here
    list1.append(e[0].real) # store each value in a separate list

# your code to plot list1 and list 2 against I goes here

7.3.1. Question

In what range of \(I\) are the real parts of eigenvalues positive?

7.3.2. Question

Compare this with your earlier result for the critical \(I\). What does this imply for the stability of the fixed point? What has become stable in this system instead of the fixed point?