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.
get_trajectory()
get_fixed_point()
plot_flow()
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]\).
Note
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.
Note
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:
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
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.
Note
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\).
Note
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
list2.append(e[1].real)
# 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?