# 6. 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 neurodynex3.phase_plane_analysis import fitzhugh_nagumo
fitzhugh_nagumo.plot_flow()
fixed_point = fitzhugh_nagumo.get_fixed_point()
print("fixed_point: {}".format(fixed_point))
plt.figure()
trajectory = fitzhugh_nagumo.get_trajectory()
plt.plot(trajectory[0], trajectory[1])
```

## 6.1. Exercise: Phase plane analysis¶

We have implemented the following Fitzhugh-Nagumo model.

### 6.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.

### 6.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.

### 6.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
```

### 6.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?

## 6.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.

### 6.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

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
```

### 6.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`

.

### 6.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\).

## 6.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 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
list2.append(e[1].real)
# your code to plot list1 and list 2 against I goes here
```

### 6.3.1. Question¶

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

### 6.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?