4. 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

The phase_plane_analysis.fitzhugh_nagumo module contains all code required for this exercise. At the beginning of your exercise solutions, import the contained functions by running

from neurodynex.phase_plane_analysis.fitzhugh_nagumo import *

You can then simply run the exercise functions by executing the functions, e.g.


4.1. Exercise: Phase plane analysis

Create a script file (e.g. answers.py) and add the following header:

from neurodynex.phase_plane_analysis.fitznagumo import *

# your code here ..

You will type the code for your answers right below, adding more code with each exercise. When you want to execute your code, open ipython in a terminal and type

run answers.py

Or alternatively, from the command line execute

python answers.py

For some exercises you will have to plot and analyze data. This can be done by importing matplotlib and numpy:

import matplotlib.pyplot as plt
import numpy as np

4.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.

4.1.2. Question

Get the lists t, u and w by calling t, u, w = 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 current \(I\). Plot the nullclines for this given current and the trajectories into the \(u-w\) plane.

4.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

4.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?

4.2. Exercise: Jacobian & Eigenvalues

Consider the following two-dimensional 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}\]

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.

Set \(\varepsilon=.1\). Create the variable \(I\) and set it to zero for the moment.

4.2.1. Question

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

4.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.

4.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\).

4.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 use this example loop to help you getting started

list1 = []
list2 = []
currents = 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

4.3.1. Question

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

4.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?