# The Hodgkin-Huxley model¶

Here we examine the canonical Hodgkin-Huxley model (Hodgkin and Huxley, 1952). An uncertainty analysis of this model has been performed previously (Valderrama et al., 2015), and we here we repeat a part of that study using Uncertainpy.

The here used version of the Hodgkin-Huxley model has 11 parameters:

Parameter Value Unit Meaning
$$V_{0}$$ -10 mV Initial voltage
$$C_\mathrm{m}$$ 1 $$\text{F}/\text{cm}^2$$ Membrane capacitance
$$\bar{g}_{\mathrm{Na}}$$ 120 $$\text{mS/cm}^2$$ Sodium (Na) conductance
$$\bar{g}_{\mathrm{K}}$$ 36 $$\text{mS/cm}^2$$ Potassium (K) conductance
$$\bar{g}_{\mathrm{I}}$$ 0.3 $$\text{mS/cm}^2$$ Leak current conductance
$$E_\mathrm{Na}$$ 112 mV Sodium equilibrium potential
$$E_\mathrm{K}$$ -12 mV Potassium equilibrium potential
$$E_\mathrm{I}$$ 10.613 mV Leak current equilibrium potential
$$n_0$$ 0.0011   Initial potassium activation gating variable
$$m_0$$ 0.0003   Initial sodium activation gating variable
$$h_0$$ 0.9998   Initial sodium inactivation gating variable

As in the previous study, we assume each of these parameters have a uniform distribution in the range $$\pm 10\%$$ around their original value.

We use uncertainty quantification and sensitivity analysis to explore how this parameter uncertainty affect the model output, i.e., the action potential response of the neural membrane potential $$V_m$$ to an external current injection. The model was exposed to a continuous external stimulus of $$140 \mu \mathrm{A/cm}^2$$ starting at $$t = 0$$, and we examined the membrane potential in the time window between $$t$$ = 5 and 15 ms

As in the cooling coffee cup example, we implement the Hodgkin-Huxley model as a Python function (found in /examples/valderrama/valderrama.py):

import uncertainpy as un

import numpy as np
from scipy.integrate import odeint

# External stimulus
def I(time):
return 140 # micro A/cm**2

def valderrama(V_0=-10,
C_m=1,
gbar_Na=120,
gbar_K=36,
gbar_L=0.3,
E_Na=112,
E_K=-12,
E_l=10.613,
m_0=0.0011,
n_0=0.0003,
h_0=0.9998):

# Setup time
end_time = 15          # ms
dt = 0.025             # ms
time = np.arange(0, end_time + dt, dt)

# K channel
def alpha_n(V):
return 0.01*(10 - V)/(np.exp((10 - V)/10.) - 1)

def beta_n(V):
return 0.125*np.exp(-V/80.)

def n_f(n, V):
return alpha_n(V)*(1 - n) - beta_n(V)*n

def n_inf(V):
return alpha_n(V)/(alpha_n(V) + beta_n(V))

# Na channel (activating)
def alpha_m(V):
return 0.1*(25 - V)/(np.exp((25 - V)/10.) - 1)

def beta_m(V):
return 4*np.exp(-V/18.)

def m_f(m, V):
return alpha_m(V)*(1 - m) - beta_m(V)*m

def m_inf(V):
return alpha_m(V)/(alpha_m(V) + beta_m(V))

# Na channel (inactivating)
def alpha_h(V):
return 0.07*np.exp(-V/20.)

def beta_h(V):
return 1/(np.exp((30 - V)/10.) + 1)

def h_f(h, V):
return alpha_h(V)*(1 - h) - beta_h(V)*h

def h_inf(V):
return alpha_h(V)/(alpha_h(V) + beta_h(V))

def dXdt(X, t):
V, h, m, n = X

g_Na = gbar_Na*(m**3)*h
g_K = gbar_K*(n**4)
g_l = gbar_L

dmdt = m_f(m, V)
dhdt = h_f(h, V)
dndt = n_f(n, V)

dVdt = (I(t) - g_Na*(V - E_Na) - g_K*(V - E_K) - g_l*(V - E_l))/C_m

return [dVdt, dhdt, dmdt, dndt]

initial_conditions = [V_0, h_0, m_0, n_0]

X = odeint(dXdt, initial_conditions, time)
values = X[:, 0]

# Only return from 5 seconds onwards, as in the Valderrama paper
values = values[time > 5]
time = time[time > 5]

# Add info needed by certain spiking features and efel features
info = {"stimulus_start": time[0], "stimulus_end": time[-1]}

return time, values, info


We use this function when we perform the uncertainty quantification and sensitivity analysis (found in /examples/valderrama/uq_valderrama.py). We first initialize our model:

# Initialize the model
model = un.Model(run=valderrama,
labels=["Time (ms)", "Membrane potential (mV)"])


Then we create the set of parameters:

# Define a parameter dictionary
parameters = {"V_0": -10,
"C_m": 1,
"gbar_Na": 120,
"gbar_K": 36,
"gbar_L": 0.3,
"m_0": 0.0011,
"n_0": 0.0003,
"h_0": 0.9998,
"E_Na": 112,
"E_K": -12,
"E_l": 10.613}

# Create the parameters
parameters = un.Parameters(parameters)


We use set_all_distributions() and uniform() to give all parameters a uniform distribution in the range $$\pm 10\%$$ around their fixed value.

# Set all parameters to have a uniform distribution
# within a 20% interval around their fixed value
parameters.set_all_distributions(un.uniform(0.2))


set_all_distributions sets the distribution of all parameters. If it receives a function as input, it gives that function the fixed value of each parameter, and expects to receive Chaospy functions. uniform is a closure. It takes interval as input and returns a function which takes the fixed_value of each parameter as input and returns a Chaospy distribution with this interval around the fixed_value. Ultimately the distribution of each parameter is set to interval around their fixed_value:

cp.Uniform(fixed_value - abs(interval/2.*fixed_value),
fixed_value + abs(interval/2.*fixed_value)).


We can now use polynomial chaos expansions with point collocation to calculate the uncertainty and sensitivity of the model. We also set the seed to easier be able to reproduce the result.

# Perform the uncertainty quantification
UQ = un.UncertaintyQuantification(model,
parameters=parameters)
# We set the seed to easier be able to reproduce the result
data = UQ.quantify(seed=10)


The complete code for the uncertainty quantification and sensitivity becomes:

import uncertainpy as un
import chaospy as cp

from valderrama import valderrama

# Initialize the model
model = un.Model(run=valderrama,
labels=["Time (ms)", "Membrane potential (mV)"])

# Define a parameter dictionary
parameters = {"V_0": -10,
"C_m": 1,
"gbar_Na": 120,
"gbar_K": 36,
"gbar_L": 0.3,
"m_0": 0.0011,
"n_0": 0.0003,
"h_0": 0.9998,
"E_Na": 112,
"E_K": -12,
"E_l": 10.613}

# Create the parameters
parameters = un.Parameters(parameters)

# Set all parameters to have a uniform distribution
# within a 20% interval around their fixed value
parameters.set_all_distributions(un.uniform(0.2))

# Perform the uncertainty quantification
UQ = un.UncertaintyQuantification(model,
parameters=parameters)
# We set the seed to easier be able to reproduce the result
data = UQ.quantify(seed=10)