The dynamical systems view of the Hodgkin-Huxley model

hodgkin_huxley_animated.py

# -*- coding: utf-8 -*-
"""
Voltage spikes, also known as action potentials, are considered the key unit of information flow in
the nervous system. One of the first quantitative descriptions of the action potential is the
Hodgkin-Huxley Model, named after nobel prize winning physiologists Alan Lloyd Hodgkin and Andrew
Fielding Huxley. This matplotlib animation shows the result of simulating the Hodgkin-Huxley model
in python with parameters that are set in a way, that the model resembles a cortical pyramidal cell.
At 50ms a current is injected, making the model spike. The upper panel shows the dynamical systems
view of the action potential. Before the current injection, the model stays at the resting
potential. When the current is injected a bifurcation occurs. The system finds a limit cycle.
This limit cycle is formed by an interaction of the current injection and the potassium activation
variable n. The current injection brings the voltage towards positive values. The potassium
activation variable on the other hand, increases with the voltage but drives the voltage towards
the initial resting potential. The lower panel shows the same voltage over time. At 250ms the
current injection stops and the voltage returns to resting equilibrium. Animations like this,
created with matplotlib, capture the aspect of time in a way that static plots can’t and are
therefore great a taking a more dynamic look at the action potential

@scidanm
"""

import numpy as np
from matplotlib.animation import FuncAnimation
import matplotlib.pyplot as plt

# Load the data
data = np.loadtxt("corticalpc_data.txt")
results_v = data[0]
results_u = data[1]
t = data[2]
dt = 0.1

# Animate the trajectory
fig = plt.figure()
fig.set_size_inches(16,9)
fig.tight_layout()
# Making the figure fullscreen might fail depending on the backend

ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
#line, = ax1.plot([], [], lw=0.5)

xticklabels = ax1.get_xticklabels()
for x in xticklabels:
    x.set_size(12)
yticklabels = ax1.get_yticklabels()
for y in yticklabels:
    y.set_size(12)
xlabel1 = ax1.xaxis.get_label()
ylabel1 = ax1.yaxis.get_label()
xlabel1.set_text('potassium activation n')
xlabel1.set_size(12)
ylabel1.set_text('voltage (mV)')
ylabel1.set_size(12)
xlabel2 = ax2.xaxis.get_label()
ylabel2 = ax2.yaxis.get_label()
xlabel2.set_text('time (ms)')
xlabel2.set_size(12)
ylabel2.set_text('voltage (mV)')
ylabel2.set_size(12)

fig.show()
ax2.set_xlim((t[0], t[-1]))
ax2.set_ylim((results_v.min()-5, results_v.max()+5))

def update_plot(i):
    ax1.clear()
    if i<100:
        ax1.plot(results_u[0:i], results_v[0:i], color='k', alpha=0.9, lw=1)
    else:
        ax1.plot(results_u[0:i-100], results_v[0:i-100], color='r', alpha=0.5, lw=1)
        ax1.plot(results_u[0+i-100:i], results_v[0+i-100:i], color='k', alpha=0.9, lw=1)
    ax2.plot(t[0:i], results_v[0:i], color='k', lw=1)
    ax1.set_xlim(results_u.min()-0.01, results_u.max()+0.01)
    ax1.set_ylim(results_v.min()-2, results_v.max()+2)
    txt = ax1.text(0.089,13, "t=" + "{0:.1f}ms".format(i * dt), size=12)
    xticklabels = ax1.get_xticklabels()
    for x in xticklabels:
        x.set_size(12)
    yticklabels = ax1.get_yticklabels()
    for y in yticklabels:
        y.set_size(12)
    xlabel = ax1.xaxis.get_label()
    ylabel = ax1.yaxis.get_label()
    xlabel.set_text('potassium activation n')
    xlabel.set_size(12)
    ylabel.set_text('voltage (mV)')
    ylabel.set_size(12)
    fig.show()

anim = FuncAnimation(fig, update_plot, frames=np.arange(0, len(t),8), interval=1)
anim.save('hh-dynamics.mp4', dpi=150, fps = 30, writer='ffmpeg')