CoCalc Public FilesLabs / Neuron / The FitHugh-Nagumo Model.ipynb
Views : 28
Compute Environment: Ubuntu 20.04 (Default)

This lab is based on a lab by Bard Ehrmentrout at the University of Pittsburgh (http://math.bu.edu/people/mak/MA666/11_FitzHugh_Nagumo.pdf).

Brains are composed of staggeringly complex networks of cells called neurons (pictured below). Neurons have two sides. The dendrites receive inputs from other neurons in the form of ions that alter the electrical background of the fluid the cell sits in. When this background passes a threshhold, the neuron fires, sending an electrical signal down the axon and triggering the axon terminals to emit ions that will trigger firing in the next neurons in the network.

This allows the network of cells to function like components of a computer that send zeros (not firing) and ones (fire!) throughout the network. Alan Hodgkin and Andrew Huxley first developed a mathematical model for the firing of a neuron, specifically the squid giant axon. They were awarded the Nobel Prize in Physiology or Medicine in 1963, the year after Watson and Crick.

The Hodgkin-Huxley model is a system of four differential equations. This is a bit complicated to visualize. In 1961 Richard FitzHugh proposed a simplified version of this model that had only two variables. The following year, Japanese scientist J. Nagumo built an electrical circuit with those exact dynamics (SageMath didn't exist yet, so they had no easy way to simulate solutions!)

# The Model

The FitzHugh-Nagumo model has two independent variables: the membrane potential $V$ and the recovery variable $W$. The goal was to capture the excitation and propagation phenomena from the chemical and electrical phenomena that triggered them. The membrane potential has a cubic nonlinearity that allows it to rapidly self-excite itself through positive feedback (as $V$ grows larger, $\frac{dV}{dt}$ grows larger still). It is limited only by the recovery variable that linearly provides negative feedback. The recovery variable grows when the membrane potential is high, allowing it to eventually ease the spike back to normal levels. The membrane potential is also subject to a constant forcing $I$ that can be understood as the inputs to the dendrites. Together, we get:

where $f(V)$ is a cubic function and $a$, $b$, and $c$ are positive parameters. For this lab, we will let $a = 0.08$, $b = 0.8$, $c = 0.7$.

1. When V is large and positive, what happens? When V is very negative, what happens? When V is close to zero, what happens? Write your answers in the Markdown cell below.

When $V$ is large and positive, $\frac{dV}{dt}$ will be very negative because you are subtracting a cubed term and $\frac{dW}{dt}$ is also going to be a bigger positive number. When $V$ is large and negative, $\frac{dV}{dt}$ will be very positive because you are adding a cubed term and $\frac{dW}{dt}$ is also negative. When $V$ is close to 0, $\frac{dV}{dt}$ can be approximated by the difference between $W$ and the constant forcing $I$ and $\frac{dW}{dt}$ will decrease to some fixed level because the $V$ term is not there to counteract the decay of $W$.

# Analysis

1. Let's begin by setting $I = 0$, indicating no external stimulus.

a. How many equilibrium points are there?

b. What is the stability of each?

c. Convince me of the number and stability of each equilibrium point in two ways.

a) There is 1 equilibria point at $(-1.199, -.624)$ b) From the Jacobian matrix, we found the trace $T = 0.628$ and the determinant $D = -0.0013$. Because $T>0$ and $D<0$, the equilibrium point is a saddle. This means the equilibrium point is unstable. c) We found the equilibrium point algebraically by setting the equations of the two nullclines together. We have also graphed the intersection of the nullclines, which shows that there is an equilibrium point at $(-1.199, -.624)$. We used trace-determinant analysis to determine that the equilibrium point is unstable, and you can also see this in the behavior of the phase portrait for the system.

In [36]:
# plotting the nullclines to see where they intersect
V,W,t = var('V,W,t')
p = plot((V + 0.7)/0.8,(V,-2,1), color = "red", thickness = 1, legend_label = "W-nullcline") #plots the W-nullcline
p += plot(V - (V^3)/3,(V,-2,1), color = "green", thickness = 1, legend_label = "V-nullcline") #plots the V-nullcline on the same graph
a = 0.08
b = 0.8
c = 0.7 #constructing our parameters
I = 0 #set I equal to 0

F = [V - (V^3)/3 - W + I, a*(V - b*W +c)]
normalizingFactor = sqrt(F[0]^2 + F[1]^2)
F_unit = [ F[0]/normalizingFactor, F[1]/normalizingFactor]
p += plot_vector_field(F_unit, (V,-2,1), (W,-2,2))

p


1. Find formulae for the V-nullclines and the W-nullclines when $I = 0$ and plot them together in different colors over the phase portrait. Place $V$ on the x-axis and W on the y-axis.

When I = 0, the V-nullcline is V - $\frac{V^3}{3}$ - W = 0 which gives W = V - $\frac{V^3}{3}$. The W-nullcline is W = $\frac{V + c}{b}$.

In [28]:
#plotting the nullclines over the phase portrait:
a = 0.08
b = 0.8
c = 0.7 #constructing our parameters
I = 0 #set I equal to 0

V, W,t = var('V,W,t')
F = [V - (V^3)/3 - W + I, a*(V - b*W +c)]
normalizingFactor = sqrt(F[0]^2 + F[1]^2)
F_unit = [ F[0]/normalizingFactor, F[1]/normalizingFactor]
p = plot_vector_field(F_unit, (V,-2,1), (W,-2,2))

p += plot((V + 0.7)/0.8,(V,-2,1), color = "red", thickness = 2, legend_label = "W-nullcline")
p += plot(V - (V^3)/3,(V,-2,1), color = "green", thickness = 2, legend_label = "V-nullcline") #plotting our respective nullclines
p


1. Use plots or mathematical demonstrations to describe what happens to the system as $I$ increases from zero. How does the phase portrait change? Are there any bifurcations? Do not find any solution trajectories (yet).
In [22]:
#using plots:

#plotting for I = 0:
a = 0.08
b = 0.8
c = 0.7 #constructing our parameters
I = 0 #set I equal to 0

V, W,t = var('V,W,t')
F = [V - (V^3)/3 - W + I, a*(V - b*W +c)]
normalizingFactor = sqrt(F[0]^2 + F[1]^2)
F_unit = [ F[0]/normalizingFactor, F[1]/normalizingFactor]
p = plot_vector_field(F_unit, (V,-2,1), (W,-2,2))
p += plot((V + 0.7)/0.8,(V,-2,1), color = "red", thickness = 2, legend_label = "W-nullcline")
p += plot(V - ((V^3)/3) + I,(V,-2,1), color = "green", thickness = 2, legend_label = "V-nullcline") #plotting our respective nullclines

p


In [23]:
#plotting for I not equal to 0:
a = 0.08
b = 0.8
c = 0.7 #constructing our parameters
I = 0.5 #set I equal to 0.5

V, W,t = var('V,W,t')
F = [V - (V^3)/3 - W + I, a*(V - b*W +c)]
normalizingFactor = sqrt(F[0]^2 + F[1]^2)
F_unit = [ F[0]/normalizingFactor, F[1]/normalizingFactor]
p = plot_vector_field(F_unit, (V,-2,1), (W,-2,2))
p += plot((V + 0.7)/0.8,(V,-2,1), color = "red", thickness = 2, legend_label = "W-nullcline")
p += plot(V - ((V^3)/3) + I,(V,-2,1), color = "green", thickness = 2, legend_label = "V-nullcline") #plotting our respective nullclines

p

In [2]:
#plotting for I not equal to 0:
a = 0.08
b = 0.8
c = 0.7 #constructing our parameters
I = 1 #set I equal to 1

V, W,t = var('V,W,t')
F = [V - (V^3)/3 - W + I, a*(V - b*W +c)]
normalizingFactor = sqrt(F[0]^2 + F[1]^2)
F_unit = [ F[0]/normalizingFactor, F[1]/normalizingFactor]
p = plot_vector_field(F_unit, (V,-2,1), (W,-2,2))
p += plot((V + 0.7)/0.8,(V,-2,1), color = "red", thickness = 2, legend_label = "W-nullcline")
p += plot(V - ((V^3)/3) + I,(V,-2,1), color = "green", thickness = 2, legend_label = "V-nullcline") #plotting our respective nullclines

p

In [33]:
#plotting for I not equal to 0:
a = 0.08
b = 0.8
c = 0.7 #constructing our parameters
I = 4 #set I equal to 4

V, W,t = var('V,W,t')
F = [V - (V^3)/3 - W + I, a*(V - b*W +c)]
normalizingFactor = sqrt(F[0]^2 + F[1]^2)
F_unit = [ F[0]/normalizingFactor, F[1]/normalizingFactor]
p = plot_vector_field(F_unit, (V,-2,1), (W,-2,5))
p += plot((V + 0.7)/0.8,(V,-2,1), color = "red", thickness = 2, legend_label = "W-nullcline")
p += plot(V - ((V^3)/3) + I,(V,-2,1), color = "green", thickness = 2, legend_label = "V-nullcline") #plotting our respective nullclines

p

In [35]:
V,W = var('V,W')
p = plot((V + 0.7)/0.8,(V,0,20), color = "red", thickness = 1, legend_label = "W-nullcline") #plots the W-nullcline
p += plot(V - (V^3)/3 +1,(V,0,20), color = "green", thickness = 1, legend_label = "V-nullcline") #plots the V-nullcline on the same graph

show(p)
#the nullclines will only ever intersect at one point, so there will only ever be one equilbrium point.


As $I$ increases, the phase portrait goes through several changes. In general, as $I$ incFor very low values of $I$, the equilibrium point seems to be a stable equilibrium. As $I$ gets bigger, the phase portrait shows that the solution becomes unstable because the arrows right near the equilibrium are pointing away from the equilibrium, not toward it. For very large values of $I$, the solution appears to turn back into a stable equilibrium. So, it appears that there is a Hopf bifurcation. The value of $I$ where the equilibrium goes from stable to unstable appears to be around $I = 0.4$ or $I = 0.5$.

the phase portrait stays pretty similar. The equilibrium point is shifted upward as a result of the increase in I, which causes a variation in the behavior of the solution depending on the initial condition, but the overall picture doesn't change dramatically. As I increases from zero we can see the phase portraits, the W-nullcline does not change. However, the y-intercept of the V-nullcline can change. There are two nullclines however because we have a system that goes from a stable system to an unstable system, then back to a stable one again.

1. Choose a fixed initial condition. For several choices of $I$, solve the system and plot the solution trajectory in the phase plane (alongside the nullclines). Also plot the $V(t)$ solution on its own axis. Be sure to demonstrate several different patterns of behavior such as no spike, one spike, and a "spike train."
In [69]:
#No spike

#Initial Conditions: V = 1, W = 1
a = 0.08
b = 0.8
c = 0.7 #constructing our parameters

#I = 0:
I = 0 #set I equal to 0

V, W,t = var('V,W,t')
F = [V - (V^3)/3 - W + I, a*(V - b*W +c)]
normalizingFactor = sqrt(F[0]^2 + F[1]^2)
F_unit = [ F[0]/normalizingFactor, F[1]/normalizingFactor]
p = plot_vector_field(F_unit, (V,-2,1), (W,-2,2))
p += plot((V + 0.7)/0.8,(V,-2,1), color = "red", thickness = 2, legend_label = "W-nullcline")
p += plot(V - ((V^3)/3) + I,(V,-2,1), color = "green", thickness = 2, legend_label = "V-nullcline") #plotting our respective nullclines

# plot the initial condition ()
solution = desolve_system_rk4(F, [V,W], ivar = t, ics = [0,1,1], end_points = 50, step = 0.1)
# produces a list of lists where solution[0] is [0,1,1] and solution[1] is [0.1, V(0.1), W(0.1)]
VSol = [ [time,VVal] for time,VVal,WVal in solution] # extract (t,V) as a list of ordered pairs
WSol = [ [time,WVal] for time,VVal,WVal in solution] # extract (t,W) as a list of ordered pairs

p2 = line(VSol, color = "green", legend_label = "V(t)") # plot (t,V)

trajectory = [ [VVal,WVal] for time,VVal,WVal in solution] # extract (V,W) as a list of ordered pairs

p += line(trajectory, color = "blue", thickness = 2) # add (V,W) to vector field
p += arrow(trajectory[int(len(trajectory)/5)] , trajectory[int(len(trajectory)/5)+1])

show(p)
show(p2)


In [51]:
#one spike

#Initial Conditions: V = 1, W = 1
a = 0.08
b = 0.8
c = 0.7 #constructing our parameters

#I = 2:
I = 2 #set I equal to 2

V, W,t = var('V,W,t')
F = [V - (V^3)/3 - W + I, a*(V - b*W +c)]
normalizingFactor = sqrt(F[0]^2 + F[1]^2)
F_unit = [ F[0]/normalizingFactor, F[1]/normalizingFactor]
p = plot_vector_field(F_unit, (V,-2,2), (W,-2,2))
p += plot((V + 0.7)/0.8,(V,-2,1), color = "red", thickness = 2, legend_label = "W-nullcline")
p += plot(V - ((V^3)/3) + I,(V,-2,1), color = "green", thickness = 2, legend_label = "V-nullcline") #plotting our respective nullclines

# plot the initial condition ()
solution = desolve_system_rk4(F, [V,W], ivar = t, ics = [0,1,1], end_points = 25, step = 0.1)
# produces a list of lists where solution[0] is [0,1,1] and solution[1] is [0.1, V(0.1), W(0.1)]
VSol = [ [time,VVal] for time,VVal,WVal in solution] # extract (t,V) as a list of ordered pairs
WSol = [ [time,WVal] for time,VVal,WVal in solution] # extract (t,W) as a list of ordered pairs

p2 = line(VSol, color = "green", legend_label = "V(t)") # plot (t,V)

trajectory = [ [VVal,WVal] for time,VVal,WVal in solution] # extract (V,W) as a list of ordered pairs

p += line(trajectory, color = "blue", thickness = 2) # add (V,W) to vector field
p += arrow(trajectory[int(len(trajectory)/5)] , trajectory[int(len(trajectory)/5)+1])

show(p)
show(p2)