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)