Hodgkin-Huxley

hh_aljoscha_cns.py

# Implementation of Hodgkin-Huxley model
# v2

import matplotlib.pyplot as plt

from numpy import *
from scipy import *

# Spike counting

def find_peaks(data, minimum = 0):
    """Takes numpy array and determines peaks based on a simple two-criterion
    heuristic of dA/dt = 0 and peak over minimum. Returns list of peak positions
    in data."""

    peaks = []

    # Ignore first and last element of array
    for i in xrange(1, len(data) - 1):
        point = data[i]
        if point > minimum:
            if (data[i-1] < point) and (data[i+1] < point):
                peaks.append(i)

    return peaks

def determine_firing_rate(trace, start = 0, timestep = 0.01):
    """Calculates approximate firing rate between specified start and
    trace end using non-weighted average."""

    # Remove data at beginning
    clean_trace = trace[(start/timestep):]

    # Gather peaks
    peak_table = find_peaks(clean_trace)

    distance_sum = 0

    for i in xrange(1, len(peak_table)):
        distance_sum = distance_sum + (peak_table[i] - peak_table[i-1])

    # Hack-ish. Assumes millisecond as unit for start and timestep!
    if distance_sum > 0:
        return 1000 / (timestep * (distance_sum / len(peak_table)))
    else:
        return 0

# HH parameter routines
def voltage_delta(v, params, c):
    """Vectorized function for dV/dt depending on m, n, h in params."""

    # Fixed settings for parameters
    e_l = -54.4
    e_k = -77
    e_na = 50
    
    g_l = 0.3
    g_k = 36
    g_na = 120

    # Calculate membrane current
    leak = g_l * (v - e_l)
    potassium = g_k * params[1]**4 * (v - e_k)
    sodium = g_na * params[0]**3 * params[2] * (v - e_na)

    membrane_current = leak + potassium + sodium

    # Calculate delta
    return ((-1 * membrane_current) + c)

def gating_delta_m(v, m):
    
    v_term = 0.1 * (v + 40)
    
    alpha = v_term / (1 - exp(-1 * v_term))
    beta = 4 * exp(-0.0556 * (v + 65))

    inf = alpha / (alpha + beta)
    tau = 1 / (alpha + beta)

    return (inf - m) / tau

def gating_delta_n(v, n):

    v_term = 0.01 * (v + 55)

    alpha = v_term / (1 - exp(-10 * v_term))
    beta = 0.125 * exp(-0.0125 * (v + 65))

    inf = alpha / (alpha + beta)
    tau = 1 / (alpha + beta)

    return (inf - n) / tau

def gating_delta_h(v, h):

    alpha = 0.07 * exp(-0.05 * (v + 65))
    beta = 1 / (1 + exp(-0.1 * (v + 35)))

    inf = alpha / (alpha + beta)
    tau = 1 / (alpha + beta)

    return (inf - h) / tau

# Simulation

def hh(currents = array([0]), length = 600, timestep = 0.01, offset = 50,
       v_init = -65, m_init = 0, n_init = 0, h_init = 1):
    """Simulates spiking neuron following Hodgkin-Huxley for a specified
    range of currents and modeling parameters. Returns a numpy array containing
    a trace for each input current."""

    steps = length / timestep
    off_currents = zeros_like(currents)
    
    # Initialize arrays for v and parameters
    v = zeros((steps, len(currents)))
    params = zeros((3, len(currents)))

    # Initialize voltages to init
    v[0] = v_init

    params[0] = m_init
    params[1] = n_init
    params[2] = h_init

    for i in xrange(0, int(steps) - 1):
        
        c = currents if i * timestep > offset else off_currents

        # Euler

        # Update parameters
        params[0] = params[0] + gating_delta_m(v[i], params[0]) * timestep
        params[1] = params[1] + gating_delta_n(v[i], params[1]) * timestep
        params[2] = params[2] + gating_delta_h(v[i], params[2]) * timestep

        # Update voltage
        v[i+1] = v[i] + voltage_delta(v[i], params, c) * timestep

    return v